{'fraction_non_alphanumeric': 0.05189119065958568, 'fraction_numerical': 0.014999797300036486, 'mean_word_length': 4.371951219512195, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stability of gap soliton complexes in the nonlinear Schr\"odinger equation with periodic potential and repulsive nonlinearity

The work is devoted to numerical investigation of stability of stationary localized modes ("gap solitons") for the one-dimentional nonlinear Schr\"odinger equation (NLSE) with periodic potential and repulsive nonlinearity. Two classes of the modes are considered: a bound state of a pair of in-phase and out-of-phase fundamental gap solitons (FGSs) from the first bandgap separated by various number of empty potential wells. Using the standard framework of linear stability analysis, we computed the linear spectra for the gap solitons by means of the Fourier collocation method and the Evans function method. We found that the gap solitons of the first and second classes are exponentially unstable for odd and even numbers of separating periods of the potential, respectively. The real parts of unstable eigenvalues in corresponding spectra decay with the distance between FGSs exponentially. On the contrary, we observed that the modes of the first and second classes are either linearly stable or exhibit weak oscillatory instabilities if the number of empty potential wells separating FGSs is even and odd, respectively. In both cases, the oscillatory instabilities arise in some vicinity of upper bandgap edge. In order to check the linear stability results, we fulfilled numerical simulations for the time-dependent NLSE by means of a finite-difference scheme. As a result, all the considered exponentially unstable solutions have been deformed to long-lived pulsating formations whereas stable solutions conserved their shapes for a long time.


Introduction
In the past few decades the nonlinear Schrödinger equation (NLSE) with an additional non-autonomous linear term has become one of the actual and challenging physical problems which has been studied worldwide by different scientific groups. This equation can be written in the following three-dimensional form: From a physical viewpoint, equation (1) is related to the models of plasma physics [8,12], nonlinear optics [15,21,20,26,27] and Bose-Einstein condensation theory [13]. It is worth mentioning that in the latter context, equation (1), also called the Gross-Pitaevskii equation (GPE), has become especially relevant after experimental observation of the Bose-Einstein condensate in 1995 [6,11,14]. Such a state of matter has been predicted in 1924 by Einstein and Bose [10,16]. In the context of the meanfield theory of the Bose-Einstein condensate, the term |U | 2 describes the local density of the condensate. The function V (x, y, z) has the meaning of an external potential which allows one to confine the condensate spatially. Parameter σ describes interparticle interactions in the condensate: the value σ = +1 corresponds to repulsive interactions between atoms whereas σ = −1 corresponds to attractive interactions. The one-dimensional version of (1): describes an elongated (cigar-shaped) cloud of the condensate. An important set of solutions of (2) is a class of so-called stationary modes having the following form: where µ has the meaning of the chemical potential of the condensate. Stationary modes u(x) satisfy the localization condition: lim x→±∞ u(x) = 0.
In what follows, it is assumed that u(x) is a real-valued function [5]. Besides, the parameter σ is chosen according to the repulsive case of interactions, σ = +1, and the potential V (x) is fixed in the model form V (x) = −V 0 cos 2x. Here, the parameter V 0 is related to the depth of potential wells. Substituting (3) into (2), one arrives at the following time-independent equation: Due to the localization condition (4), for |x| ≫ 1 equation (5) can be replaced by the linear Schrödinger equation: It is well-known [7] that the spectrum of problem (6) has a band-gap structure which is depicted in Figure 1.
There is a countable number of bands separated by bandgaps. Such a band-gap structure of linearized equation imposes several restrictions for existence of localized solutions of (5). If the point (µ, V 0 ) lies in a band, equation (6) has no solutions that tend to zero at +∞ or at −∞. It means that the localization condition (4) for equation (5) for these values µ and V 0 cannot be satisfied, except for the trivial solution u(x) ≡ 0. Therefore, the stationary localized modes (3) can be obtained only in bandgaps. For this reason, the localized solutions of (5) are also called gap solitons.
The simplest localized solution of (5) is called fundamental gap soliton (FGS). This solution represents a single density hump with decaying tails. It was proven in [2] that more complex solitons can be considered as compositions of several FGSs situated at different potential wells. More precisely, all solutions of (5) situated in Region 1 in Figure 1 can be coded by means of bi-infinite symbolic sequences of the finite alphabet A = {+, 0, −}. The symbol "+" (or "−") at the n-th entry of the bi-infinite sequence (n runs over the all integers from −∞ to +∞) indicates that the n-th potential well is occupied by the FGS (or by the FGS taken with negative sign). Respectively, the zero symbol means that the corresponding potential well is "empty". Thus, the code of the FGS has the form (. . . , 0, 0, +, 0, 0, . . .), and the code of the FGS taken with negative sign has the form (. . . , 0, 0, −, 0, 0, . . .). More generally, the code of an arbitrary gap soliton "starts" and "ends" with the infinite number of zero symbols, i.e. has the form (. . . , 0, 0, s 1 , s 2 , . . . , s n , 0, 0, . . .), 0 ∈ A, s i ∈ A, s 1 = 0, s n = 0. These infinite sequences of zero symbols describe the asymptotically decaying tails of the soliton. For the sake of brevity, in what follows we omit the zero symbols situated at the soliton tails. Then the code of the FGS can be written as (+). The code (++) corresponds to a soliton composed of two neighbor FGSs, both taken with positive sign, and the code (+−) describes a composition of two FGSs taken with different signs. The code (+0+) describes a bound state of two FGSs separated by an "empty" potential well.
An important property of a gap soliton is its stability which indicates the robustness of the soliton against perturbations. Various articles have reported analytical and numerical results about stability of gap solitons of (5). Let us list the known results about stability of gap solitons of (5) situated in Region 1 in Figure 1. In was claimed in [24] that the small-amplitude FGS having the code (+) is linearly stable. However, in [24] it was also suggested that the FGS of greater amplitudes may suffer weak oscillatory instabilities caused by complex eigenvalues in the linear spectrum [see the spectral problem (9)]. These instabilities were found in the same paper for the case of attractive nonlinearity [i.e., the case σ = −1 in (2)]. Turning to the large amplitudes, the stability of the FGS was examined numerically in [22] where the authors concluded that it is a stable mode as well. Besides, the authors of [22] mentioned their finding of complex eigenvalues with a small non-zero real parts for (9), but they attributed them to the insufficient accuracy of their computations. Recently, in [19] the stability of the FGS was studied using the Evans function approach, and it was concluded that the spectrum of the FGS does include complex eigenvalues with small real parts which leads to weak instabilities. Therefore, the FGS typically is either stable or exhibits weak oscillatory instabilities.
The next basic solution of (5) has the code (+−). It was studied analitically in [18,24] using the smallamplitude approximation and turned out to be unstable (the spectrum includes a pair of real eigenvalues which implies so-called exponential instability). For large amplitudes, this solution was studied numerically [22] and turned out to be exponentially unstable as well.
The stability of more complex modes of (5), double-humped and triple-humped gap solitons with the codes (++) and (+ + +) (also called Truncated-Bloch-wave solitons), was checked numerically in [28] where these solutions were reported to be stable. However, it was shown recently in [19] that these solutions also may suffer from small oscillatory instabilities.
Two families of separated FGSs having the codes (+0+) and (+0−) have been studied numerically in [22] where they were found to be oscillatory unstable. However, the stability of the gap solitons with codes (+0 . . . 0+) and (+0 . . . 0−), i.e. having the symbols "+" or "−" separated by some number of symbols "0", up to the moment have not been studied in detail. An interesting analysis of bifurcations of these solutions can be found in [1]. In the present study we focus our attention on compexes of gap solitons of (5) having the codes (+0 . . . 0+) and (+0 . . . 0−) with different number of separating periods. We do not consider the codes (−0 . . . 0−) and (−0 . . . 0+) due to the symmetry of (5). We study numerically the linear stability of these modes and check it by direct integration of equation (2) using an appropriate numerical scheme. In Section 1 we describe the linear stability problem in detail and describe the required numerical methods. In Section 2 we present the main outcomes of this paper and in Section 3 we summarize the results of the work. x

Linear stability problem
From a physical viewpoint, the stability is an important property of solutions of equation (2). We study the linear (spectral) stability following traditional approach [29]. Let u(x) be an arbitrary solution of (5) satisfying the boundary conditions (4). This solution can be regarded as the initial condition U (x, 0) for equation (2). Adding a small perturbation to the function u(x) and substituting it to (3), one has the following perturbed solution of (2): where a(x) and b(x) are real and small enough, |a|, |b| ≪ 1. Substituting (7) into (2) and omitting all high-order terms with respect to a(x) and b(x), we come to the following system: Introducing the linear operators (8) can be rewritten in the form of the spectral problem: Generally, the spectrum of (9) consists of an essential and a discrete parts. Since the linear stability of considered modes is related to the latter part, let us mention some properties of discrete eigenvalues of (9). The eigenvalue λ = 0 is double, with eigenvector (0, u) T and generalized eigenvector (∂u/∂µ, 0) T that obey: Re λ respectively. This fact allows one to check the accuracy of numerical techniques in a simple way. We use it to check the Evans function method which is described in [9] in detail. Besides, if the spectrum of (9) contains an isolated eigenvalue λ, it also contains the values −λ, λ * and −λ * where the asterisk means the complex conjugation. Solving the problem (9) one concludes whether the regarding mode (3) is stable or unstable. In the case the whole spectrum is pure imaginary, the solution (3) is said to be linearly stable. On the other hand, if there is at least one eigenvalue with non-zero real part, then the perturbation grows exponentially and the solution (3) is regarded to be unstable.

Numerical methods
In order to construct the gap solitons of (5), we applied the modified shooting algorithm [4]. This algorithm allows one to compute the profile of a required mode by its symbolic representation. It includes solving of a Cauchy problem for equation (5) which was implemented using the Runge-Kutta 4-order scheme.
Concerning the numerical solution of the spectral problem (9), we used two alternative algorithms. One of them is the well-elaborated Fourier collocation method (FCM, see [29] for details). It is well-suited for fast and accurate detecting of relatively strong instabilities. Hovewer, its accuracy may not be sufficient for tracing relatively weak instabilities associated with eigenvalues with small nonzero real parts. For accurate detection of such weak oscillatory instabilities, we take advantage of the second algorithm based on the Evans function (EFM) that allows one to compute complex eigenvalues with tiny real parts with high precision. On the other hand, the realization of EFM is quite sophisticated (see [9]).
The linear stability results for gap solitons of (5) obtained by both FCM and EFM allow us to predict the dynamics of corresponding solution of (2) in time. To check such a prediction, we integrated (2) by means of a finite difference scheme described in [25]. In the case the solution U (x, t) undergoes any type of instabilities, the profile of initial condition U (x, 0) is expected to undergo deformation due to the growth of the perturbation whereas the profile of a stable solution is expected to approximately preserve its shape for any period of time.

Results
Using the numerical schemes described above we fulfilled two series of numerical experiments for equations (2) and (9). The Sections 2.1 and 2.2 below are devoted to two classes of gap solitons consisting of two separated in-phase and out-of-phase FGSs, respectively. Notice that since the profile of each gap soliton depends on the parameters µ and V 0 , for a fixed V 0 there is one-parametric family of gap solitons having the same code. Since only two gap solitons with the codes (+) and (+−) can be continued to the lower boundary of the first bandgap, µ = µ − , the solutions with all other codes "die" in bifurcations for some greater values µ > µ − .

The stability of in-phase FGSs
We started our investigation with bound states consisting of two in-phase FGSs separated by an arbitrary number of empty potential wells. Such solutions have the codes of the form (+0 . . . 0+) and (−0 . . . 0−) with various numbers of zero symbols "0" in the middle. The simplest representatives of this class are depicted in Figure 2(a)-(c) for parameter µ = 1. On the one hand, they all are situated in Region 1 in Figure 1, far away from the upper bandgap edge, which implies the absence of weak oscillatory instabilities in their spectra [19]. On the other hand, they feature different stability properties. Looking at corresponding spectra presented in Figure 3 N )]. Therefore we choose another characteristic, which is free of this disadvantage and allows us to display the families properly: Notice that if a solution u(x) is either even or odd then its characteristic M is equal to zero. Therefore we shift the solutions shown in Figure 2 in such a way that one of two FGSs is situated in zero potential well, x = 0. The corresponding families on the plane (µ, M ) are depicted in Figure 5(a) by black curves. Thin and thick segments of a curve are related to regions of linear stability and instability of gap solitons, respectively. These curves are situated in the white spacious region of Figure 5 which represents the first bandgap. We found that the whole families (a) and (c) are linearly unstable due to the presence of one pair of real eigenvalues in spectrum of (9). In the same time, the family (b) is stable for µ − < µ < µ * where the threshold value is µ * ≈ 1.57. If µ is greater than µ * , the family (b) undergoes weak oscillatory instabilities which are similar to those described in [19]. For instanse, the spectrum of solution (b) includes the quartet of complex eigenvalues with real part Re λ ∼ 10 −5 in the case µ = 1.7 and Re λ ∼ 10 −4 as µ = 2.3. In addition, the inset presented in the left top corner in Figure 5(a) zooms in the bifurcation edges of considered families.
To summarize, our numerical observations allow one to conjecture that if the number of empty potential wells between two in-phase FGSs (zero symbols in middle of the code) is odd, the composed gap soliton is unstable due to a strong exponential instability. Otherwise, it is either stable or exhibits weak oscillatory instability.

