Stability analysis on dark solitons in quasi-1D Bose–Einstein condensate with three-body interactions

The stability properties of dark solitons in quasi-one-dimensional Bose–Einstein condensate (BEC) loaded in a Jacobian elliptic sine potential with three-body interactions are investigated theoretically. The solitons are obtained by the Newton-Conjugate Gradient method. A stationary cubic-quintic nonlinear Schrödinger equation is derived to describe the profiles of solitons via the multi-scale technique. It is found that the three-body interaction has distinct effect on the stability properties of solitons. Especially, such a nonlinear system supports the so-called dark solitons (kink or bubble), which can be excited not only in the gap, but also in the band. The bubbles are always linearly and dynamically unstable, and they cannot be excited if the three-body interaction is absent. Both stable and unstable kinks, depending on the physical parameters, can be excited in the BEC system.


Model
At ultra-low temperatures, the dynamical behaviors of BECs with two-and three-body interactions can be described by the following three-dimensional (3D) nonlinear GPE 11 with two-body and three-body interaction where � = �(r, t) labels the condensate wave function, V (r) is an experimentally generated macroscopic potential, r = (x, y, z) the Cartesian coordinate vector, ∇ 2 the Laplacian operator, the reduced Planck constant, m the mass of the atom. g 1 = 4π 2 a s m denotes the strength of two-body interaction. a s is the s-wave scattering length ( a s > 0 and a s < 0 respectively represents the repulsive and attractive interaction), which can be tuned to any desired value by using the "Feshbach resonance" technique 42 . g 2 is the strength of the three-body interaction. In Ref. 43 , g 2 is given by a universal formula g 2 = 12π 2 a 2 s m d 1 + d 2 tan s 0 ln |a s | |a 0 | + π 2 , where the numerical values of the universal constants d 1 , d 2 , a 0 and s 0 are given in Refs. 43,44 . From Ref. 43 , we know that g 2 can be tuned from −∞ to +∞ . The total atom particles are N = | | 2 d 3 r.
In experiments, the BEC atoms are usually confined in a harmonic potential V (r) = 1 2 m(ω 2 x x 2 + ω 2 y y 2 + ω 2 z z 2 ) with ω x , ω y and ω z being the trap frequencies along x, y and z-directions. In the disk-shaped condensates, i.e., ω x ≈ ω y and ω z ≫ ω x , the 3D GPE can be reduced to 2D GPE. In the cigar-shaped condensates, i.e., ω y , ω z ≫ ω x , the 3D GPE can be reduced to 1D GPE. In this paper, we only consider the cigar-shaped condensate. By introducing the dimensionless variables t = ω x t , x = x/a h0 , � = a 3 h0 /n 0 � , where a h0 = √ /mω x is the characterized length of harmonic oscillator, n 0 is a given particle density. In general, one can take n 0 = N (This is not necessary.). By using the approach presented in Ref. 45 and omitting the tilde ′ ∼ ′ above all the variables, one can obtain the following quasi-1D dimensionless GPE Under which, the total particle number can be expressed as N = n 0 |�(x, t)| 2 dx . We choose the periodic external potential as V (x) = V 0 sn 2 (x, q) , where sn(x, q) is the Jacobian elliptic sine function with modulus q(0 ≤ q < 1) . This potential can be regarded as a generalization to the trigonometry function. In experiments, such a potential can be well approximated by using only two laser beams 46,47 .
Solitary waves of Eq. (2) are sought in the form where µ is the chemical potential, ψ(x) is a real-valued function, which satisfies the equation is an exact solution of Eq. (4), then so does ψ = −ϕ(x) . When ψ(x) is infinitesimal, the terms ψ 3 and ψ 5 in Eq. (4) can be neglected, which results in a linear equation The bounded solutions of Eq. (5) are called linear Bloch modes, and the corresponding constant µ forms linear Bloch bands. Since V(x) is a periodic function, Eq. (5) is a generalized form of Mathieu's equation. Its bounded solution can be written as 48 where p(x; µ) has the same period as the potential V(x), µ = µ(k) is the 1D dispersion relation. Both the linear Bloch bands and the dispersion relations can be gotten by solving the obtained eigenvalue problem. However, unfortunately, it is rather difficult to solve the eigenvalue problem exactly and analytically. Here, we solve it numerically by the Fourier collocation method 49 . The graph of Bloch band structures is omitted for it is similar as that shown in Ref. 46 .

