Oscillatory instabilities of gap solitons in a repulsive BoseEinstein condensate
Abstract
The paper is devoted to numerical study of stability of nonlinear localized modes (“gap solitons”) for the spatially onedimensional GrossPitaevskii 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 (nonfundamental) 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),
(1) 
describes the meanfield dynamics of a quasionedimensional BoseEinstein condensate (BEC) confined in the potential PS03 (). In Eq. (1), is the complexvalued macroscopic wavefunction of the condensate. The squared amplitude of the wavefunction describes the local density of the BEC, while the gradient describes the velocity of atoms of condensate. The nonlinear term takes into account the interactions between the particles. The case corresponds to repulsive interparticle interactions, while describes attractive interactions. Both these cases are of physical relevance, but in what follows, we mainly focus on the repulsive case, . It is assumed that the potential is periodic, which corresponds to the optical confinement of the BEC BK04 (); MO06 ().
An important class of solutions of the GPE (1) corresponds to the stationary modes which can be represented in the form , where is a real parameter having the meaning of the chemical potential of the BEC. Function satisfies the conditions of the spatial localization
(2) 
Without loss of generality one can assume that is realvalued AKS02 (). Then can be found from the stationary GPE
(3) 
If the potential is periodic, the nonlinear modes satisfying (2)–(3) are called gap solitons BK04 (); MO06 (); EK2003 (); HOM02 (); Kiv2003 (); PSK04 (), since values of corresponding to these solutions lie in the spectral gaps of the linear Schrödinger equation BS91 ()
(4) 
The simplest class of the gap solitons are fundamental gap solitons (FGSs) EK2003 (); Kiv2003 (); Zhang09_1 (); Zhang09_2 (); MM2006 (). 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 . A variety of more complex (or higherorder) solitons includes truncated Bloch waves TrunkW09 (); Kiv2006 () (which consist of several inphase peaks placed in a row), various asymmetric states, and complex bound states of two (or more) wellseparated waves AHY12 (), 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 AA13 (). Specifically, if the lattice is deep enough, then all the gap solitons in the first gap can be put into a onetoone correspondence with the set of biinfinite sequences of symbols from a threesymbol 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 . For instance, the truncated Bloch waves TrunkW09 (); Kiv2006 () consisting of several inphase peaks placed in a row can be viewed as complexes of singlehump 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 HOM02 (); Kiv2003 (); PSK04 (); Pelinovsky11 (); Zhang09_2 (); MM2006 (); TrunkW09 (); Yang10 (); HAY11 (). The most relevant for our study outcomes for the case of repulsive interactions () can be summarized as follows. Significant part of the studies concluded that the singlehump FGSs are stable in the first gap HOM02 (); Kiv2003 (); Zhang09_1 (); Zhang09_2 (); MM2006 () and in the second gap Zhang09_1 (); Zhang09_2 (); MM2006 (). Regarding the higherorder states consisting of two or three inphase peaks, they have been reported unstable near the upper band edge in Kiv2003 () in the first gap. However, these states have been found to be stable both in the first Zhang09_2 (); TrunkW09 () and in the second Zhang09_2 () 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 PSK04 (); Yang10 (), the accurate tracing of OIs is much more challenging Kiv2003 (); PSK04 (). 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 onehump gap solitons in GPE is especially remarkable in view of wellknown OIs of Bragg gap solitons in nonlinear Dirac equations BarPelZem98 (); BarZem00 (); Derk05 (). The latter system can be deduced from the GPE with a shallow periodic potential using an asymptotic multiplescale expansion Pelinovsky11 (); SM04 (), and therefore the results about OIs of the solitons for Dirac system seem to be not consistent with the stability results for the singlehump 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 PSK04 (); Yang10 (). 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 PSK04 () to trace OIs of gap solitons in the attractive condensates. It was further demonstrated in Blank () 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 Derk05 (); Blank (). 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 Derk05 (); Blank ().
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 higherorder gap solitons in the repulsive BEC. We focus on the first and second spectral gaps. In each gap, we consider the singlehump FGS and two higherorder solitons bearing two or three inphase humps. The main outcomes of our numerical study can be outlined as follows.