The stability of out-of-phase FGSs
Now let us consider the second class of bound states consistsing of two separated out-of-phase FGSs with the codes of the form (+0 . . . 0−) and (−0 . . . 0+). Let us fix again the depth of the potential, V 0 = 6. One can find the simplest states from this class in Figure 2(d)-(f) for parameter µ = 1. Using FCM and EFM, we found that the gap soliton (d) and (f) are linearly stable, see Figure 3(d), (f). On the contrary, the gap soliton (e) was found to be unstable, see Figure 3(e).
The families of bound states (d)-(f), see Figure 5(b), can be extended numerically from the value µ = 1 to the smaller and greater values. But none of these families can be continued to the lower bandgap edge. The bifurcation points of solutions (d)-(f) are µ ≈ −2.774, µ ≈ −2.765 and µ ≈ −2.777, respectively. Our numerical finding here is that the families (d) and (f) are linearly stable for all values of µ up to some threshold µ * ≈ 1.57. Notice that it is similar to the case of two in-phase FGSs. Moreover, in [19] the equality of the stability thresholds among another families of gap solitons with codes (+), (++) and (+ + +) was also observed. The existance of such threshold may be related to the FGS itself which is an elementrary entity for more complex gap solitons. For greater values of µ, the families of two out-of-phase FGSs exhibit weak oscillatory instabilities similar to [19]. The orders of the real parts of the unstable eigenvalues are approximately the same as in the case of two separated in-phase FGSs described above. The family (e) is unstable for the entire interval of its existence. Looking to these linear stability results, one can conjecture that the solutions from the second class are linearly stable if two FGSs are separated by an odd number of empty potential wells and the corresponding value µ is situated far enough from the upper bandgap edge. If the number of empty sites is even, then these solutions are supposed to be unstable.
To study the dependence between instability of a gap soliton and the number of empty potential wells in its middle, we computed the spectra of five solutions with codes (+0+), (+00−), (+000+), (+0000−) and (+00000+). All these gap solitons exhibit a strong exponential instability caused by a pair of real eigenvalues. In Figure 4 one can find the results of this experiment. Real parts of unstable eigenvalues decay exponentially as the number of empty potential wells grows. It may be related to the coupling and uncoupling of two separated FGSs due to short and long distance of separating, respectively. If the separated FGSs are uncoupled, then the corresponding gap soliton is linearly stable.

