{'fraction_non_alphanumeric': 0.046786551170322, 'fraction_numerical': 0.01827229404175772, 'mean_word_length': 4.5086294416243655, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From nonlocal gap solitary waves to bound states in periodic media

Solitary waves in one-dimensional periodic media are discussed employing the nonlinear Schr\"odinger equation with a spatially periodic potential as a model. This equation admits two families of gap solitons that bifurcate from the edges of Bloch bands in the linear wave spectrum. These fundamental solitons may be positioned only at specific locations relative to the potential; otherwise, they become nonlocal owing to the presence of growing tails of exponentially-small amplitude with respect to the wave peak amplitude. Here, by matching the tails of such nonlocal solitary waves, higher-order locally confined gap solitons, or bound states, are constructed. Details are worked out for bound states comprising two nonlocal solitary waves in the presence of a sinusoidal potential. A countable set of bound-state families, characterized by the separation distance of the two solitary waves, is found, and each family features three distinct solution branches that bifurcate near Bloch-band edges at small, but finite, amplitude. Power curves associated with these solution branches are computed asymptotically for large solitary-wave separation, and the theoretical predictions are consistent with numerical results.


Introduction
Nonlinear wave phenomena in periodic media are currently attracting a great deal of research interest in nonlinear optics, Bose-Einstein condensation and applied mathematics (see Christodoulides et al. 2003 & Yang 2009; Yang 2010 for reviews). Apart from scientific curiosity, this research activity is also driven by various potential applications, ranging from light routing in lattice networks (Christodoulides & Eugenieva 2001) to image transmission through nonlinear media .
A characteristic feature of periodic media is the existence of bands in the linear spectrum where linear disturbances, the so-called Bloch modes, may propagate. Between these Bloch bands are bandgaps in which linear disturbances are evanescent but nonlinear localized modes, commonly known as gap solitons, are possible. Since the first theoretical prediction (Christodoulides & Joseph 1988) and experimental observation (Eisenberg et al. 1998) of fundamental solitons in the semiinfinite bandgap of one-dimensional (1D) periodic waveguides, various types of gap solitons in one-and multi-dimensions have been reported theoretically and experimentally. Examples include 2D fundamental gap solitons, vortex solitons, dipole solitons, reduced-symmetry solitons, vortexarray solitons, truncated-Bloch-wave solitons and arbitrary-shape gap solitons (see Yang 2010 for a review). Multi-soliton bound states in periodic media have also been constructed in the framework of a discrete nonlinear Schrödinger (NLS) model (Kevrekidis et al. 2001). In addition, the stability of some of these gap solitons has been examined (Pelinovsky et al. 2004(Pelinovsky et al. , 2005Shi et al. 2008;Yang 2010;Hwang et al. 2011).
The plethora of gap solitons in periodic media calls for systematic identification and classification of the various types of solutions. On this issue, some progress has been made, particularly for gap solitons that bifurcate from linear Bloch modes at the edges of Bloch bands. Specifically, in 1D, only two gap-soliton families bifurcate from each edge of a Bloch band under self-focusing or selfdefocusing nonlinearity (Neshev et al. 2003;Pelinovsky et al. 2004;Hwang et al. 2011), and the positions of these solitons relative to the periodic medium (or potential) are determined by a certain recurrence relation (Hwang et al. 2011). In 2D, four gap-soliton families (or a multiple of four families) bifurcate from each edge of a Bloch band under self-focusing or self-defocusing nonlinearity (Yang 2010).
However, for the majority of gap solitons, that do not bifurcate from band edges, systematic theoretical treatment is lacking at present. Numerical results indicate that, typically, those solution families feature multiple branches which bifurcate near band edges at small (but finite) amplitude (see Yang 2010 and references therein), yet there has been no analytical explanation for this phenomenon.
In this article, we employ the NLS equation with a spatially periodic potential to analytically construct and classify a broad class of gap solitons that bifurcate away from band edges in 1D periodic media. Our approach is based on the observation that the two families of fundamental gap solitons that bifurcate from a band edge can be located only at specific positions relative to the underlying potential; otherwise, they would be nonlocal owing to the presence of growing tails of exponentially-small amplitude with respect to the peak wave amplitude. However, by piecing together such nonlocal solitary waves, it is possible to construct higher-order gap solitons, or bound states. The proposed asymptotic theory is illustrated by working out details for bound states involving two nonlocal solitary waves. We show that a countable set of bound-state families can be constructed. Each family is characterized by the separation distance of the two solitary waves involved, and features three distinct solution branches that bifurcate near band edges. The analytical predictions are verified by comparison with numerical results.

Preliminaries
Our study of nonlinear wave phenomena in periodic media is based on the 1D NLS equation, with a periodic potential V (x) and self-focusing (σ = 1) or self-defocusing (σ = −1) cubic nonlinearity. This equation is the appropriate mathematical model for Bose-Einstein condensates loaded in optical lattices (Dalfovo et al. 1999;Morsch & Oberthaler 2006) and laser beam transmission in photonic lattices under the paraxial approximation (Yang 2010). Although it is possible to consider a general periodic potential V (x) as in Hwang et al. (2011), here, for simplicity, we shall work with the sinusoidal potential which is π-periodic and also symmetric in x, V 0 being the potential depth. This potential frequently arises in nonlinear optics and Bose-Einstein condensates.
Solitary-wave solutions of (1) are sought in the form where µ is the propagation constant, and the amplitude function ψ(x) is real-valued and localized in space. Inserting (3) into (1), we find that ψ(x) satisfies For infinitesimal solutions ψ(x), the nonlinear term in the above equation drops out. The resulting linear version of (4) is a Mathieu-type equation and, by the Bloch-Floquet theory, its wave spectrum features a band-gap structure. Specifically, the linear version of (4) admits two linearly independent solutions in the form withp(x; µ) being π-periodic in x. The wave character of these so-called Bloch modes hinges on whether the corresponding wavenumber k is real or complex. By requiring k to be real, one then obtains an infinite number of Bloch bands for µ in which the Bloch modes (5) propagate. These propagation bands are separated by gaps, where k turns out to be complex, implying evanescent behaviour.
At the edges of Bloch bands, where the modes (5) switch from propagating to evanescent, two Bloch modes in the band collide at either k = 0 or k = ±1, and a single real-valued Bloch mode that is either π-or 2π-periodic, arises there. This suggests that edges of Bloch bands are possible bifurcation points of solitary-wave solutions of the nonlinear equation (4). These solitary waves reside inside band gaps and are the so-called gap solitary waves (or gap solitons).
The bifurcation of 1D gap solitons from band edges was discussed by Pelinovsky et al. (2004) using a multiple-scale expansion in powers of the wave amplitude, along with a certain constraint obeyed by locally confined solutions of (4). They identified two families of gap solitons that bifurcate from each band edge, namely 'on-site' and 'off-site' solitons, depending on whether the soliton is centred at a minimum or maximum of the sinusoidal potential (2).
More recently, Hwang et al. (2011) re-visited the problem of small-amplitude gap solitons near band edges, taking a different approach, that is applicable for a general potential V (x). Rather than the constraint utilized in Pelinovsky et al. (2004), Hwang et al. (2011) focused on the behaviour of the tails of Bloch-wave packets. These tails are controlled by the coupling of the wave envelope to the periodic Bloch mode at the band edge, an effect that lies beyond all orders of the usual multiple-scale expansion in powers of the wave amplitude. Hwang et al. (2011) carried this expansion beyond all orders using techniques of exponential asymptotics (Yang & Akylas 1997), for the case of Bloch-wave packets whose envelope features a single hump. It turns out that the tails of such wave packets decay as x → ±∞, and hence gap solitons arise, only for two specific locations of the wave envelope relative to the underlying periodic potential. Both these soliton families bifurcate from the linear (infinitesimal) periodic Bloch mode at a band edge, and they coincide with the on-site and off-site gap solitons found by Pelinovsky et al. (2004) in the case of a symmetric periodic potential.
In the present work, making use of the asymptotic expressions derived in Hwang et al. (2011) for the tails of a single-hump Block-wave packet, we shall construct small-amplitude gap-soliton families, in the form of bound states, that comprise two or more such wave packets. In contrast to the so-called fundamental gap solitons found in Pelinovsky et al. (2004) and Hwang et al. (2011), these higher-order soliton families bifurcate at small but finite amplitude close to a band edge, and they feature multiple branches that do not connect to fundamental soliton families or band edges.
In preparation for the ensuing analysis, we now summarize the main results of the multiplescale perturbation procedure for a single-hump Bloch-wave packet (Pelinovsky et al. 2004;Hwang et al. 2011). Close to a band edge µ = µ 0 , say, gap solitons are expected to be weakly nonlinear Bloch-wave packets. The solution to equation (4) is expanded in powers of an amplitude parameter, where η = ±1, p(x) ≡ p(x; µ 0 ), and X = ǫx is the 'slow' variable of the envelope function A(X). By imposing the appropriate solvability condition at O(ǫ 3 ), it turns out that A(X) satisfies the stationary NLS equation where The well-known soliton solution of (9) is where a = 2/α , β = |D| , and X 0 = ǫx 0 denotes the position of the peak of the envelope. This solution exists provided Dη < 0 and Dσ > 0; the first of these conditions requires µ to lie in the interior of the band gap, while the second condition can be met in the presence of self-focusing (σ = 1) nonlinearity if D > 0, or self-defocusing (σ = −1) nonlinearity if D < 0. It is important to note that the envelope equation (9) is translation-invariant, and X 0 is a free parameter in the solution (11). As a result, gap solitons of equation (4) could be obtained regardless of the position of the envelope (11) relative to the underlying periodic potential. This conclusion would seem rather suspicious, though, given that the original amplitude equation (4) is This issue was recognized by Pelinovsky et al. (2004), who pointed out that locally confined solutions of (4) must obey the constraint which can be readily obtained by multiplying (4) with ψ x and integrating with respect to x. Upon inserting the perturbation solution (6) and the potential (2) into (13), this condition then restricts the peak of the envelope to be at either a minimum (x 0 = 0) or a maximum (x 0 = π/2) of the potential (2), corresponding to on-site or off-site gap solitons, respectively. It is also worth noting that the so-called Melnikov function M (x 0 ) in (13), which depends on the shift x 0 of the envelope relative to the potential V (x), is exponentially small with respect to ǫ; hence, the constraint (13) for possible locations of gap solitons hinges upon information beyond all orders of the two-scale expansion (6).

Nonlocal solitary waves
Another way of reconciling the two-scale expansion (6) with the fact that gap solitons can be placed only at specific locations relative to the potential V (x), is by examining the behaviour of the tails of Bloch-wave packets near a band edge. Assuming that they have infinitesimal amplitude, these tails are governed by the linear version of equation (4). Hence, for µ inside a band gap and close to the edge µ 0 , using the same notation as in (7) and (10), the asymptotic representation of ψ(x) away from the solitary-wave core is, generically, a linear combination of two evanescent Bloch modes to leading order in |µ − µ 0 | = ǫ 2 ≪ 1. As expected, the slow exponential growth/decay of these modes is consistent with the behaviour of the tails of the envelope function A(X) according to the steady NLS equation (9). Now, for the purpose of constructing a solitary-wave solution of (4), the left-hand tail must For generic values of C + , when one integrates equation (4) from this left-hand tail to the right-hand side, the right-hand tail would involve both ψ − (x) and ψ + (x), so ψ(x) is not locally confined. However, the amplitude E + of the growing tail in (15) is exponentially small with respect to ǫ (see below), and hence this tail cannot be captured by an expansion in powers of ǫ, such as (6). Only when C + takes certain special values can the growing-tail amplitude E + vanish, resulting in a true solitary wave solution. This restriction on C + is consistent with the integral constraint (13) (6) beyond all orders of ǫ for a general periodic potential V (x), utilizing an exponential-asymptotics procedure in the wavenumber domain (Akylas & Yang 1995;Yang & Akylas 1997). According to this revised perturbation analysis, E + , which indeed is exponentially small with respect to ǫ, does vanish for two specific positions of the envelope (11) relative to the underlying potential, thus furnishing two families of gap solitons that bifurcate at the band edge. For the symmetric periodic potential (2), in particular, E + vanishes when the peak of the envelope is placed at x 0 = 0 or x 0 = π/2, corresponding to the on-site and off-site gap solitons, as was obtained earlier by Pelinovsky et al. (2004).
As the details of the exponential-asymptotics procedure are rather involved, we shall only quote the main results. When the peak of the envelope (11) is not at a minimum or maximum of the potential (x 0 = 0, π/2), the resulting 'solitary' waves are nonlocal owing to the presence of both growing and decaying solutions (14) in at least one of the tails. Specifically, assuming that the left-hand tail is locally confined in accordance with (6) and (11), the right-hand tail features a growing component of exponentially-small amplitude with respect to ǫ, in addition to a decaying component analogous to (16a): Note that the amplitude of the growing tail in (16b) is proportional to a constant C. As explained in Hwang et al. (2011), this constant is determined from solving a certain recurrence relation that contains information beyond all orders of expansion (6). In particular, C > 0 for the sinusoidal potential (2). As expected, the growing component of the right-hand tail (16b) is absent when x 0 = 0, π/2. By utilizing the symmetries of equation (4) as noted below, it is straightforward to deduce from (16) the asymptotic behaviour of nonlocal solitary waves whose right-hand tail is locally confined, but the left-hand tail comprises a growing and a decaying component: For the symmetric potential (2), equation (4) is invariant with respect to reflection in x (x → −x), and the Bloch wave p(x) at a band edge is either symmetric or antisymmetric in x, p(−x) = ±p(x).
Since (4) is also invariant with respect to ψ → −ψ, the reflected solution ψ may thus be written as (17) for both symmetric and antisymmetric p(x). The asymptotic expressions (16) and (17) for the tails of nonlocal solitary waves are key to constructing new families of locally confined solutions, in the form of bound states, as discussed below.

Bound states
As indicated above, gap solitons in the form of Bloch-wave packets with the 'sech' envelope (11) arise only when the peak of the envelope is at a minimum (x 0 = 0) or a maximum (x 0 = π/2) of the potential (2). Apart from these fundamental soliton states, it is possible, however, to construct other gap-soliton families by piecing together two or more of the nonlocal solitary waves found in §3. The overall strategy for determining these so-called bound states is similar to that followed in Yang & Akylas (1997) for constructing multi-packet solitary-wave solutions of the fifth-order KdV equation. Here, we shall work out the details of finding bound states involving two Bloch-wave packets.

Matching of tails
Consider two nonlocal solitary-wave solutions, ψ + (x) and ψ − (x), whose 'sech'-profile envelope functions are centred at x + 0 and −x − 0 , respectively, with x ± 0 > 0. In addition, let ψ ± (x) → 0 as x → ±∞, so the right-hand tail of ψ − (x) and the left-hand tail of ψ + (x) are nonlocal. Let us define the separation between the solitary-wave cores of the two constituent nonlocal waves as Then, assuming S ≫ 1/ǫ, (16b) and (17b) are legitimate asymptotic expressions for these tails in Specifically, according to (16b), the right-hand tail of ψ − (x) is given by and, according to (17b), the left-hand tail of ψ + (x) reads Note that the upper sign in (19) corresponds to the case where the envelopes of the two solitary waves have the same polarity (sign), while the lower sign pertains to the case of opposite polarity. Now, in order for these two nonlocal solitary waves to form a bound state, their growing and decaying tail components must match smoothly in the overlap region between the two main cores. Based on (19) and (20), this requires and, for this matching condition to be met, it is necessary that a 2πC Clearly, owing to the growing exponential e πβ/ǫ on the left-hand side of (22), this constraint cannot be satisfied for ǫ → 0 when S is finite. As a result, solution families of bound states are expected to bifurcate at a finite amplitude ǫ = ǫ c , say, depending on the separation distance S. Below, we shall verify this claim and compute ǫ c for the various solution branches.

Solution branches
Attention is now focused on the matching condition (21), in order to delineate the solution branches of bound states that bifurcate close to a band edge. Consider first the case where the two nonlocal solitary waves have envelopes with the same polarity. Then, the upper sign applies in (21), and sin 2x − 0 = sin 2x + 0 < 0. This condition can be satisfied in two ways: where m, n are integers and 0 < d < π/2. Accordingly, the solitary-wave separation (18) corresponding to these two possibilities is where N = m + n. Moreover, the matching condition (21) reads sin 2d = a 2πC with S given by (i) and (ii) in (24). It is convenient to label families of bound states by the integer N which, in view of (24), controls the separation between the two nonlocal solitary waves that form the bound states. As noted earlier, the present theory assumes that these solitary waves are well-separated, namely, S ≫ 1/ǫ, which in turn requires N ≫ 1/ǫ; the validity of this condition will be verified later in this section.
In preparation for tracing the solution branches of the family of bound states corresponding to a given N , we write d = π/4 + δ, with −π/4 < δ < π/4, so that from (24), we have where S 0 = (N + 1 2 )π. The matching condition (25), with S given by (i) and (ii) above, then leads to the following two conditions, respectively, where Based on equations (28), given ǫ, one may find the values of δ, and from (27) the corresponding separations S, for which bound states are possible. Specifically, from equation (28a), it is easy to see that two values of δ arise for each ǫ above a certain threshold, ǫ (1) c , which is the bifurcation point of the bound states obeying (28a). This threshold is associated with a double root, δ = δ From this equation and (28a), it follows that Inserting this result into (28a), ǫ c then is found by solving −π/4 π/4 c , it is clear from (28a) that cos 2δ → 0, and hence the two branches approach δ = ±π/4, so d = π/4 + δ → π/2, 0. Therefore, in this limit, in view of (23a), the two solution branches that bifurcate at ǫ (1) c approach limiting configurations of bound states, with both solitary waves being on-site on one branch (as δ → π/4) and off-site on the other branch (as δ → −π/4). In other words, as the solution curve in figure 1 is traced from one branch to the other, the two solitary waves making up the bound state transform simultaneously from on-site to off-site. In this transition, δ decreases from π/4 to −π/4, and according to (i) in (27), the separation distance S of the two solitary waves decreases from (N + 1)π to N π. This transition behaviour will be verified numerically in §6; see figure 2.
Turning next to the second possibility of bound states in (27), equation (28b) furnishes two values of δ, for each ǫ, as long as W (ǫ, N ) ≤ 1. Hence, ǫ must exceed a certain threshold, ǫ ≥ ǫ (2) c , where W (ǫ (2) c , N ) = 1. In view of (29), this bifurcation point is determined from the equation Here, the two values of δ in (33) actually correspond to equivalent bound-state configurations (under mirror reflection of the x-axis), as can be easily verified from (23b) with d = 1 4 π + δ. Therefore, (28b) describes only one solution branch. For ǫ ≫ ǫ (2) c , in particular, δ → ±π/4, and this solution branch approaches a limiting bound-state configuration that, according to (23b), involves an on-site and an off-site solitary wave.
We remark that the double root δ = δ (2) c = 0 of (28b) for ǫ = ǫ (2) c is also a solution of (28a) for the same value of ǫ. Hence, the solution branch (28b) bifurcates from the point δ = 0 on the lower branch of the solutions (28a), which themselves bifurcate at ǫ (1) c , as illustrated schematically in figure 1. As a consequence, the bound-state family corresponding to a certain N comprises three distinct solution branches which are connected with each other. The bifurcation points, ǫ (1) c and ǫ (2) c , do not coincide; but when N is large, they are close to each other. Indeed, in the limit N ≫ 1, it can be readily shown, using (28), (29) and (34), that ǫ c ∼ βN −1/2 . Thus, for large solitary-wave separation (N ≫ 1), bound-state families bifurcate near a band edge, but never exactly at the band edge itself. Moreover, this scaling of ǫ (1) c and ǫ (2) c in terms of N is consistent with the condition N ≫ 1/ǫ imposed earlier for the purpose of matching the tails (19) and (20) of nonlocal solitary waves forming a bound state.
We now consider the case when the two nonlocal solitary waves participating in a bound state have envelopes with opposite polarity, and the lower sign in (21) applies. Hence, sin 2x + 0 = sin 2x − 0 > 0. This condition can be satisfied in two ways: for some integers m, n and 0 < d < π/2. As a result, the expressions (24) for the solitary-wave separation S = x + 0 + x − 0 are also valid here, so the matching condition (21) gives rise to the same equations (28) that describe the solution branches δ = δ(ǫ) of bound states belonging to the family labelled by N = m + n, where d = π/4 + δ as before. Specifically, for a given N , the three solution branches found earlier arise in this instance as well: two of these branches bifurcate at ǫ = ǫ (1) c and, in the limit ǫ ≫ ǫ (1) c , represent bound states with both solitary waves being on-site (δ → −π/4) or off-site (δ → π/4). Marching along this solution curve from the on-site branch to the off-site branch, δ increases from −π/4 to π/4, and, in view of (35a), the separation distance S of the two nonlocal solitary waves increases from N π to (N + 1)π. (This transition behaviour is also brought out by the numerical results presented in §6; see figure 3.) This transition contrasts with the case of same envelope polarity examined earlier, where the separation distance S of the two nonlocal solitary waves decreases from (N + 1)π to N π when marching along the solution curve from the on-site branch to the off-site branch.

Power curves
In applications, gap-soliton solution branches are often described through which is commonly referred to as the soliton power. This quantity is also key to understanding the stability properties of gap solitons (Vakhitov & Kolokolov 1973;Shi et al. 2008;Yang 2010). Here, we shall derive an analytical expression for the power of small-amplitude bound states near band edges, and trace the solution branches found earlier in terms of P . The theoretical predictions will be compared against numerical results in §6. By virtue of perfect matching of the tails (19) and (20) in the overlap region −x − 0 ≪ x ≪ x + 0 , one may write the following uniformly valid approximation to bound states involving two nonlocal solitary waves (to the leading order): here again the upper (lower) sign corresponds to the case where the envelopes of the solitary waves have the same (opposite) polarity. Inserting (37) into (36), the power associated with a bound state may then be approximated as where The integrals above can be evaluated by substituting the Fourier series for the even π-periodic function p 2 (x), into (39) and (40) and then integrating term-by-term. Specifically, in the limit ǫ ≪ 1, ǫS ≫ 1, we find From (25), however, it is clear that ǫ 2 Se −ǫS/β ≫ e −πβ/ǫ ; hence, the second term in (42) is subdominant in comparison to (43), and (38) yields the following asymptotic expression for P , Consider first the case of bound states comprising two solitary waves with envelopes of the same polarity, where the upper sign in (44) applies. Corresponding to the first possibility for S in (27), there are two bound-state solution branches, δ = δ(ǫ), governed by (28a), and the associated power according to (44) is where S 0 = (N + 1 2 )π as before. On the other hand, corresponding to the second possibility for S in (27), there is only one solution branch, given by (33), and its power is In both (46) and (47), the second term in the brackets derives from the interaction of the tails of the nonlocal solitary waves. Moreover, it should be kept in mind that these expressions are valid when ǫ ≪ 1 and ǫS ≫ 1; i.e., close to the band edge and for wide separation of the two solitary waves forming a bound state. Based on the asymptotic formulae (46) and (47), it is possible to deduce the qualitative behaviour of the power curves. We recall that, for bound states of solitary waves with the same envelope polarity, δ decreases from π/4 to −π/4 as the solution curve (28a) is traversed from the on-site to the off-site branch. Since, for given ǫ, the power (46) is a decreasing function of δ (for ǫS > β), the power on the off-site branch is then always higher than the power on the on-site branch. Similarly, one can see that the power (47) of the mixed-site solution branch (28b) always lies between the powers of the on-site and off-site braches. Hence, the power curve associated with a bound-state solution family involving solitary waves of the same polarity, is such that the on-site branch is always the lower-power branch, the mixed-site branch lies in the middle and the off-site branch is the higher-power branch. Also, in the (ǫ, δ) plane, the middle solution branch bifurcates from the lower (off-site) branch (see figure 1); thus, on the power curve, the intermediate-power (mixed-site) branch bifurcates from the higher-power (off-site) branch. These general features of the power curves will be verified numerically in §6 (see figures 2 and 5).
Next, in the case of bound states with nonlocal solitary waves having envelopes of opposite polarity, where the lower sign in (44) applies, a similar calculation based on (35), along with (27), yields for the two solution branches bifurcating at ǫ (1) c , and for the single solution branch bifurcating at ǫ (2) c . Based on these power formulae, it is again possible to show that, on the associated power curve, the on-site branch is the lower-power branch, the mixed-site branch is the immediate-power branch, and the off-site branch is the higher-power branch, just as in the case of same envelope polarity above. However, the intermediate-power branch now bifurcates from the lower-power branch, which is the opposite of the conclusion reached earlier in the same-envelope-polarity case. These theoretical predictions are confirmed by numerical results below (see figure 3).

Numerical results
We now turn to a numerical investigation of bound states in order to make a comparison of numerical results against the predictions of the asymptotic theory discussed earlier. To this end, equation (4) is solved numerically by the Newton-conjugate-gradient method (Yang 2010), using the sinusoidal potential (2) with V 0 = 6. For this value of the potential depth, when σ = 1 (self-focusing nonlinearity), gap solitons bifurcate from the lower edge µ 0 = 2.061318 of the first Bloch band, where the diffraction coefficient D is positive, and the parameters (12) with C = 3.0133. We begin with a few qualitative comparisons between the analysis and numerics. For this purpose, we assume self-focusing nonlinearity (σ = 1) and consider bound states involving solitary waves with envelopes of the same polarity for N = 10. The numerical results for this family of bound states are displayed in figure 2. From the power curve in figure 2a, it is seen that this family indeed exhibits three distinct solution branches that bifurcate near the band edge, as predicted by the theory. On the lower-power branch, the bound state at point c comprises two on-site gap solitons which are separated by (N + 1)π = 11π (see figure 2c); on the higher-power branch, the bound state at point e comprises two off-site gap solitons which are separated by N π = 10π (see figure 2e). Thus, as one moves from the lower branch to the upper branch along the power curve, the two solitary waves forming a bound state transform from on-site to off-site, and their separation decreases by one lattice period (π). On the middle branch, however, the bound state at point d comprises an on-site and an off-site gap soliton which are separated by (N + 1/2)π = 10.5π (see figure 2d), confirming that this is the mixed-site branch. At the bifurcation point near the band edge (tip of the power curve), the bound-state profile is displayed in figure 2f. This bound state comprises two low-amplitude Bloch-wave packets which are well separated, as assumed in our asymptotic analysis. These features of the power curve and the corresponding bound-state profiles are in perfect qualitative agreement with the asymptotic theory. Furthermore, from the amplification of the power curve near the bifurcation point, shown in figure 2b, it is clear that the middle branch bifurcates from the higher-power branch, again in agreement with the previous analysis in the case of bound states with solitary waves of the same envelope polarity.
Next, we turn to quantitative comparison between the analysis and numerics. For this purpose, we choose to consider bound states comprising solitary waves with opposite envelope polarity (the results of comparison for bound states with same envelope polarity are very similar). Specifically, once again, we assume self-focusing nonlinearity (σ = 1) and take N = 10, but the envelopes of the two solitary waves in the bound state now have opposite polarity. Figure 3a  power curves are also in reasonable agreement, especially considering that the analytical bifurcation point here is ǫ (1) c = 0.1368, which is not all that small; moreover, for this value of ǫ and N = 10, ǫS ∼ 4, which is not all that large. Figure 3c,d shows quantitative comparison of numerics and analysis for the bound-state profiles on the lower-and higher-power branches at µ = 2.0438 (near the bifurcation point). The corresponding analytical bound-state profiles (red dashed curves) are given by (37). Here, the locations of the two solitary waves involved in the bound state, namely, x ± 0 , are determined by (35); these values are dependent on δ(ǫ), which is found by solving equations (28) with ǫ = |µ 0 − µ|. The analytical solution profiles show good agreement with the numerical ones.
As another quantitative comparison, we examine the dependence of the bifurcation point ǫ  solitons with the same as well as opposite envelope polarities and self-focusing nonlinearity (σ = 1). Here, the analytical bifurcation point ǫ (1) c is obtained from solving equation (32), and this value is the same for both cases of envelope polarity (with the same N ). For both envelope polarities, the numerical values for ǫ (1) c approach the analytical ones when N → ∞, confirming the asymptotic accuracy of our analysis.
As a final comparison, we assume self-defocusing nonlinearity (σ = −1) and consider bound states involving solitary waves with the same envelope polarity. For N = 10, the numerical results for this family of bound states are displayed in figure 5. The power curve features three branches (figure 5a), and the middle branch bifurcates from the upper branch (figure 5b), in agreement with the asymptotic analysis. In addition, the bound state at point c on the lower branch comprises two on-site solitons which are separated by (N + 1)π = 11π (see figure 5c), and the bound state at point d on the upper branch comprises two off-site solitons which are separated by N π = 10π, again consistent with the previous analysis. It may seem puzzling at first sight that the two on-site solitons in figure 5c are of opposite signs, given that the envelopes of these solitons ought to have the same polarity. To explain this feature, note that at the band edge relevant here, µ 0 = 2.266735, which is the upper edge of the first Bloch band, the Bloch mode p(x) is 2π-periodic, and its adjacent peaks have opposite signs, similar to the function cos x (see Yang 2010, §6.1.1). Since the two onsite gap solitons in figure 5c are separated by (N + 1)π = 11π, the Bloch function p(x) at the centres of these solitons then has opposite sign. Thus, even though the envelopes of the two gap solitons have the same polarity, the overall gap-soliton profiles, which are products of the Bloch function and the envelope function, have opposite sign, as in figure 5c.

Concluding remarks
In this paper, an asymptotic theory was developed for stationary gap-soliton bound states consisting of two fundamental gap solitons in 1D periodic media. It was shown that there is a countable set of such bound state families, characterized by the separation distance of the two solitary waves involved, and each family features three distinct solution branches that bifurcate near band edges at small, but finite, amplitude. Of these three solution branches, one branch contains bound states whose fundamental solitons are both on-site; the second branch contains bound states whose fundamental solitons are both off-site; and the third branch represents bound states in which one fundamental soliton is on-site and the other off-site. The power curves associated with these solution branches were computed asymptotically for large solitary-wave separation. On the power diagram, the on-site solution branch is always the lower-power branch, the mixed-site branch the middle branch, and the off-site branch is the higher-power branch; in addition, the middle branch bifurcates from the upper (lower) branch when the envelopes of the two fundamental solitons in the bound state have the same (opposite) polarity. These analytical results were compared with numerical results and good agreement was obtained.
There are several common features between the results of this article and those obtained by Yang & Akylas (1997) for two-wavepacket bound states in the fifth-order KdV equation. In that case, a countable set of bound states also bifurcate at finite values of the wave amplitude, and each solution family also features multiple branches, with asymmetric-wave branches bifurcating off symmetric-wave branches. Furthermore, the techniques for constructing bound states in these two models follow along the same lines. Even though the current lattice model (4) is not translation invariant, while the fifth-order KdV equation is, remarkably, both models admit very similar classes of solitary-wave solutions.
In addition to being of fundamental interest, the bound states discussed here could also prove useful in applications as they allow increased flexibility in the profiles of localized nonlinear modes within band gaps. In a recent demonstration of image transmission in photonic lattices, in fact, Yang et al. (2011) utilized 2D higher-order gap solitons of various shapes, and those gap solitons are closely related to the 1D bound states constructed in the present paper. For the purpose of assessing the potential usefulness of these 1D bound states, of course, it is necessary to examine the stability properties of the various bound-state families, a question that will be taken up in future studies.
In this paper, attention was focused on bound states consisting of only two fundamental gap solitons; extending the analysis to more complicated bound states involving three or more fundamental solitons, as well as to bound states in 2D periodic media, will be left to future studies.
Our final remark is that, since the bound states constructed in this paper bifurcate near Blochband edges at finite amplitude, their power curves always terminate before band edges and never reach them (see figures 2a and 5a). This feature of the power curve is shared by some other types of gap solitons, such as the truncated-Bloch-wave solitons (Alexander et al. 2006;Yang 2010). Whether the present analysis can be applied to those gap solitons needs further investigation.