In the first gap, all the considered solitons (including the singlehump 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.

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 in the first and in the second gap. On the other hand, our results advance the current understanding of the higherorder modes TrunkW09 (), 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.
Ii Families of gap solitons
In our study, we use a prototypical example of the periodic potential in the form
(5) 
where real characterizes the depth of the lattice.
The spectrum of the linear eigenvalue problem (4) with the potential (5) consists of one semiinfinite gap and a countable set of finite gaps BS91 () (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 semiinfinite 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 (the singlehump FGS), and two higherorder solitons, with two and three dominating peaks. In the terminology of TrunkW09 (), these higherorder solitons correspond to the truncated Bloch waves, and in terms of the coding approach AA13 () they are the bound states of two or three inphase singlehump FGSs with codes and . Equation (3) also supports families of outofphase multihump solitons (with the codes , , etc. However, these solutions have been reported to suffer relatively strong exponential instabilities Yang10 () 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 Kiv2003 (); PSK04 (); Zhang09_2 (); MM2006 (); TrunkW09 (); Yang10 (). These families can be visualized in the plane vs. , where the squared norm has the physical meaning of the number of particles in the condensate. The dependencies vs. for are presented in Fig. 3(a). In the first gap, one of the families (corresponding to the singlehump 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 Yang10 (); AHY12 (). 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 multihump 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 Zhang09_2 () (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) PSK04 (); HAY11 (). One of them is exponentially unstable PSK04 (), and another one (called subfundamental soliton in MM2006 ()) is unstable as the number of particles exceeds a certain threshold Zhang09_2 (). Since instabilities of these solutions are wellknown, in what follows we do not incorporate them in our study and focus on the simplest FGS and solitons consisting of two and three inphase FGSs whose instabilities have not been seen before.
Iii Stability of gap solitons
iii.1 Statement of the problem and the numerical method
Following to the standard procedure of the linear stability analysis, we consider a perturbed solution
(6) 
where describes a profile of a gap soliton, while and , , describe real and imaginary parts of a smallamplitude perturbation. After substitution of (6) into the GPE (1) and neglecting higherorder in and terms, one arrives at the following linear eigenvalue problem
(7) 
where linear operators are defined as
and . If the spectrum of is purely imaginary, then the amplitudes of perturbations and 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 . The corresponding eigenvector and the generalized eigenvector read
(8) 
Then and . In (8), 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 , 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 PSK04 ().
In spite of recent advances in the rigorous theory for analysis of the eigenvalue problem (7) Pelinovsky11 (); KP13 (), 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 Kiv2003 (); PSK04 (). One of wellelaborated numerical approaches to this task is based on the the Evans function PSK04 (); KP13 (); Derk05 (); Blank () which is especially useful if the eigenfunctions associated with unstable eigenvalues of (7) are poorly localized. The Evans function is an analytic function defined in the complex domain except the points of the essential spectrum of the linear stability operator . 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 KP13 (), 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 Blank () 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 Derk05 (); Blank (). Redefined in this way, the Evans function was used to study the stability of the surface gap solitons Blank (). 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 Blank (), for the sake of selfcontainment 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 wellelaborated tool for computing the instabilities of gap solitons Yang10 (). 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.
iii.2 Linear stability results
Using the Evans function approach, we have examined in details stability of fundamental (singlehump) and multihump solitons in the first and the second spectral gaps of a lattice (5) with the fixed depth . 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 () and a deeper () 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 . 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 eigenvalues. 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 [Fig. 3(b)]. The other two coexist in a sufficiently wide instability window which starts at at continues up to the upper gap edge. In the interval 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 do not exceed and typically lie between and .
Turning to multihump solitons in the first gap [Fig. 3(f,g)], we observe that the stability picture is remarkably similar to that for the FGSs: no unstable eigenvalues is observed for sufficiently small , but for large values of the chemical potential the solutions are unstable due to three quartets of unstable eigenvalues which cause two windows of instability. For the family of twohump solitons, the onset of instability is detected at . The instability persists in the narrow interval [see Fig. 3(c)]. No unstable eigenvalues is found for , but the instability again takes place for all between and the upper gap edge. In a similar way, the threehump solitons have two windows of instability: the narrow window is situated at [Fig. 3(d)], and the wide interval starts at and continuous up to gap edge. The typical instability rates of the multihump solitons are comparable with those for the fundamental solitons and do not exceed .
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 onehump 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 . However, the instabilities become stronger for multihump solutions in the second gap [Fig. 3(i, k)]: for these solutions reaches . 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 () and a deeper () 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 . For the fundamental family the first window of instability is situated at . In this region, the spectrum includes two quartets of unstable eigenvalues. The second window of instability starts at and continues to the upper gap edge. In this case, the instability is again associated with two quartets of eigenvalues. For twohump and threehump 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, [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 . 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 and in the narrow window . Thus the increase of enhances the stability of FGSs in the second gap.
iii.3 Solution of the timedependent 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 semiimplicit finite difference scheme TroPes2009 (). The equation was spatially discretized in a sufficiently large computational window with the grid step equal to . We have additionally introduced absorbing (dissipative) boundary layers by adding nonzero imaginary part to the potential for and . 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 multiplicative perturbation: . The temporal step was .
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 wellvisible in the plot with , 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 as in Fig. 4(b)] with an unstable FGS with the same but different . Again, the striking difference in their behaviors is much better pronounced on the plot which shows the density at [Fig. 7(c)]. Looking at the longterm 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 smallamplitude oscillations around some localized states. These oscillations persist for indefinitely long time.
The results for dynamics of multihump 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 BoseEinstein 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 solitons. 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 , but displays windows of stability for .
Besides the fundamental singlehump solitons, we have also examined stability of twohump and threehump 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 GrossPitaevskii 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 AA13 ()) and the spectrum of the corresponding stability problem. For the solitons from the first bandgap, the possibility of this relation has been discussed in Yang10 (), 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 FS14 () 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 ABK04 () 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 Blank (), we rewrite the linear stability eigenvalue problem (7) as follows
(9) 
where and . Further, Eq. (9) can be rewritten as an ODE system, where plays the role of a parameter:
(10) 
Due to the localization , in the limit Eq. (10) decouples into a pair of Mathieu equations
(11)  
(12) 
Equation (11) has the stable manifold for and the unstable manifold . The stable and unstable manifolds can be found numerically using the monodromy matrix associated with Eq. (11). In a similar way, one can compute the monodromy matrix associated with Eq. (12) and identify its stable and unstable manifolds. Using vectors and , we can find an unstable manifold of Eq. (10) and its stable manifold :
(13)  
(14) 
Using (13)–(14) the Evans function is defined as the determinant of a matrix Blank ()
(15) 
The definition (15) 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 , , solve the ODE system from to and from to in order to obtain the necessary values at , 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 (10) in the exterior algebra. Following to Blank (), we introduce a matrix function
(16) 
and instead of solving system (10), we introduce a new ODE system
(17) 
We compute two solutions of (17), and , which are correspond to initial conditions
Finally, the Evans function is redefined as
(18) 
where
(19) 
and stands for the standard inner product in with the complex conjugation of the first argument.
\backslashbox  

\backslashbox  

Appendix B Details of numerical implementation and accuracy issues
The evaluation of the Evans function involves two parameters of the numerical method: (step of the RungeKutta method) and (halfwidth of the computational domain). Besides, it requires a profile of the soliton . Let us fix parameters of the model and 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 , we use the shooting method which requires integration of the equation (3). The RungeKutta method with the accuracy was used. In order to compute the Evans function, we solve the ODE Eq. (17) using the RungeKutta method with the same accuracy. In order to check the accuracy of evaluation of the Evans function, we use the fact that is a double eigenvalue of Eq. (9), and therefore
(20) 
In Table 1 we show the numerical value of the Evans function at . It converges to the exact value as , which corresponds to the accuracy of the RungeKutta method.
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 2 we show the convergence to a simple zero computed for the same solution but with different and . For sufficiently small the order of convergence is close to .
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 , [for computation of solitons and the Evans function in the first spectral gap] and [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 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) Yang10 () 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 fixing the parameters of the numerical method as , . We consider three close values of the chemical potential , namely, and . Using the FCM (implemented with spectral harmonics using the computer codes from Yang10 ()) we have detected a pair of OIs ( and ) at . However, the FCM does not indicate any instability for and . Moreover, if we try to check the obtained results by repeating the computation in a smaller () or in a larger () computational window, the FCM does not find any unstable eigenvalue at all (for each of the three considered values of ). In a similar way, the results of the FCM are sensitive to the choice of the spatial step . Hence, in this particular situation the results of the FCM cannot be considered as reliable.
On the other hand, the Evans function indicates a pair of OIs for each of the chosen values of . Namely, for the Evans function finds unstable eigenvalues , , for it gives , , and at it gives , [see Fig. 3(e)]. Moreover, these instabilities are reproduced with good accuracy if one repeats the numerical procedure with different and . Thus the Evans function robustly traces the instability.
References
 (1) L. Pitaevskii, S. Stringari, BoseEinstein Condensation, Clarendon Press, Oxford, 2003.
 (2) V. A. Brazhnyi, V. V. Konotop, Theory of Nonlinear Matter Waves in Optical Lattices, Mod. Phys. Lett. B 18, 627651 (2004).
 (3) O. Morsch and M. Oberthaler, Dynamics of BoseEinstein condensates in optical lattices, Rev. Mod. Phys. 78, 179 (2006).
 (4) G. L. Alfimov, V. V. Konotop and M. Salerno, Matter solitons in BoseEinstein Condensates with optical lattices, Europhys. Lett. 58 (2002) 7–13.
 (5) N.K. Efremidis and D.N. Christodoulides, Lattice solitons in BoseEinstein condensates, Phys. Rev. A 67, (2003) 063608.
 (6) K. M. Hilligsoe, M. K. Oberthaler and K. P. Marzlin, Stability of gap solitons in a BoseEinstein condensate, Phys. Rev. A 66, 063605, (2002).
 (7) P. J. Y. Louis, E.A. Ostrovskaya, C.M. Savage and Yu.S. Kivshar, BoseEinstein condensates in optical lattices: Bandgap structure and solitons, Phys. Rev. A 67 (2003) 013602.
 (8) D.E. Pelinovsky, A.A. Sukhorukov, and Yu.S. Kivshar, Bifurcations and stability of gap solitons in periodic potentials, Phys. Rev. E, 70 (2004), 036618.
 (9) F.A. Berezin and M.A. Shubin, The Shrödinder Equation, Kluwer, Dordrech, 1991.
 (10) Yo. Zhang and B. Wu, Composition Relation between Gap Solitons and Bloch Waves in Nonlinear Periodic Systems, Phys.Rev.Lett. 102 (2009) 093905.
 (11) Yo. Zhang, Zh. Liang and B. Wu, Gap solitons and Bloch waves in nonlinear periodic systems, Phys. Rev. A 80 (2009), 063815.
 (12) T. Mayteevarunyoo, and B.A. Malomed, Stability limits for gap solitons in a BoseEinstein condensate trapped in a timemodulated optical lattice, Phys. Rev. A 74 (2006) 033616.
 (13) J. Wang, J. Yang, T.J. Alexander, and Yu.S. Kivshar, TruncatedBlochwave solitons in optical lattices, Phys. Rev. A 79 (2009) 043610.
 (14) T.J. Alexander, E.A. Ostrovskaya and Yu.S. Kivshar, Selftrapped nonlinear matter waves in periodic potentials, Phys. Rev. Lett. 96 (2006) 040401.
 (15) T.R. Akylas, G. Hwang, and J. Yang, From nonlocal gap solitary waves to bound states in periodic media, Proc. R. Soc. A 468 (2012) 116135.
 (16) G.L. Alfimov and A.I. Avramenko, Coding of nonlinear states for the GrossPitaevskii equation with periodic potential, Physica D 254 (2013) 2945.
 (17) D.E. Pelinovsky, Localization in periodic potentials: from Schrödinger operators to the GrossPitaevskii equation, University Press, Cambridge, 2011.
 (18) H. Sakaguchi and B.A. Malomed, Dynamics of positive and negativemass solitons in optical lattices and inverted traps J. Phys. B: At. Mol. Opt. Phys. 37 (2004) 1443.
 (19) J. Yang, Nonlinear Waves in Integrable and Nonintegrable Systems, SIAM, Philadelphia, 2010.
 (20) G. Hwang, T. R. Akylas, J. Yang, Gap solitons and their linear stability in onedimensional periodic media, Physica D, 240 (2011) 10551068.
 (21) I.V. Barashenkov, D.E. Pelinovsky, and E.V. Zemlyanaya, Vibrations and Oscillatory Instabilities of Gap Solitons, Phys. Rev. Lett. 80 (1998) 51175120.
 (22) I.V. Barashenkov and E.V. Zemlyanaya, Oscillatory instabilities of gap solitons: a numerical study, Comp. Phys. Comm. 126 (2000) 2227.
 (23) G. Derks, G.A. Gottwald, A robust numerical method to study oscillatory instability of gap solitary waves, SIAM J. Appl. Dyn. Sys. 4 (2005) 140158.
 (24) E. Blank and T. Dohnal, Families of Surface Gap Solitons and Their Stability via the Numerical Evans Function Method, SIAM J. Appl. Dyn. Sys. 10(2) (2011) 667706.
 (25) T. Kapitula and K. Promislow, Spectral and Dynamical Stability of Nonlinear Waves, Applied Mathematical Sciences 185, Springer Science+Business Media New York 2013.
 (26) V.A. Trofimov, N.V. Peskov, Mathematical Modelling and Analysis 14 (2009) 109126.
 (27) R. Fukuizumi, A. Sacchetti, Stationary States for Nonlinear Schrödinger Equations with Periodic Potentials, J. Stat. Phys (2014) 156, 707–738.
 (28) G. L. Alfimov, V. A. Brazhnyi, V. V. Konotop, On classification of intrinsic localized modes for the discrete nonlinear Schrödinger equation, Physica D 194 (2004) 127–150.