The dynamics of time-dependent GPE
In order to check the predictions of linear stability analysis, we integrated numerically the time-dependent NLSE (2) with initial conditions depicted in Figure 2(a) and (d) with different parameters µ and V 0 . As we mentioned above, we used a finite difference scheme with spatial step h ≈ 0.1 and time step τ = 0.01 in the domain [−10π; 10π] × [0; 40000]. In order to eliminate the reflection from the boundaries, we added artificial absorbing layers for x > 20 and x < −20. Besides, we perturbed the initial condition by means of 3-percent perturbation, i.e. U (x, 0) → 1.03 × U (x, 0). The results of this experiment are shown in Figure 6. Panels (a) and (b) correspond to mode (+0−) for parameter µ = 1 and µ = 2, respectively. One can see that the results of numerical integration of (2) are in a good agreement with linear stability analysis considered in the previous subsections. The gap soliton computed in the panel (a) preserved its shape for the whole time of computation while that in the panel (b) has been deformed due to weak oscillatory instabilities. Panel (c) exhibits the evolution of the mode (+0+) for µ = 1. It has been confirmed by the simulations that this mode is unstable. Notice that the characteristic time necessary for the growing perturbation to manisfest itself is much smaller for the exponential instability than for the weak oscillatory instability.
One more observation is as follows. The temporal evolution of the unstable mode (+0+) results not in its destroying, but in forming of some pulsating formation resembling a breather. Simulations for the codes (+00−), (+000+) and different parameters µ and V 0 yield the similar result. We conjecture that this phenomenon is typical and every exponentially unstable bound state, e.g. one with the code (+0 . . .

Conclusions
In the paper, we studied the stability of gap solitons for the one-dimentional nonlinear Shrödinger equation (NLSE) with cosine potential (2). It is known, that under certain conditions [2] each gap soliton of this equation can be coded by means of a bi-infinite sequence of symbols of some finite alphabet. In out study, the parameters µ and V 0 (µ is the chemical potential, see the anzats (3), V 0 is the depth of the potential) were chosen from the Region 1 in Figure 1 where three-symbols alphabet, "−", "0" and "+", for the coding is relevant. The gap solitons of interest are the bound states of a pair of fundamental gap solitons (FGSs) separated by some number of periods of the potential. In-phase bound states correspond to the codes (+0 . . . 0+) and (−0 . . . 0−) (the top row in Figure 2) and out-of-phase ones have the codes (+0 . . . 0−) and (−0 . . . 0+) (the bottom row of Figure 2).
In order to study linear stability of above-mentioned classes of gap solitons, we solved numerically the spectral problem (9) using the Fourier collocation method (FCM) and the Evans function method (EFM). It follows from our computations that some of the considered gap solitons are stable, some of them exhibit strong (exponential) instabilities and other ones undergo weak oscillatory instabilities. The modes with the codes (+0−), (+00+) and (+000−) are linearly stable for all values of µ below a certain threshold. As µ exceeds this threshold, these modes exhibit weak oscillatory instabilities. On the contrary, the gap solitons with the codes (+0+), (+00−) and (+000+) turned out to be unstable due to the presence of a pair of real eigenvalues in their spectra for all considered parameters µ. As the result, we conjecture that two separated in-phase (out-of-phase) FGSs are linearly stable (unstable) for an even (odd) number of empty potential wells between them. Also, two in-phase (out-of-phase) FGSs are linearly unstable (stable) for an odd (even) number of empty sites between them. The direct integration of (2) by means of a finite-difference scheme supported the linear stability results. While all the results presented above correspond to V 0 = 6, we have also checked that the main conclusions of our study remain valid for several other values of V 0 (such as V 0 = 5 and V 0 = 7).
It is interesting that the results of this study confirm the known results on stability/instability of intrinsic localized modes (ILM) for the discrete nonlinear Schrödinger equation (DNLSE). It is known that the DNLSE can be regarded as an approximation for the NLSE with periodic potential and ILM can be used for qualitative description of gap solitons in the first bandgap [3,17] (so-called tight-binding limit). Rigorous statement about the number of unstable eigenvalues for ILM can be found in [23]. Theorem 3.6 of [23] connects the number of unstable eigenvalues for an ILM with the number of flips in its code. In [23] this result was formulated for the case of the DNLSE with focusing nonlinearity and is valid near so-called anticontinuous limit only. Translating this statement to case of the defocusing DNLSE can be made by means of the staggering transformation. However, using the analogy between the DNLSE and the periodic NLSE, one can expect instability of the bound states for the codes (+0 . . . 0+) and (−0 . . . 0−) with odd number of zero symbols between nonzero symbols and for (+0 . . . 0−) and (+0 . . . 0−) with even number of zero symbols, exactly as it is claimed in this paper. At the same time, the analogy with the discrete model fails to reveal quite delicate oscillatory instabilities described above that were not found in the case of ILM.
A further generalization of this study may be associated with the following question: does there exist a general relation between the stability of a gap soliton and its code? We believe that the numerical experience summarized in [29] and in the present paper may be extended to gap solitons with more complex codes, as well as to the gap solitons from higher bandgaps. However, this interesting issue lies beyond the scope of this paper.

Acknowledgment
Author is grateful to Prof. Georgy L. Alfimov and Dr. Dmitry A. Zezyulin for a careful reading of the manuscript, significant remarks and discussions.