{'fraction_non_alphanumeric': 0.03897785779715089, 'fraction_numerical': 0.016582962190846207, 'mean_word_length': 4.507044703426683, 'pattern_counts': {'":': 0, '<': 2, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Oscillatory instabilities of gap solitons in a repulsive Bose-Einstein condensate

The paper is devoted to numerical study of stability of nonlinear localized modes ("gap solitons") for the spatially one-dimensional Gross-Pitaevskii equation (1D GPE) with periodic potential and repulsive interparticle interactions. We use the Evans function approach combined with the exterior algebra formulation in order to detect and describe weak oscillatory instabilities. We show that the simplest ("fundamental") gap solitons in the first and in the second spectral gaps can undergo oscillatory instabilities for certain values of the frequency parameter (i.e., the chemical potential). The number of unstable eigenvalues and the associated instability rates are described. Several stable and unstable more complex (non-fundamental) gap solitons are also discussed. The results obtained from the Evans function approach are independently confirmed using the direct numerical integration of the GPE.


I. INTRODUCTION
The Gross-Pitaevskii equation (GPE), describes the meanfield dynamics of a quasi-onedimensional Bose-Einstein condensate (BEC) confined in the potential V (x) [1]. In Eq. (1), Ψ = Ψ(t, x) is the complex-valued macroscopic wavefunction of the condensate. The squared amplitude of the wavefunction |Ψ(t, x)| 2 describes the local density of the BEC, while the gradient (arg Ψ(t, x)) x describes the velocity of atoms of condensate. The nonlinear term σ|Ψ| 2 Ψ takes into account the interactions between the particles. The case σ = 1 corresponds to repulsive interparticle interactions, while σ = −1 describes attractive interactions. Both these cases are of physical relevance, but in what follows, we mainly focus on the repulsive case, σ = 1. It is assumed that the potential V (x) is periodic, which corresponds to the optical confinement of the BEC [2,3]. An important class of solutions of the GPE (1) corresponds to the stationary modes which can be represented in the form Ψ(t, x) = e −iµt u(x), where µ is a real parameter having the meaning of the chemical potential of the BEC. Function u(x) satisfies the conditions of the spatial localization lim x→±∞ u(x) = 0. ( Without loss of generality one can assume that u(x) is real-valued [4]. Then u(x) can be found from the stationary GPE If the potential V (x) is periodic, the nonlinear modes satisfying (2) values of µ corresponding to these solutions lie in the spectral gaps of the linear Schrödinger equation [9] u xx + (µ − V (x))u = 0.
The simplest class of the gap solitons are fundamental gap solitons (FGSs) [5,7,[10][11][12]. Under the repulsive nonlinearity, in the first gap there exists one family of FGSs which are characterized by the presence of a single dominating peak localized in one well of the potential V (x). A variety of more complex (or higher-order ) solitons includes truncated Bloch waves [13,14] (which consist of several in-phase peaks placed in a row), various asymmetric states, and complex bound states of two (or more) well-separated waves [15], etc. In spite of their rich diversity, under certain (not very restrictive) conditions all possible gap solitons in a repulsive BEC can be viewed as complexes of FGSs and classified using an alphabet consisting of a few symbols [16]. Specifically, if the lattice is deep enough, then all the gap solitons in the first gap can be put into a one-to-one correspondence with the set of bi-infinite sequences of symbols from a three-symbol alphabet. In simple terms, these symbols denote the presence or the absence of the FGS (taken with plus or minus sign) in a potential well situated on the period of the potential V (x). For instance, the truncated Bloch waves [13,14] consisting of several in-phase peaks placed in a row can be viewed as complexes of single-hump FGSs. For classification of the gap solitons in the second spectral gap, an alphabet of five symbols is necessary, and so on.
An important property of a gap soliton is its stability, since only dynamically stable modes are likely to be experimentally feasible. The literature about stability of gap solitons is rather abundant [6-8, 11-13, 17, 19, 20]. The most relevant for our study outcomes for the case of repulsive interactions (σ = 1) can be summarized as follows. Significant part of the studies concluded that the single-hump FGSs are stable in the first gap [6,7,[10][11][12] and in the second gap [10][11][12]. Regarding the higherorder states consisting of two or three in-phase peaks, they have been reported unstable near the upper band edge in [7] in the first gap. However, these states have been found to be stable both in the first [11,13] and in the second [11] gap if the lattice depth is large enough.
The results listed above have been obtained on the basis of numerical studies of stability. In the meanwhile, it is recognized that the numerical analysis of stability of the gap solitons is quite a delicate problem. A standard approach to the stability relies on the linear (or spectral) stability technique which reduces the stability question to a study of the spectrum of a certain linear operator. Depending on the character of unstable eigenvalues, the instability typically manifests itself either as a purely exponential instability (when the unstable eigenvalues have zero imaginary parts) or as an oscillatory instability (OI) (when the unstable eigenvalues are complex with nonzero imaginary parts). While the instabilities of the former type are relatively simple to detect [8,19], the accurate tracing of OIs is much more challenging [7,8]. As a result, the information about OIs of gap solitons in optical lattices is rather scarce. The absence of information on OIs for the simplest one-hump gap solitons in GPE is especially remarkable in view of well-known OIs of Bragg gap solitons in nonlinear Dirac equations [21][22][23]. The latter system can be deduced from the GPE with a shallow periodic potential using an asymptotic multiple-scale expansion [17,18], and therefore the results about OIs of the solitons for Dirac system seem to be not consistent with the stability results for the single-hump FGS mentioned above.
The numerical difficulties arising in the analysis of the OIs of the gap solitons are related to several issues. First, the rates of OIs are typically quite small [8,19]. Another difficulty results from poor localization of the gap soliton and (or) of the eigenfunction associated with an unstable eigenvalue. This situation typically takes place when the chemical potential µ is close to the gap edge. It requires unpractically wide computational windows or a particularly accurate treatment of the boundary conditions. Some of these difficulties can be overcome using the Evans function approach which was employed in [8] to trace OIs of gap solitons in the attractive condensates. It was further demonstrated in [24] that the numerically accurate evaluation of the Evans function requires a careful treatment of the stiffness issue which arises for some values of the complex argument of the Evans function. The stiffness problem can be fixed if one redefined the Evans function using the exterior algebra formalism [23,24]. This idea has been developed into a robust numerical technique which was demonstrated to provide reliable results even for relatively weak instabilities of gap solitons [23,24].
In the present paper, we use the Evans function approach combined with the exterior algebra formulation in order to reveal and describe weak OIs of FGS and higher-order gap solitons in the repulsive BEC. We focus on the first and second spectral gaps. In each gap, we consider the single-hump FGS and two higher-order solitons bearing two or three in-phase humps. The main outcomes of our numerical study can be outlined as follows.
1. In the first gap, all the considered solitons (including the single-hump FGS) are stable far from the upper band edge, but undergo OIs in the region near the upper band edge. The width of this instability region is quite significant: it occupies about 15%-20% of the width of the first gap.
2. In the second gap, all the considered solitons (including the FGSs) are, in general, unstable due to OIs. However, in a sufficiently deep potential, there exist intervals of µ where FGS are stable.
To the best of our knowledge, our results constitute the first explicit demonstration and detailed description of OIs for FGSs in the case of repulsive interactions σ = 1 in the first and in the second gap. On the other hand, our results advance the current understanding of the higherorder modes [13], since we show that they undergo OIs even in a deep potential.
In order to confirm the linear stability results, we have also performed a series of direct simulations of temporal behaviour of the solitons in the GPE (1). The results of these studies agree with the conclusions obtained from the linear stability analysis and display the slow decay of unstable gap solitons and the persistent evolution of the stable ones.
The rest of the paper is organized as follows. In Sec. II we briefly describe the families of gap solitons whose stability is the main subject of the present study. In Sec. III we present our main results on linear and nonlinear stability of the gap solitons. Section IV concludes the paper.

II. FAMILIES OF GAP SOLITONS
In our study, we use a prototypical example of the periodic potential in the form where real V 0 > 0 characterizes the depth of the lattice. The spectrum of the linear eigenvalue problem (4) with the potential (5) consists of one semi-infinite gap and a countable set of finite gaps [9] (see Fig. 1). In our study, we consider the gap solitons in the first and in the second gaps (no gap solitons exist in the semi-infinite gap in repulsive BECs). In each of the gaps, we consider three families of nonlinear modes: the fundamental gap soliton with the single dominating peak at x = 0 (the single-hump FGS), and two higher-order solitons, with two and three dominating peaks. In the terminology of [13], these higher-order solitons correspond to the truncated Bloch waves, and in terms of the coding approach [16] they are the bound states of two or three in-phase single-hump FGSs with codes (++) and (+ + +). Equation (3) also supports families of out-of-phase multi-hump solitons (with the codes (+−), (+ − +), etc. However, these solutions have been reported to suffer relatively strong exponential instabilities [19] and hence are not considered in our study.
Representative spatial profiles of the considered solutions for two different values of the chemical potential µ are shown in Fig. 2. Changing value of µ inside the gap, it is possible to construct numerically continuous families of gap solitons [7,8,[11][12][13]19]. These families can be visualized in the plane N vs. µ, where the squared L 2 -norm N = ∞ −∞ |u| 2 dx has the physical meaning of the number of particles in the condensate. The dependencies N vs. µ for V 0 = 6 are presented in Fig. 3(a). In the first gap, one of the families (corresponding to the single-hump FGS) bifurcates from the edge of the spectral band. The two other curves do not bifurcate from the band edge, but appear as a result of bifurcations which take place at a small but finite distance from the gap edge [15,19]. All the three curves can be continued to the upper edge of the first gap. Inside the band, the solitons do not exist. However, fundamental and multi-hump solitons can be found again when µ lies in the second gap. The considered families do not bifurcate from the edge of the second gap, but exist only as µ exceeds a certain finite threshold [11] (the values of these thresholds for the three considered families turn out to be very close, see the light gray shading in Fig. 3 which shows the interval in the second gap where the families do not exist.) The shapes of the solitons in the second gap are similar to their counterparts in the first gap, see Fig. 2.
In the second gap, one can also find two families of gap solitons bifurcating from the gap edge (not shown in Fig. 3) [8,20]. One of them is exponentially unstable [8], and another one (called subfundamental soliton in [12]) is unstable as the number of particles exceeds a certain threshold [11]. Since instabilities of these solutions are well-known, in what follows we do not incorporate them in our study and focus on the simplest FGS and solitons consisting of two and three in-phase FGSs whose instabilities have not been seen before.

III. STABILITY OF GAP SOLITONS
A. Statement of the problem and the numerical method Following to the standard procedure of the linear stability analysis, we consider a perturbed solution where u(x) describes a profile of a gap soliton, while a and b, |a|, |b| ≪ 1, describe real and imaginary parts of a small-amplitude perturbation. After substitution of (6) into the GPE (1) and neglecting higher-order in a and b terms, one arrives at the following linear eigenvalue problem where linear operators L ± are defined as and P = µ + V 0 cos 2x. If the spectrum of L is purely imaginary, then the amplitudes of perturbations a(x) and b(x) do not grow, and the soliton is said to be linearly stable. On the other hand, if at least one eigenvalue λ has nonzero real part, then the soliton is unstable. The largest real part of the eigenvalues characterizes the instability growth rate. The eigenvalue problem (7) has a double eigenvalue λ = 0. The corresponding eigenvector v 1 and the gener- Then Lv 1 = 0 and L 2 v 2 = 0. In (8), ∂ µ u = ∂u(x; µ)/∂µ is the partial derivative of the solution with respect to the chemical potential µ (provided that this derivative exists). It is also easy to see that if λ belongs to the spectrum of the operator L, then −λ and ±λ * also belong to the spectrum. Thus there are two typical scenarios of instability: (i) the instability may take place due to the presence of a pair of real eigenvalues λ and −λ and (ii) it may be caused by a quartet of unstable eigenvalues (λ, λ * , −λ, −λ * ), where real and imaginary parts of λ are nonzero. The instabilities of the type (ii) are known as oscillatory instabilities [8].
In spite of recent advances in the rigorous theory for analysis of the eigenvalue problem (7) [17,25], in the general case the description of its spectrum can be pursued only numerically. Numerical solution of the linear stability problem (7) is recognized to be a sufficiently challenging problem [7,8]. One of well-elaborated numerical approaches to this task is based on the the Evans function [8,[23][24][25] which is especially useful if the eigenfunctions associated with unstable eigenvalues of (7) are poorly localized. The Evans function f (λ) is an analytic function defined in the complex domain except the points of the essential spectrum of the linear stability operator L. The set of zeros of the Evans function coincides with the set of the isolated eigenvalues of the linear stability operator. Moreover, the order of the zero of the Evans function is equal to the algebraic multiplicity of the corresponding isolated eigenvalue.
According to the standard definition [25], the Evans function is a λ-dependent determinant whose columns are the vectors of stable and unstable manifolds for the linear stability problem (7). In this case the problem (7) is treated as a system of ODEs depending on the complex parameter λ, and these vectors should be computed solving the Cauchy problem for (7). However, it was demonstrated [24] that for certain values of the complex argument λ the system (7) is stiff and the direct evaluation of the Evans function cannot be fulfilled with the sufficient accuracy. The stiffness issue can be overcome by redefinition of the Evans function using the exterior algebra formalism [23,24]. Redefined in this way, the Evans function was used to study the stability of the surface gap solitons [24]. In the present study, we use the exterior algebra formulation of the Evans function to describe oscillatory instabilities of the fundamental and higher order gap solitons in the repulsive BEC. While the detailed explanation of the numerical approach can be found in [24], for the sake of self-containment of our work we describe briefly the main ingredients of the method in the Appendix A. Appendix B addresses technical details of implementation of the method. Additionally, in Appendix B we briefly compare the stability results obtained with the Evans function approach and with the Fourier collocation method which is another well-elaborated tool for computing the instabilities of gap solitons [19]. As follows from our results, the accuracy of the Fourier collocation method may be not sufficient to compute weak OIs. This confirms that the intricate Evans function technique is essential for the accurate tracing of OIs of the FGSs.

B. Linear stability results
Using the Evans function approach, we have examined in details stability of fundamental (single-hump) and multi-hump solitons in the first and the second spectral gaps of a lattice (5) with the fixed depth V 0 = 6. The results of this study are presented in Fig. 3. Then, in order to check that the obtained results are generic and to understand the effect of the depth of the lattice on the found instabilities, we have repeated the computations in a more shallow (V 0 = 3) and a deeper (V 0 = 9) lattices, see Fig. 4 and Fig. 5. The chosen depths of the potential are shown with vertical dashed lined in Fig. 1. Let us now proceed to the detailed presentation of the main outcomes of our work. Figure 3 shows the results of the first series of numerical experiments. Let us start from the FGSs in the first gap. As shown in Fig. 3(e), we have not found any unstable eigenvalue of FGS for µ < 1.586. Therefore, we conjecture that the FGSs are stable for the corresponding values of the chemical potential. However, for larger µ, we have identified three quartets of unstable eigen- values. Imaginary parts of the unstable eigenvalues are typically nonzero, which confirms to the oscillatory character of the instabilities [see Fig. 3(k)]. The first quartet of unstable eigenvalues exists only in a relatively narrow window 1.586 µ 1.625 [ Fig. 3(b)]. The other two coexist in a sufficiently wide instability window which starts at µ ≈ 1.633 at continues up to the upper gap edge. In the interval 1.625 µ 1.633 no instability has been found [see Fig. 3(b)]. Looking at the instability rates associated with the found OIs [ Fig. 3(e)], we notice that the real parts of the found eigenvalues Reλ do not exceed 10 −2 and typically lie between 10 −5 and 10 −3 .
Turning to multi-hump solitons in the first gap [ Fig. 3(f,g)], we observe that the stability picture is re- Proceeding to the solitons in the second gap, we observe that all three considered solutions are unstable in the entire domain of their existence. The one-hump FGSs [ Fig. 3(h)] are unstable due to the presence of two (or more, depending on the particular value of µ) quartets of unstable eigenvalues. Imaginary parts of these eigenvalues are presented in Fig. 3(l). Again, the instability rate is below 10 −2 . However, the instabilities become stronger for multi-hump solutions in the second gap [ Fig. 3(i, k)]: for these solutions |Reλ| reaches 10 −1 . In general, multihump solitons are unstable because of the presence of multiple coexisting quartets of unstable eigenvalues.
In order to examine the role of the lattice depth, we performed the second series of numerical runs extending the study onto solitons in a more shallow (V 0 = 3) and a deeper (V 0 = 9) potentials. Figure 4 presents the stability results for the shallow potential. We observe that the solitons in the first gap also feature two windows of instability, similar to their counterparts with V 0 = 6. For the fundamental family the first window of instability is situated at 1.394 µ 1.542. In this region, the spectrum includes two quartets of unstable eigenvalues. The second window of instability starts at µ ≈ 2.052 and continues to the upper gap edge. In this case, the instability is again associated with two quartets of eigenvalues. For two-hump and three-hump solitons the picture is qualitatively the same. The considered families of gap solitons do not exist in the second gap of the shallow lattice, and therefore the corresponding panels with the stability results are absent Fig. 4.
Turning to the deep lattice, V 0 = 9 [see Fig. 5], we again observe two intervals of instability of solitons in the first gap. The instability rate remains small and typically does not exceed 10 −3 . More interestingly, we can conjecture that solitons in the second gap are stable if µ belongs to certain narrow stability windows [ Fig. 5(h-j)]. For example, for the FGSs in the second gap, no unstable eigenvalues has been found near the lower gap edge µ 3.072 and in the narrow window 5.793 µ 5.949. Thus the increase of V 0 enhances the stability of FGSs in the second gap.

C. Solution of the time-dependent GPE
In order to check the linear stability results obtained with the Evans function method, we have also tested the nonlinear stability of the gap solitons by means of the direct integration of the GPE (1) using a semi-implicit finite difference scheme [26]. The equation was spatially discretized in a sufficiently large computational window x ∈ [−20π, 20π] with the grid step equal to π/256. We have additionally introduced absorbing (dissipative) boundary layers by adding nonzero imaginary part to the potential U (x) for x > 50 and x < −50. These absorbing layers mimic a setup in which the perturbations are removed from the system as they escape from the central region and reach the boundaries. In order to boost the development of eventual instabilities, the initial conditions have been taken in the form of slightly perturbed gap solitons. Specifically, we have used a 3%multiplicative perturbation: ψ → 1.03ψ. The temporal step was τ = 10 −3 .
In general, the results of the evolutional simulations are in the complete agreement with the predictions of the linear stability analysis. For stable solutions, the amplitude of the introduced perturbation does not grow, and the solution persists undistorted for indefinitely long time. For unstable solutions, the introduced perturbation grows slowly, which eventually leads to the breakup of the solution. In order to exemplify this pattern, in Fig. 6(a,b) we compare the dynamics of a stable and an unstable FGSs from the first gap. The drastic difference in their evolutions is especially well-visible in the plot with |Ψ(t, x = 0)| 2 , see Fig. 6(c). To illustrate the dynamics in the second gap, in Fig. 7 we compare the behavior of a stable FGS which exists in the second gap if the potential is deep enough [specifically, we chose V 0 = 9 as in Fig. 4(b)] with an unstable FGS with the same V 0 but different µ. Again, the striking difference in their behaviors is much better pronounced on the plot which shows the density at x = 0 [ Fig. 7(c)]. Looking at the long-term dynamics of solutions in Figs. 6 and 7, we observe that amplitude of unstable solutions decreases, and after a long transient period the shape of solutions features small-amplitude oscillations around some localized states. These oscillations persist for indefinitely long time.
The results for dynamics of multi-hump solitons also agree with the predictions of the linear stability analysis.

IV. CONCLUSION
To conclude, we have performed a systematic study of instabilities of gap solitons in a repulsive Bose-Einstein condensate trapped in a periodic potential. While this topic has been addressed in several previous studies, the complete picture of stability of gap solitons have not been completely understood yet. One of the reasons for this is related to the fact that the accurate study of instabilities of gap solitons can be a challenging problem which requires a sufficiently advanced numerical tool. In our study we used the Evans function approach combined with the external algebra formalism. The main advance made by our work is related to the explicit presentation of oscillatory instabilities of the fundamental soli-tons. While the possibility of existence of these instabilities has been hypothetically suggested, we present the first, to the best our knowledge, explicit demonstration of unstable eigenvalues, describe the regions of stability and instability and discuss the typical instability rates. Next, we have demonstrated that the found oscillatory instabilities are rather generic and can be also found in a more shallow and a deeper optical lattices. On the other hand, it was in particular found that the increase of the lattice depth can enhance the stability of FGSs in the second gap: the entire branch of solutions is unstable for the depth equal to V 0 = 6, but displays windows of stability for V 0 = 9.
Besides the fundamental single-hump solitons, we have also examined stability of two-hump and three-hump solutions and demonstrated that such solutions from the first gap suffer oscillatory instabilities even in a deep potential. Finally, we have performed a series of numerical runs with the direct integration of the Gross-Pitaevskii equation and observed how the predicted oscillatory instabilities manifest themselves during the temporal evolution.
In our opinion, the results of this work may be regarded as a first step in the study of the hypothetical relation between the code of a gap soliton (in sense of the paper [16]) and the spectrum of the corresponding stability problem. For the solitons from the first bandgap, the possibility of this relation has been discussed in [19], focusing mainly on the exponential instabilities of solitons coded by heuristically constructed sequences of "+" and "−" symbols. Another potential subject for further studies is the analysis of bifurcations of gap solitons, in terms of their codes, which can be also strongly related to their stability properties. Here we would mention the recent result of [27] about the connection between the modes of Discrete Nonlinear Schrodinger equation and gap solitons of the GPE. This connection allows to identify the codes of solitons that "annihilate" each other, see [28] for the table of such bifurcations.

Acknowledgments
The work of DAZ was supported by the FCT (Portugal) through the grant No. UID/FIS/00618/2013.

Appendix A: Definition of the Evans function
Following to [24], we rewrite the linear stability eigenvalue problem (7) as follows where p = (a+ib)/2 and q = (a−ib)/2. Further, Eq. (A1) can be rewritten as an ODE system, where λ plays the role of a parameter: Due to the localization lim x→±∞ u(x) = 0, in the limit x → ±∞ Eq. (A2) decouples into a pair of Mathieu equations T . The stable and unstable manifolds can be found numerically using the monodromy matrix associated with Eq. (A3). In a similar way, one can compute the monodromy matrix associated with Eq. (A4) and identify its stable T manifolds. Using vectors {v ± } and {w ± }, we can find an unstable manifold E u of Eq. (A2) and its stable manifold E s : Using (A5)-(A6) the Evans function is defined as the determinant of a 4 × 4 matrix [24] f The definition (A7) is, in principle, suitable for numerical evaluation of the Evans function for different values of the complex argument λ. To this end, one can start from a sufficiently large x = ±L, L ≫ 0, solve the ODE system from −L to 0 and from L to 0 in order to obtain the necessary values at x = 0, and compute the necessary determinant. However, this procedure suffers from a stiffness problem which manifests itself when the growth/decay rates of the two vectors in the manifold of interest (stable or unstable) are largely different. This issue can be overcome by a reformulation of (A2) in the exterior algebra. Following to [24], we introduce a 6 × 6 matrix function and instead of solving system (A2), we introduce a new ODE system We compute two solutions of (A9), U u and U s , which are correspond to initial conditions Finally, the Evans function is redefined as and ·, · stands for the standard inner product in C 6 with the complex conjugation of the first argument.

Appendix B: Details of numerical implementation and accuracy issues
The evaluation of the Evans function involves two parameters of the numerical method: h (step of the Runge-Kutta method) and L (half-width of the computational domain). Besides, it requires a profile of the soliton u(x). Let us fix parameters of the model µ = 1.615 and V 0 = 6 and perform the accuracy tests for the single-hump FGS which exists with these parameters [see Fig. 3(a)]. In order to find the solution u(x), we use the shooting method which requires integration of the equation (3). The Runge-Kutta method with the accuracy O(h 4 ) was used. In order to compute the Evans function, we solve the ODE Eq. (A9) using the Runge-Kutta method with the same accuracy. In order to check the accuracy of evaluation of the Evans function, we use the fact that λ = 0 is a double eigenvalue of Eq. (A1), and therefore In Table I  In order to locate the zeros of the Evans function, we use the Newton method. The derivative of the Evans function with respect to λ was approximated numerically using the first difference. As an example, in Table II we show the convergence to a simple zero computed for the same solution but with different h and L. For sufficiently small h the order of convergence is close to O(h 4 ).
While the Newton method is suitable for detecting and tracing of unstable eigenvalues, it, strictly speaking, does not ensure that all unstable eigenvalues have been found; we therefore cannot give a "computational proof" of stability; we can only conjecture on stability of a solution if no unstable eigenvalues has been found.
Most of the stability results presented above in Fig. 3  and 4 have been obtained with typical values h = π/256, L = 10π [for computation of solitons and the Evans function in the first spectral gap] and L = 20π [for computation of solitons and the Evans function in the second gap. We have also selectively tested several particular unstable solitons by repeating the computation with smaller step h and confirmed that the found instabilities are reproduced.
Finally, we notice that the detection of OIs is indeed a complex problem, and an advanced numerical tool is necessary for its solution. To this end, let us compare the results of our analysis with the results obtained with the Fourier collocation method (FCM) [19] which is another common tool for study of spectral stability of gap solitons. The FCM is known to provide excellent results for detection of relatively strong and exponentials instabilities, but as follows from our numerical experiments it may not be sufficiently robust to trace weak OIs. To illustrate this, let us consider FGSs in the first gap with V 0 = 6 fixing the parameters of the numerical method as L = 15π, h = π/4096 ≈ 10 −3 . We consider three close values of the chemical potential µ, namely,   On the other hand, the Evans function indicates a pair of OIs for each of the chosen values of µ. Namely, for µ = 2.22 the Evans function finds unstable eigenvalues λ 1 = 0.000299 + 4.924631i, λ 2 = 0.000566 + 4.925482i, for µ = 2.23 it gives λ 1 = 0.000298 + 4.931506i, λ 2 = 0.000601 + 4.932380i, and at µ = 2.25 it gives λ 1 = 0.000294 + 4.944917i, λ 2 = 0.000678 + 4.945838i [see Fig. 3(e)]. Moreover, these instabilities are reproduced with good accuracy if one repeats the numerical procedure with different L and h. Thus the Evans function robustly traces the instability.