Gap solitons for the repulsive Gross-Pitaevskii equation with periodic potential: coding and method for computation

The paper is devoted to nonlinear localized modes ("gap solitons") for the spatially one-dimensional Gross-Pitaevskii equation (1D GPE) with a periodic potential and repulsive interparticle interactions. It has been recently shown (G. L. Alfimov, A. I. Avramenko, Physica D, 254, 29 (2013)) that under certain conditions all the stationary modes for the 1D GPE can be coded by bi-infinite sequences of symbols of some finite alphabet (called"codes"of the solutions). We present and justify a numerical method which allows to reconstruct the profile of a localized mode by its code. As an example, the method is applied to compute the profiles of gap solitons for 1D GPE with a cosine potential.


1.
Introduction. The Gross-Pitaevskii equation (GPE) is a commonly recognized model for description of nonlinear matter waves in Bose-Einstein condensates (BECs) [20]. In the dimensionless form spatially one-dimensional (1D) GPE, describes the dynamics of a cigar-shaped BEC cloud in the so-called mean-field approximation. In (1) the variables are properly scaled, Ψ = Ψ(t, x) is the complexvalued macroscopic wavefunction of the condensate, and its squared amplitude |Ψ(t, x)| 2 describes the local density of BEC. The nonlinear term σ|Ψ| 2 Ψ takes into account the interactions between the particles: the case σ = 1 corresponds to the repulsive interparticle interactions, σ = −1 corresponds to the attractive interactions. Both these cases are of physical relevance. The real-valued potential U (x) describes a trap which confines the BEC. In the case of optical confinement of BEC, the potential U (x) is a periodic function. An important particular class of solutions of the GPE corresponds to the stationary modes. These modes are of the form Ψ(t, x) = e −iωt ψ(x) and the wavefunction ψ(x) solves the stationary GPE ψ xx + (ω − U (x))ψ − σ|ψ| 2 ψ = 0.
(2) Figure 1. The diagram of gaps and bands of the Matheiu equation ψ xx + (ω − A cos 2x)ψ = 0. The regions where according to [3] the coding of bounded solutions for the nonlinear equation (4) is possible are shown in the first and the second gaps: if ω and A belong to Region 1, then the bounded solutions can be coded using an alphabet of N = 3 symbols; in Region 2 an alphabet of N = 5 symbols is necessary.
If the potential U (x) is periodic, the modes satisfying (3) and (4) are called gap solitons. This terminology can be understood if one considers the linear version of (4), i.e. ψ xx + (ω − U (x))ψ = 0. (5) Due to the conditions (3), (5) approximates (4) for |x| ≫ 1. Equation (5) is a linear eigenvalue problem for the Schrödinger operator with periodic potential U (x) and eigenvalue ω. It is well-known [11], that the spectrum of such a problem typically consists of a countable set of continuous bands separated by gaps. Zero boundary conditions at x → +∞ and x → −∞ for (5) cannot be satisfied if ω lies in a band of the spectrum. Therefore the solutions of (4) satisfying (3) can exist only if ω belongs to one of the gaps.
It is known that apart from the gap solitons [2,15,16,19], (4) supports nonlinear periodic structures (nonlinear Bloch waves) [16,21], domain walls [13], gap waves [1] and so on. The authors of [23,24] observed that more complex structures described by (4), such as nonlinear Bloch waves, can be built from elementary entities called fundamental gap solitons (FGS) by means of the so-called composition relation. This idea was further developed in [3], where it was suggested to code complex nonlinear modes of (4) using the methods of symbolic dynamics. The approach of [3] is consistent with the recent results of [12] where it was proved that in the semiclassical limit a certain class of gap solitons for (4) can be mapped to a set of localized modes of Discrete Nonlinear Schrödinger Equation (DNLS) which admit the description in terms of coding sequences [4].
The "coding" approach of [3] exploits the following fact. In the case of repulsive interactions (σ = 1) the "most" of the solutions of (4) collapse (i.e., tend to infinity) at some finite point of the real axis. The complementary set of non-collapsing solutions (including all gap soliton solutions) can be associated with a set on the plane of initial conditions for (4). This set has fractal structure and under certain conditions (Hypotheses 1-3 of [3]) it can be described in terms of symbolic dynamics and mapped into the set of bi-infinite sequences of symbols of some finite alphabet A. The relation between the solutions and the codes is a homeomorphism: each bounded in R solution of (4) corresponds to unique bi-infinite code of the type and each code of this type corresponds to unique bounded in R solution of (4). It was argued in [3], that in the case of the cosine potential U (x) = A cos 2x (6) this coding is possible for vast areas in the parameter plane (ω, A). These areas are situated in the gaps of linear equation (5) which becomes the Mathieu equation for the cosine potential (6) (see Figure 1). For Region 1 (situated in the first gap), the alphabet A consists of three symbols, while in the Region 2 (situated in the second gap of the Mathieu equation) the alphabet of five symbols must be used. In order to transform the coding approach into a convenient tool for analysis and computation of gap solitons, one has to answer the following practical questions: Question 1. How can one get a code of a given gap soliton? Question 2. How can one compute a profile of gap soliton by its code?
The answer to the Question 1 is relatively simple (see subsection 4.1) and follows immediately from the results of [3]. The answer to the Question 2 is not simple, and it is the focal point of this paper. From the practical point of view, the answer to this question offers a method for numerical computation of gap solitons. While several numerical approaches to this task have been described in the literature (see, for instance, a survey in the monograph [22]), the peculiar advantage of a codingbased numerical algorithm is that it yields a nonlinear mode which corresponds to the specific code chosen beforehand. Applying the algorithm to different codes it is possible to compute and classify gap solitons with shapes of different complexity assigned in advance.
Looking ahead, we note that a high-accuracy algorithm for computation of nonlinear modes is a crucial ingredient of the numerically accurate stability analysis. It is known that even the simplest FGSs from lower gaps undergo oscillatory instabilities under both attractive [19] and repulsive [14] nonlinearities (instabilities of this type have been also found in the coupled mode limit [9,10,17]). The oscillatory instabilities of gap solitons can be extremely weak and can be revealed by a quite thorough study only [14]. This, in particular, requires to compute the profile of nonlinear mode with high accuracy.
In what follows, the consideration is restricted by the repulsive version of GPE (σ = 1). Section 2 contains a brief description of the coding approach. Section 3 includes some theoretical results necessary for construction of the algorithm. Section 4 contains the description of the main numerical algorithm whose listing is presented in Appendix B. In section 5, the method is illustrated for the GPE with cosine potential (6). Section 6 concludes the paper and presents a short outlook for the future work.
2. Review of the coding approach. Let us briefly describe the method for coding of nonlinear modes for (4) (see [3] for the detailed presentation). The basis of this method, the coding theorem, is formulated in subsection 2.3. In order to formulate the coding theorem, in subsection 2.1 and 2.2 we introduce several definitions, assumptions, and auxiliary concepts.
2.1. Collapsing solutions. For σ = 1 the equation (4) for the nonlinear modes reads The potential U (x) is assumed to be bounded, π-periodic and even. A prototypical example is the cosine potential (6). The following definitions are necessary, [3].
Evidently, ∆ 0 consists of the points that have both T -image and T -pre-image. The features of ∆ 0 are as follows: 1. ∆ 0 is bounded (Theorem 2.1 in [3]) and includes the origin (0, 0). Definition 2.6. The sets ∆ ± n , n = 0, 1, . . . are defined by the following recurrence rule Evidently, if p ∈ ∆ + n then there exist Repeating the procedure n times one arrives at the relations Let us illustrate the construction of the sets ∆ ± n by an example. Figure 3 shows the action of T on ∆ 0 for the case of potential U (x) = −3 cos 2x and ω = 1. The set GAP SOLITONS: CODING AND METHOD FOR COMPUTATION 7 ∆ 0 consists of three connected components, marked D −1 , D 0 and D 1 . T -images of these components are three infinite strips. Each of these infinite strips intersects each of the connected components D −1,0,1 , and therefore the set ∆ − 1 = T ∆ 0 ∩ ∆ 0 consists of nine connected components. Further, the set ∆ − 2 consists of 27 connected components and so on. The sets ∆ + n are symmetric to ∆ − n with respect to ψ-axis, so they have the same properties.