Solitons and their stability properties
When ψ(x) is not infinitesimal, the terms ψ 3 and ψ 5 in Eq. (4) can not be neglected. Considering that the chemical potential µ enters into the band-gap, where the linear Bloch waves no longer exist, from the boundary of the band k = k 0 , µ = µ(k 0 ) , completely localized solitary waves, namely gap solitons 32 , can be excited. The gap solitons can be found numerically by the NCG method 28 , which is an effective numerical method for seeking the solitary wave solutions of nonlinear evolution equations. It is based on Newton iteration and conjugate-gradient iteration method to solve the resulting linear equation. This method can be applied to compute both the ground states and excited states in various physical systems. More importantly, this method usually converges faster than the other numerical methods and it is easy to implement in Matlab. A detailed description about the NCG method can be found in Ref. 28 .

Bright solitons.
In literatures, the bright solitons are often defined as the solutions with vanishing boundary conditions, i.e. ψ(±∞) = 0 . As we had done in Ref. 46 , we numerically find that virous of bright solitons (for examples, the on-site and off-site gap solitons) still exist when the three-body interaction is taken into account.
Here we omitted the profiles of the bright solitons as they have the similar structures as those illustrated in Ref. 46 . We also noted that the amplitude of the bright solitons decreases when the three-body interaction strength |g 2 | increases.
To see it more clearly, we define the amplitude of bright soliton as A = max(|ψ|) and the particle number P = |ψ| 2 dx , which implies that the total particles is N = n 0 P . In nonlinear optics, P is often called the power. Figure 1 exhibits the amplitude and power of on-site and off-site solitons versus the chemical potential µ for different three-body interaction strength g 2 . As can be seen from Fig. 1a and b, both A and P decrease when µ moves toward the first band for fixed g 2 . In the semi-infinite gap, P linearly decreases when µ increases but far away from the band edge. However, P decreases rapidly if µ moves near the band. This is opposite in the first gap, i.e., P decreases linearly when µ decreases but far away from the first band edge. We also noticed from Fig. 1 that the larger the three-body interaction |g 2 | , the smaller the amplitude A and the power P. This conclusion will be explained theoretically later. Figure 2 shows the amplitude of on-site gap solitons versus the three-body interaction strength g 2 under different two-body interaction strength g 1 . Other parameters are V 0 = 2 and q = 0.1 . For a fixed g 2 , it is obvious that the larger the |g 1 | is, the smaller the amplitude of the solitons. In the semi-infinite gap, see Fig. 2a, the amplitude of solitons increases as g 2 increases for a given g 1 . However, in the first gap, this is somewhat opposite. That is, see Fig. 2b, the amplitude of solitons decreases when g 2 increases. On the other hand, when the two-body interaction is stronger, i.e. |g 1 | is relatively larger, one can see from Fig. 2 that the amplitude of gap solitons nearly invariant with the increasing of g 2 . This is because the amplitude of gap solitons is very small when |g 1 | is larger, so that the term ψ 5 in Eq. (5) can be neglected. It is also noted that gap solitons indeed exist when g 1 g 2 < 0 for a relatively narrow interval of g 2 . We will make a theoretical explanation later.
It is more difficult to find dark solitons than bright ones when the NCG method is applied. One must choose the initial guess carefully, or the iteration will be divergent or convergent to an unwanted result. In practice, we take the initial guess www.nature.com/scientificreports/ for seeking the kink and for the bubble, where T x is the periodicity of the external potential ( T x ≈ π when q = 0.1 ). One can choose appropriate values for a 1 , a 2 and a 3 to obtain the wanted results. Of course, the coefficient in front of the cosine function also can be adjusted if necessary. Figure 3 shows the profiles of kinks ( Fig. 3a and c) and bubbles ( Fig. 3d and f) given by the NCG method. The amplitudes of the dark solitons are defined as illustrated in Fig. 3, where the upper green dashed lines denotes the average value of the Bloch waves in one periodicity. One can see that the kink is odd symmetric, while the   www.nature.com/scientificreports/ bubble is even symmetric. When |x| is large enough, the matter waves oscillate with the same periodicity as the external potential. This is because of the periodical carrier Bloch wave has been modulated by a kink or a bell-like soliton. Figure 3b and e show the residual error, measured as the maximum of the residue in Eq. (4), when the NCG method is applied to seek the dark solitons shown in Fig. 3a and d, respectively. It is seen that the residual drops below 10 −12 with less than or about 200 times iterations, implying that the NCG method can quickly gain the dark solitons with very high accuracy. When the numerical calculation is performing, the values of a 1 , a 2 and a 3 for initial guesses are taken as a 1 = 1.8 for Fig. 3a, a 1 = 0.45 for Fig. 3c, a 2 = −a 3 = −0.25 for Fig. 3d and a 2 = −a 3 = −0.45 for Fig. 3f. The solitons given in Fig. 3 can be used as the initial guesses for NCG method to seek the dark solitons for other parameters. It is also can be seen easily from Fig. 3 that the amplitude of dark solitons decrease when three-body interaction strength g 2 increases. Figure 4a and b exhibit the amplitude of the kinks and bubbles versus the chemical potential µ under different nonlinear strength g 1 < 0 and g 2 > 0 . From Fig. 4, we see that the dark solitons can be excited not only in the gap, but also in the band. This is quite different to the bright soliton discussed previously. The bright solitons exist in the gaps. On the other hand, in the semi-infinite gap, the dark solitons can be exited only for µ is greater than a certain value. Take an example, when g 1 = −1, g 2 = 0.5 , the kinks exist for µ > 0.355 and the bubbles for µ > 0.532 , while the bright solitons exist for the entire gap. Whenµ moves toward the band edge, the amplitude of the bright soliton becomes smaller and smaller. However, this is not the fact for the dark solitons. From Fig. 4a and b, we see that the amplitude of the dark solitons increases monotonously as µ increases. It is known from "Bright solitons" section that the amplitude of the bright soliton decreases when the two-body interaction strength |g 1 | increases (other parameters are fixed). To our surprise, the amplitude of the dark solitons increases as |g 1 | increases.  www.nature.com/scientificreports/ Figure 5a and b display the amplitude of the kinks and bubbles versus the chemical potential µ under different nonlinear strength g 1 > 0 and g 2 < 0 . Here, we start computing from the points P 1 (corresponding to µ = 1 ). The parameters for initial guess are a 1 = 0.45, a 2 = −a 3 = −0.25 . The corresponding wave profiles at the marked points P 1 − P 4 are shown in the subplots. The dashed lines denote the unstable solitons while the solid line denotes stable ones (The stable property will be discussed later.). Again, we see that the dark solitons can be excited in the band. It can be seen from Fig. 5a that the kink solitons exist in the first band and first gap, but it can not be excited in the semi-infinite gap. The amplitude of the kink increases as chemical potential µ increases. When µ moves from P 1 to P 2 , the amplitude becomes smaller distinctly. Both the amplitude of kink (modulation wave) and the amplitude of the Bloch wave (carrier wave) tend to zero as µ moves toward the lower edge of the first band. The kink soliton and the Bloch wave no longer exist when µ falls into the semi-infinite gap. From Fig. 5a, we suspect that the kink soliton may be regarded as a kink function, such as the tanh function, multiply by a Bloch wave function when the chemical potential µ near the lower edge of the first band.
When it comes to the bubbles (see Fig. 5b), the numerical results are quite different to the kink ones. One can see that the amplitude of the bubbles decreases as µ increases. When µ falls into the semi-infinite gap, the gray soliton reduces to an on-site-like bright soliton (see the wave profile at point P 3 ). This is because the linearly periodic Bloch wave is prohibited in the semi-infinite gap. The profiles of this kind of solitons are similar to the bright solitons discussed before. However, they have different dynamical behaviors. The numerical results indicate that the on-site-like solitons are linearly unstable and dynamically unstable if µ is near the band edge. There is a critical value for µ . The on-site-like solitons are linearly and dynamically stable only when µ is less than this critical value. For the parameters used in Fig. 5b ( g 1 = 1, g 2 = −0.5 ), the critical value for µ is about 0.438. This critical value depends on the nonlinear strength g 1 and g 2 . From Fig. 5b, we suspect that the bubble soliton may be regarded as a superposition of a completely localized function and a Bloch wave function when µ near the lower edge of the first band.
Multi-scale method. We now make a theoretical analysis on the solitons obtained in the above section.
Considering that µ falls into the band-gap from the band edge k = k 0 , µ 0 = µ(k 0 ) , then ψ(x) and µ can be expanded with multi-scale X 0 = x , X 1 = ε 1/2 x , that is where ε = k − k 0 is a small quantity. Suppose that g 1 = O(ε), g 2 = O(ε) . Substituting the above expansions into Eq. (4), one can get (9) is similar as Eq. (5), thus it possesses solution ψ 0 = B(X 1 )p(X 0 ) , where p(X 0 ) is the carrier wave and B(X 1 ) is the modulating wave. From Eq. (10), we then have ψ 1 = dB dX 1 H(X 0 ) , where H(X 0 ) is a periodic function, which satisfies L 0 H = dp dX 0 . It's easy to verify that the Fredholm condition is satisfied automatically. Substituting ψ 0 and ψ 1 into Eq. (11) yields Applying the Fredholm condition to Eq. (12) leads to the following stationary nonlinear Schrödinger equation for B(X 1 ) 49 . When the three-body interaction can be neglected, i.e. g 2 = 0 , then Eq. (13) reduces to the stationary cubic NLSE. Note that α 1 and α 2 are always positive.
Equation (13) has a localized solitary wave solution, which reads . Solution (14) is just the so-called "bright soliton". If the three-body interaction is negligible, i.e. g 2 = 0 , this soliton reduces to the bell-like soliton of the cubic NLSE.
Equation (13) also possesses the kink (the hole center has zero intensity) and bubble (having a nonzero value at the hole center) solitons 39,50,51 , which can be expressed as and where β 4 = 1 D µ 2 − g 1 α 1 2β 5 , β 5 = g 1 α 1 ± √ g 2 1 α 2 1 +4g 2 α 2 µ 2 2µ 2 and β 6 = 3g 1 α 1 β 2 5 +6g 2 α 2 β 5 3g 1 α 1 β 5 +4g 2 α 2 . Figure 6 is a sketch for the profiles of kink and bubble, where B 2 is adopted rather than B(X 1 ) itself. B 2 tends to a nonzero constant 1/β 5 when X 1 → ±∞ . At the center of the soliton, i.e. X 1 = 0 , B kink is zero but B bubble is nonzero. It is also worth remarkable that the kink B dark has a kink shape and the bubble B bubble has a bell-like one. β −1 5 and β −1 5 − (β 5 − β 6 ) −1 can be regarded as the amplitude of the kink and bubble, respectively. When β 6 tends to zero, the bubbles can no longer be excited. For the kink, to avoid the singularity, it is required that both β 5 and β 6 should be positive. However, for the bubble, these conditions become β 5 > 0 and β 6 < β 5 , implying that β 6 can be negative. Note that β 6 equals to β 5 when g 2 = 0 , suggesting that this kind of bubbles do not exist for the cubic NLSE. Therefore, the emergence of bubbles is just due to the effect of three-body interaction.