Assumptions.
Recall that a function f (x) is called Lipschitz function if there exists finite κ > 0 such that for any x 1 and In order to formulate the main result of [3] let us introduce the following definitions.

Definition 2.8 (An island). We call an island an open curvilinear quadrangle
D ⊂ R 2 bounded by two pairs of curve segments, such that the segments from one pair do not intersect and are increasing, and the segments from another pair do not intersect and are decreasing.
Let us make the following assumptions about the set ∆ 0 and the map T associated with (7).
We say that one boundary of an island is opposite to another boundary if these two boundaries do not intersect.
Hypothesis+. For each island D i from Hypothesis 1, one boundary of D i consists of π-collapsing to +∞ points, and the opposite boundary consists of π-collapsing to −∞ points. Similarly, from another pair of opposite boundaries, one consists of −π-collapsing to +∞ points, and another consists of −π-collapsing to −∞ points.
In what follows, for an island D i , we will denote α ± i the boundaries consisting of π-collapsing to ±∞ points and denote β ± i the boundaries consisting of −π-collapsing to ±∞ points, see Figure 4.
Hypothesis+ naturally holds when U ± π are infinite curvilinear strips (see Figure 2 as an example). In this case the boundary ∂U + π is formed by two curves, one of them consists of π-collapsing to +∞ points and other one consists of π-collapsing to −∞ points. These curves can be mapped one into another by reflection with respect to the origin. The similar situation takes place for the boundary ∂U − π . The opposite sides α ± of each island are formed by segments of opposite boundaries of ∂U + π , as well as opposite sides β ± are formed by segments of opposite boundaries of ∂U − π . Definition 2.9 (v-and h-curves). Let D be one of the islands from Hypothesis 1 bounded by curve segments α ± and β ± . We call v-curve a Lipschitz curve β with endpoints on α − and α + which is increasing if β ± are increasing and is decreasing if β ± are decreasing. Similarly, we call h-curve a Lipschitz curve α with endpoints on β − and β + which is increasing if α ± are increasing and is decreasing if α ± are decreasing.  Hypothesis 3. The sets ∆ ± n defined by (9), (10) are such that lim For any particular choice of potential U (x) and parameter ω the validity of Hypotheses 1-3 should be supported by numerical arguments. On the basis of a systematic numerical investigation, the following conjecture was formulated in [3]: Conjecture. If the parameters ω and A for (7) with the cosine potential (6) belong to dark gray areas (Regions 1 and 2) on the parameter plane (ω, A) in Figure 1, then Hypotheses 1-3 hold.
The upper boundaries of Regions 1 and 2 coincide with upper boundary of the first and the second gap, respectively. The lower boundaries of these regions correspond to breaking of Hypothesis 1: if ω and A lie below these boundaries, then ∆ 0 is not split into a set of disjoint islands, i.e., Hypothesis 1 does not hold. Hypothe-sis+ was not introduced in [3]. However for ω and A lying in Regions 1 and 2 the sets U ± π are curvilinear strips bounded by two curves, consisting of ±π-collapsing to +∞ points and of ±π-collapsing to −∞ points correspondingly. This implies that Hypothesis+ also holds for ω and A lying in Regions 1 and 2.
The set Ω N has the structure of a topological space where the neighborhood of . . Definition 2.12. Denote by P the set of bi-infinite sequences (called orbits) where each p n = (ψ n , ψ ′ n ), n = 0, ±1, ±2, . . ., belongs to ∆ 0 . It is clear that any orbit is uniquely determined by its single entry, for instance, p 0 . The point p 0 cannot be arbitrary since, as it was mentioned before, a great part of points of the plane of initial data are collapsing, either forward or backward (or both), and a collapsing point cannot generate a bi-infinite sequence s. The set P has the structure of a metric space where the distance ρ between the elements Definition 2.13. Define a map Σ : P → Ω N as follows: i k is the number i of the component D i where the point p k lies.
The following statement was proved in [3]: Theorem 2.14. Assume that Hypotheses 1-3 hold. Then Σ is a homeomorphism between the topological spaces P and Ω N .
In the rest of the paper we assume that Hypotheses 1-3 and Hypothesis+ hold. 3.1. Indexing of h-strips in an island. Consider the action of the Poincare map T on an island D i . The points that are sent by T from the island D i to an island By analogy, consider points that visit the islands D i0 , D i1 , . . . , D in under the action of T in the given order. Any point p in D i0 whose orbit visits the islands in this order belongs to an h-strip H i0...in which can be constructed using the following recurrence rule In other words, the h-strip H i0...in coincides with the set of points p ∈ R 2 satisfying three conditions: The points (a)-(c) imply that the h-strips with different multi-indices are embedded in the following sense By definition, h-strip H i0...in is bounded by two h-curves. Points from one of them are (n + 1)π-collapsing to +∞. Denote this h-curve by α + i0...in . By construction, Points from another h-curve are (n + 1)π-collapsing to −∞. Denote this h-curve by α − i0...in and note that T n α − i0...in ⊂ α − in . The indexing of v-strips is introduced in a similar manner. The points that are sent by T −1 from an island D i to an island D j constitute a v-strip denoted as By analogy, v-strip denoted as V i−n...i0 is the set of points which belong to D i0 and under the action of T −1 visited the islands D i0 , D i−1 , . . ., D i−n in the given order.
This v-strip can be constructed using the recurrence rule The v-strips with different multi-indices are embedded according to the following rule If a point p ∈ V i−n...i0 generates an orbit, then this orbit contains a block {. . . , i −n , i −n+1 , . . . , i 0 , . . .}.

3.2.
The fixed points of T and their stable and unstable manifolds. The following statement is valid.  (7) is invariant with respect to the change ψ(x) → −ψ(x), the v-curve V 0 and h-curve H 0 are symmetric with respect to the origin.

Remark 2.
Let DT (O i * ) be the operator of linearization for T at fixed point O i * ∈ D i * . Assume that the eigenvalues λ 1 and λ 2 of DT (O i * ) are such that λ n 1,2 = 1, n = 1, 2, 3, 4 (see [8], Addition 8). Since T is an area-preserving map, λ 1 λ 2 = 1. Then, due to existence of local stable and unstable manifolds, H i * and V i * , the fixed point O i * is of hyperbolic type. Assume that |λ 1 | > 1, |λ 2 | < 1. If the corresponding eigenvectors of λ 1 and λ 2 are denoted by e 1 and e 2 , then the tangent vector to V i * is e 1 and the tangent vector to H i * is e 2 .
Each ∞-curve passes through each island exactly once. Let us introduce the order of island with respect to an ∞-curve.
The following lemma is valid.
Lemma 3.6. All ∞-curves pass through the islands in ∆ 0 in the same order.
Proof. We outline the proof as follows. First of all, let us observe that by simple geometric reasons, if two ∞-curves γ 1 = T β (1) and γ 2 = T β (2) (β (1) and β (2) are v-curves) do not intersect they pass the islands in the same order. This means that if β (1) and β (2) do not intersect then the statement of the lemma holds. If β (1) and β (2) intersect they belong to the same island D i * . Then in D i * there exists one more v-curve β (3) that intersects neither β (1) nor β (2) . By transitivity, one concludes that the statement of the lemma holds for all ∞-curves.
Due to Lemma 3.6, we can use an arbitrary ∞-curve to order the islands. Consider the ordering of island with respect to ∞-curve T V 0 . Since T V 0 is symmetric with respect to the origin, the origin (0, 0) belongs to the "central" island D 0 , and the islands D i and D −i , i = 1, 2, . . . , L are situated symmetrically with respect to the origin. Let us enumerate the islands according to the introduced order. So in what follows we assume that Definition 3.7. Let γ be an ∞-curve. Denote By construction, the points h ± i (γ) are π-collapsing to ±∞ points. Remark 4. By simple geometric arguments, one can prove that if h − i (γ) ≺ h + i (γ) for some ∞-curve γ, then the same relation h − i (γ 1 ) ≺ h + i (γ 1 ) holds for any other ∞-curve γ 1 ; vice versa if h + i (γ) ≺ h − i (γ) for some ∞-curve γ, then the same relation holds for any other ∞-curve.
Remark 5. One can also prove that for any two successive islands D i and D i+1 in (17) and for arbitrary ∞-curve γ one of the two alternatives holds: . Remark 6. In simple terms, if U + π and U − π are curvilinear strips (as in Figure 2), then the ∞-curves pass "along" the strip U − π , and the order (17) is the order in which the islands appear in U − π .      Remark 7. It is also straightforward to prove that if Hypotheses 1-3 and Hypoth-esis+ hold and the ordering of islands in ∆ 0 is given by (17) then the following implications are valid If the order of the islands and the positions of their boundaries α ± i are known (with respect to an arbitrary ∞-curve), then the established relations are sufficient to describe positions of h-strips with different indexes. The ordering scheme can be sketched as follows. It is necessary to draw the islands D −L , D −L+1 , . . ., D L from the left to the right, to draw the "global" arrow in the same direction (the upper bold arrow in Figure 5). Next, it is necessary to draw an arrow over each island such that any of those arrows is directed from α − i to α + i . According to Remark 5, the directions of any two successive arrows will be opposite. Then one can continue the procedure recursively, drawing arrows of higher levels and identifying the position of h-strips with different indexes following to directions of the arrows. In simple terms, the order of h-strips of (n+1)th level within an h-strip of nth level repeats the order of islands D i , i = −L, . . . , L within the set ∆ 0 . The peculiarity is that the last index of embedded h-strips increases in direction from α − i0...in to α + i0...in . An example of the ordering procedure is presented in Figure 5. As follows from the ordering procedure, the positions of h-strips with different indexes can be catalogued conveniently using a (2L + 1)-ary tree graph, where the root of the graph corresponds to the island D 0 and leaves of increasing levels denote h-strips with indexes of growing length. Figure 6 presents an example of such a tree for L = 1 (four levels of the tree are shown).