Stability analysis. Linear stability analysis on bright solitons.
We now numerically investigate the linear stability properties of solitons given by the NCG method. The perturbation carrier wave can be written as a Bogoliubov expansion where is the eigenvalue of the normal mode, and " * ′′ denotes complex conjugation, |υ|, |w| ≪ 1 are infinitesimal normal-mode perturbations. Substituting this perturbed solution into Eq. (2) and linearizing, we find that these normal modes satisfy the following linear eigenvalue problem General speaking, it is rather difficult to solve this linear eigenvalue problem analytically and exactly. However, we can solve it numerically and efficiently by the finite difference method or the Fourier collocation method 49 .
The numerical results indicate that all the on-site gap solitons of Eq. (4) are linearly stable. To indicate that the three-body interaction strength g 2 can indeed affect the stability of the off-site solitons, Fig. 7 shows the linear spectrum of off-site gap solitons under various of parameters. The gap soliton shown in Fig. 7a is linearly www.nature.com/scientificreports/ unstable ( g 1 = −1, g 2 = 0 ) while it is linearly stable in Fig. 7b ( g 1 = −1, g 2 = 0.6 ). Similarly, with different g 2 , the soliton in Fig. 7c ( g 1 = 1, g 2 = 0 ) is linearly unstable while it is linearly stable in Fig. 7d ( g 1 = 1, g 2 = −0.45 ) .
To deeply understand the effects of g 2 on the stability properties, Fig. 8 shows the maximum growth rate of perturbation m = max[Re( )] for the gap solitons versus g 2 under different chemical potential µ ( V 0 = 2, q = 0.1 ). The corresponding profiles of gap solitons are similar as Fig. 7b or d. From Fig. 8, one can see that the unstable gap solitons become stable if g 2 increases over a critical value when other parameters are fixed. Thus, one can change the stability property of gap solitons by adjusting the three-body interaction strength in experiments.   2), where the initial condition is taken as the random perturbed soliton obtained by the NCG method. The timesplitting Fourier Spectral method 52,53 is adopted to make the time evolution. This method has high accuracy and can guarantee that the number of particles is conserved. Figure 9c and f are the contour plots of |�(x, t)| . From which, it is clearly that the two dark solitons are dynamically unstable, which agrees well with the linear stability analysis. A natural question is whether there exists stable dark solitons (kink or bubble) in the semi-infinite gap.
To answer this question, we have made lots of numerical calculations for various values of g 1 < 0 and g 2 > 0 . Unfortunately, we failed to find the stable dark solitons for this case. In Ref. 54 , the authors found that the static bubble solitons are always unstable, which is also in well agreement with our conclusion. Are there any stable dark solitons in the BEC system? Excitingly and interestingly, we indeed found stable kinks under certain parameters. Figure 10a and d illustrate the profiles of dark solitons (red lines) given by NCG method for V 0 = 2, q = 0.8, µ = 1.2, g 1 = 1, g 2 = 0.5 , in which µ lies in the first gap. Figure 10b and e are the linear stability spectrum for the dark solitons shown in (a) and (d), respectively. One can see that the kink soliton is linearly stable, while the bubble has two unstable modes, thus it is linearly unstable. Figure 10c and f are the contour plots of |�(x, t)| . From which, it is clearly that the kink soliton is dynamically stable, while the bubble is dynamically unstable, which also agrees well with the linear stability analysis.
When the above nonlinear evolution is performed, the spacial interval is truncated from (−∞, +∞) to (-32T x , +32T x ), where T x is the periodicity of the external potential ( T x ≈ π when q = 0.1 ). Then the interval is divided into 8192 grids uniformly. We think that the length of the interval is large enough and the numerical method can give us the satisfactory results. On the other hand, we also think that the numerical error near the boundaries maybe large. So, we only pay attention on the interval (−20, 20) when plotting the results (see Figs. 9, 10c and f). We believe that the numerical results have high accuracy in this smaller interval. From Figs. 9 and 10, we have the reason to suppose that the modulus q of the external potential may be an important parameter for the stability property of the dark solitons. To confirm this conclusion, Fig. 11 shows m for the dark solitons versus the modulus q under different chemical potential µ . The corresponding profiles of dark solitons are similar as those shown in Fig. 9a or d. From Fig. 11a, we see that all the bubbles are linearly unstable. When q is near 1 (but less than 1), m decreases rapidly as q increases, implying that increasing the www.nature.com/scientificreports/ modulus of the external potential can weaken the instability of the bubble solitons. In Ref. 54 , the authors also found that the static bubble solitons are always unstable. Here, we have the same conclusion. However, for the kink solitons (see Fig. 11b), there exists a critical value of q c µ for a given µ . When q < q c µ , the kinks are linearly unstable, otherwise they are linearly stable. Take a instance, q c µ ≈ 0.7 for µ = 1.2 and q c µ ≈ 0.61 for µ = 1.4 , suggesting that q c µ depends upon µ . Thus, one can change the stability properties of kinks by adjusting the modulus of external potential in experiments.
It is worth remarkable that m oscillates as q increases in Fig. 11b, which makes the figure does not look as "pretty good" as Fig. 11a. This is because, for kink solitons, m are relatively smaller ( 10 −3 ) than those of bubble-like solitons ( ∼ 10 −1 ), implying that the perturbation increases rather slowly. The oscillation may be due to the numerical errors. In fact, it is difficult to identify that a very small m is caused by the soliton instability  www.nature.com/scientificreports/ or induced by the numerical errors. This trouble cannot be overcome even if one uses the long-time nonlinear evolution. In general, the soliton can be regarded as linearly stable within a relatively short period of time if m is small enough.

Conclusion
In summary, we have analytically and numerically investigated the stability properties of bright and dark solitons in a quasi-1D BEC with three-body interaction loaded in a Jacobian elliptic sine potential. Bright and dark solitons are numerically found by the NCG method. A stationary nonlinear Schrödinger equation is derived to describe the profiles of solitons via the multi-scale technique. Linear stability analysis indicates that the threebody interaction strength has distinct effect on the stability properties. Especially, such a nonlinear system supports the so-called dark solitons (kink or bubble), which can be excited not only in the gap, but also in the band. The bubbles cannot be excited if the three-body interaction is absent and they are always unstable. Both stable and unstable kinks, depending on the physical parameters, can be excited in the BEC system.