4.
Algorithm for constructing of a gap soliton by the given code.
This solution is a gap soliton, i.e., ψ(x) → 0 as x → ±∞. Hence for construction of the profile of a gap soliton ψ(x) by the code (18) one has to find any entry of the orbit s with sufficiently good accuracy. Then the profile of ψ(x) can be recovered from the numerical solution of the Cauchy problem (19).

4.2.
Description of the algorithm. Consider a code of the form (18). Having infinite number of zeros from the left, this code corresponds to an orbit that remains in D 0 for infinitely many iterations of T −1 . The relation between the entries of the code (18), the entries of the orbit s = {. . . p −1 , p 0 , p 1 , . . .}, p n = (ψ n , ψ ′ n ), and the values of the solution ψ(x) of (7) is given in Table 1. Note that for n ≤ 0 the points p n are situated on V 0 , the local unstable manifold of O 0 , and for n > m the points p n lie on H 0 , the local stable manifold of O 0 . Introduce Let us take a closer look at the island D 0 . It is limited by a pair of h-curves α + 0 and α − 0 and a pair of v-curves β + 0 and β − 0 (see Figure 4). The curves α ± 0 are such that if q = (ψ, ψ ′ ) ∈ α ± 0 , then the solution ψ(x) of (7) with these initial data at x = 0 satisfies the condition lim x→π ψ(x) = ±∞, (the sign + or − corresponds to the superscripts in α ± 0 ). According to Corollary 2, T V 0 is an ∞-curve that intersects α + 0 and α − 0 , and V 0 = T V 0 ∩ D 0 . According to Definition 3.7, denote , and so on. The algorithm exits with the value ε * = ε(p −M ) found with a certain (controlled) accuracy. Once ε * is found, then the shape of the soliton ψ(x) on the interval x ∈ (−M π; Kπ) is described by the solution of Cauchy problem for (7) with initial data ψ(−M π) = ε * ξ, ψ x (−M π) = ε * η. For x < −M π and x > Kπ the profile can be continued with linear asymptotics, i.e., by the solutions of (5) vanishing at −∞ and +∞, respectively, which are properly scaled to ensure the continuity of the whole solution and its derivative. 4.4. Some details. Algorithm 1 was implemented in Python using libraries numpy and scipy. It includes nested recursive procedure for processing of each index of the gap soliton code. The following details should be mentioned: 1. The Cauchy problem for (7) was solved using standard Runge-Kutta method of 4th order with constant step.
2. We found that the values M = 2 or 3 and K = 2 or 3 are enough for the most of the calculations. The accuracy was estimated by changing these parameters and comparing the results. In vicinity of the upper edge of the gaps (when the solitons are poorly localized) M and K should be increased.
3. The condition lim x→nπ ψ(x) = ±∞ that appears in the algorithm, was replaced by the condition ψ(nπ) = ψ ∞ where ψ ∞ is a large positive (or, correspondingly, great negative) number. Most of the calculations were fulfilled for |ψ ∞ | = 10 6 and checked for |ψ ∞ | = 10 8 to confirm that this choice of ψ ∞ does not affect the results.
4. The main difficulty in practical implementation of Algorithm 1 is caused by the fact that widths of h-strips H (0×M)i1i2... decrease rapidly as the number of indices grows. Therefore the difference eventually becomes beyond the computer accuracy. This hinders the calculation of solitons with large m if the program code operates with floating-point numbers with short mantissa. In order to overcome this difficulty, the Python library gmpy2 for arbitrary-precision arithmetic was used. In principle, this allows to compute gap solitons with codes with arbitrary m (the algorithm was tested for codes with m ≤ 10). 5. The algorithm does not allow to find gap solitons in the situation when Hypotheses 1-3 do not hold. In particular, it cannot be applied to compute smallamplitude gap solitons in cosine potential (6) with values (ω, A) situated close to the lower gap edges [see light-gray areas in (1)]. In this situation, a numerical continuation procedure in ω or in A can be applied.

5.
Results. The algorithm was applied to (7) with cosine potential (6). Recall that according to Conjecture made in subsection 2.   [24,23] (in [22] this solution is called "on-site" gap soliton). FGS 1 is the unique elementary entity for Region 1. The code {. . . , 0, −1, 0, . . .} corresponds to FGS 1 taken with minus sign (not shown in Figure 8). More complex gap solitons can be regarded as bound states of several FGS 1 taken with the plus or minus sign. For instance, the "off-site" gap soliton in the terminology of [22] has the code {. . . , 0, 1, −1, 0, . . .} and can be regarded as a Figure 8 . This is the so-called subfundamental soliton which was introduced in [18]; solution of this type was also discussed in [23] where it was termed to as the second fundamental gap soliton. The code {. . . , 0, 1, 0, . . .} corresponds to FGS 2 taken with minus sign. Again, more complex gap solitons can be regarded as bound states of FGS 1 and FGS 2 taken with different signs.
While the solitons illustrated in Figure 8 and Figure 9 correspond to relatively simple solutions with short codes (m = 1, 2 or 3), the algorithm is also applicable to more complex solitons and to solutions from higher gaps. For instance, in Figure 10 complex gap solitons, (m = 10) from the 2-nd and the 3-rd gaps are depicted. 6. Summary and discussion. In this paper we have proposed a numerical algorithm 1, for computation of localized modes for the spatially one-dimensional Gross-Pitaevskii equation with a periodic potential U (x). These modes, i.e., gap solitons, are of the form Ψ(x, t) = e −iωt ψ(x) where ψ(x) satisfies (7) and localization conditions (3). It is assumed that the parameters of the potential U (x) allow for coding of all these modes by bi-infinite sequences of symbols from some alphabet of finite length. It is known [3], that this coding is indeed possible, if certain assumptions (Hypotheses 1-3) hold, see Theorem 2.14. In particular, according to [3], the coding is possible for the potential U (x) = A cos 2x, if the parameters ω (the frequency of the mode) and A (the amplitude of the potential) belong to vast areas on the plane (ω, A), see Figure 1.
The coding approach exploits the fact that the "most" of solutions of (7) tend to infinity at some finite point of real axis. As a result the initial data for Cauchy problem for (7) corresponding to bounded in R solutions form a fractal subset of R 2 . The initial data corresponding to gap soliton solutions are situated in the intersection of this fractal set with the local unstable manifold V 0 of the zero fixed point O 0 = (0, 0). The numerical procedure for computation of gap solitons consists in iterative localization of the segments on V 0 where the initial data corresponding to a given code are situated. Then the profile of the nonlinear mode can be computed using the Runge-Kutta method. Algorithm 1 was implemented in Python and applied for gap solitons in cosine potential (6).
Algorithm 1 allows for various generalizations. With minor changes the algorithm can be applied in the case of more complex odd nonlinearities, both autonomous and non-autonomous. It is known that the structure of U ± π in (7) in cubic-quintic [6] and periodic [7] nonlinearities are similar to that described here, and therefore the application of Algorithm 1 in these cases is straightforward. For these problems, Algorithm 1 may be a useful tool to study bifurcations of gap solitons for (7) with  different codes, similar to the one fulfilled in [4] for the intrinsic localized modes of the DNLS .
Another possible application of Algorithm 1 is the study of stability of gap solitons. An interesting question is whether it is possible to predict stability or instability of a gap soliton by its code. The results of this study will be reported elsewhere.