## Abstract

Solitary waves in one-dimensional periodic media are discussed by employing the nonlinear Schrödinger 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 non-local 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 non-local solitary waves, high-order locally confined gap solitons, or bound states, are constructed. Details are worked out for bound states comprising two non-local 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.

## 1. 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; Kivshar & Agrawal 2003; Morsch & Oberthaler 2006; Skorobogatiy & 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 (Yang *et al.* 2011).

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 band gaps in which linear disturbances are evanescent but nonlinear-localized modes, commonly known as gap solitons, are possible. Following the first theoretical prediction (Christodoulides & Joseph 1988) and experimental observation (Eisenberg *et al.* 1998) of fundamental solitons in the semi-infinite band gap of one-dimensional periodic waveguides, various types of gap solitons in one- and multi-dimensions have been reported theoretically and experimentally. Examples include two-dimensional fundamental gap solitons, vortex solitons, dipole solitons, reduced-symmetry solitons, vortex-array 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, 2005; Shi *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. Some progress has been made, in this issue, particularly for gap solitons that bifurcate from linear Bloch modes at the edges of Bloch bands. Specifically, in one-dimensional periodic media, only two gap-soliton families bifurcate from each edge of a Bloch band under self-focusing or self-defocusing 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 two-dimensional periodic media, 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 one-dimensional 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 non-local owing to the presence of growing tails of exponentially small amplitude with respect to the peak wave amplitude. However, by piecing together such non-local solitary waves, it is possible to construct high-order gap solitons, or bound states. The proposed asymptotic theory is illustrated by working out details for bound states involving two non-local 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.

## 2. Preliminaries

Our study of nonlinear wave phenomena in periodic media is based on the one-dimensional NLS equation,
2.1
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
2.2
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 equation (2.1) are sought in the form
2.3
where *μ* is the propagation constant, and the amplitude function *ψ*(*x*) is real-valued and localized in space. Inserting equation (2.3) into equation (2.1), we find that *ψ*(*x*) satisfies
2.4
For infinitesimal solutions *ψ*(*x*), the nonlinear term in the above equation drops out. The resulting linear version of equation (2.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 equation (2.4) admits two linearly independent solutions in the form
2.5
with 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 (2.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 (2.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 (2.4). These solitary waves reside inside band gaps and are the so-called gap solitary waves (or gap solitons).

The bifurcation of one-dimensional 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 equation (2.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.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 in Pelinovsky *et al.* (2004), it 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 , 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, which 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 high-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 multiple-scale 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 (2.4) is expanded in powers of an amplitude parameter, 0<*ϵ*≪1,
2.6
along with
2.7
where *η*=±1,
2.8
*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
2.9
where
2.10
The well-known soliton solution of equation (2.9) is
2.11
where
2.12
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, whereas 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 (2.9) is translation-invariant, and *X*_{0} is a free parameter in the solution (2.11). As a result, gap solitons of equation (2.4) could be obtained regardless of the position of the envelope (2.11) relative to the underlying periodic potential. This conclusion would seem rather suspicious, though, given that the original amplitude equation (2.4) is *not* translation-invariant.

This issue was recognized by Pelinovsky *et al.* (2004), who pointed out that locally confined solutions of equation (2.4) must obey the constraint
2.13
which can be readily obtained by multiplying equation (2.4) with *ψ*_{x} and integrating with respect to *x*. Upon inserting the perturbation solution (2.6) and the potential (2.2) into equation (2.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.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 equation (2.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 (2.13) for possible locations of gap solitons hinges upon information beyond all orders of the two-scale expansion (2.6).

## 3. Non-local solitary waves

Another way of reconciling the two-scale expansion (2.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 (2.4). Hence, for *μ* inside a band gap and close to the edge *μ*_{0}, using the same notation as in equations (2.7) and (2.10), the asymptotic representation of *ψ*(*x*) away from the solitary-wave core is, generically, a linear combination of two evanescent Bloch modes
3.1
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 (2.9).

Now, for the purpose of constructing a solitary-wave solution of equation (2.4), the left-hand tail must involve only *ψ*_{+}(*x*), *ψ*(*x*)∼*ϵC*_{+}*ψ*_{+}(*x*), so that *ψ*→0 as . For generic values of *C*_{+}, when one integrates equation (2.4) from this left-hand tail to the right-hand side, the right-hand tail would involve both *ψ*_{−}(*x*) and *ψ*_{+}(*x*),
3.2
so *ψ*(*x*) is not locally confined. 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 (2.13) used by Pelinovsky *et al.* (2004). However, the amplitude *E*_{+} of the growing tail in equation (3.2) is exponentially small with respect to *ϵ* (see below), and hence this tail cannot be captured by an expansion in powers of *ϵ*, such as equation (2.6).

Hwang *et al.* (2011) computed the amplitude *E*_{+} by carrying expansion (2.6) beyond all orders of *ϵ* for a general periodic potential *V* (*x*), using 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 (2.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.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 quote only the main results. When the peak of the envelope (2.11) is not at a minimum or maximum of the potential (*x*_{0}≠0,*π*/2), the resulting ‘solitary’ waves are non-local owing to the presence of both growing and decaying solutions (3.1) in at least one of the tails. Specifically, assuming that the left-hand tail is locally confined in accordance with equations (2.6) and (2.11),
3.3a
the right-hand tail features a growing component of exponentially small amplitude with respect to *ϵ*, in addition to a decaying component analogous to equation (3.3a):
3.3b
Note that the amplitude of the growing tail in equation (3.3b) 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 (2.6). In particular, *C*>0 for the sinusoidal potential (2.2). As expected, the growing component of the right-hand tail (3.3b) is absent when *x*_{0}=0,*π*/2.

Using the symmetries of equation (2.4) as noted below, it is straightforward to deduce from equation (3.3) the asymptotic behaviour of non-local solitary waves whose right-hand tail is locally confined,
3.4a
but the left-hand tail comprises a growing and a decaying component:
3.4b
For the symmetric potential (2.2), equation (2.4) is invariant with respect to reflection in *x* , and the Bloch wave *p*(*x*) at a band edge is either symmetric or antisymmetric in *x*, *p*(−*x*)=±*p*(*x*). As equation (2.4) is also invariant with respect to *ψ*→−*ψ*, the reflected solution *ψ* may thus be written as equation (3.4) for both symmetric and antisymmetric *p*(*x*). The asymptotic expressions (3.3) and (3.4) for the tails of non-local solitary waves are key to constructing new families of locally confined solutions, in the form of bound states, as discussed below.

## 4. Bound states

As indicated above, gap solitons in the form of Bloch-wave packets with the ‘sech’ envelope (2.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.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 non-local solitary waves found in §3. The overall strategy for determining these so-called bound states is similar to that followed by Yang & Akylas (1997) for constructing multi-packet solitary-wave solutions of the fifth-order Korteweg–de Vries (KdV) equation. Here, we shall work out the details of finding bound states involving two Bloch-wave packets.

### (a) Matching of tails

Consider two non-local solitary-wave solutions, *ψ*^{+}(*x*) and *ψ*^{−}(*x*), whose ‘sech’-profile envelope functions are centred at and , respectively, with . In addition, let *ψ*^{±}(*x*)→0 as , so the right-hand tail of *ψ*^{−}(*x*) and the left-hand tail of *ψ*^{+}(*x*) are non-local. Let us define the separation between the solitary-wave cores of the two constituent non-local waves as
4.1
Then, assuming *S*≫1/*ϵ*, equations (3.3b) and (3.4b) are legitimate asymptotic expressions for these tails in the ‘overlap’ region . Specifically, according to equation (3.3b), the right-hand tail of *ψ*^{−}(*x*) is given by
4.2
and, according to equation (3.4b), the left-hand tail of *ψ*^{+}(*x*) reads
4.3
Note that the upper sign in equation (4.2) corresponds to the case where the envelopes of the two solitary waves have the same polarity (sign), whereas the lower sign pertains to the case of opposite polarity.

Now, in order for these two non-local 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 equations (4.2) and (4.3), this requires
4.4
and, for this matching condition to be met, it is necessary that
4.5
Clearly, owing to the growing exponential e^{πβ/ϵ} on the left-hand side of equation (4.5), 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.

### (b) Solution branches

Attention is now focused on the matching condition (4.4), in order to delineate the solution branches of bound states that bifurcate close to a band edge.

Consider first the case where the two non-local solitary waves forming a bound state have envelopes with the same polarity. Then, the upper sign applies in equation (4.4), and . This condition can be satisfied in two ways:
4.6a
and
4.6b
where *m*,*n* are integers and 0<*d*<*π*/2. Accordingly, the solitary-wave separation (4.1) corresponding to these two possibilities is
4.7
where *N*=*m*+*n*. Moreover, the matching condition (4.4) reads
4.8
with *S* given by (i) and (ii) in equation (4.7).

It is convenient to label families of bound states by the integer *N* which, in view of equation (4.7), controls the separation between the two non-local 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 §4*b*.

In preparing to trace the solution branches of the family of bound states corresponding to a given *N*, we write
4.9
with −*π*/4<*δ*<*π*/4, so that from equation (4.7), we have
4.10
where *S*_{0}=(*N*+1/2)*π*. The matching condition (4.8), with *S* given by (i) and (ii) above, then leads to the following two conditions, respectively,
4.11a
and
4.11b
where
4.12
Based on equations (4.11), given *ϵ*, one may find the values of *δ*, and from equation (4.10) the corresponding separations *S*, for which bound states are possible.

Specifically, from equation (4.11a), it is easy to see that two values of *δ* arise for each *ϵ* above a certain threshold, , which is the bifurcation point of the bound states obeying equation (4.11a). This threshold is associated with a double root, , of equation (4.11a), where
4.13
From equations (4.13) and (4.11a), it follows that
4.14
Inserting this result into equation (4.11a), then is found by solving
4.15
For *N*≫1, the two solution branches of *δ*(*ϵ*), obtained from equation (4.11a) for , are plotted schematically in figure 1. The bifurcation point is the turning point of this double-branch curve, and . For , it is clear from equation (4.11a) that , and hence the two branches approach *δ*=±*π*/4, so *d*=*π*/4+*δ*→*π*/2,0. Therefore, in this limit, in view of equation (4.6a), the two solution branches that bifurcate at 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 equation (4.10), the separation distance *S* of the two solitary waves decreases from (*N*+1)*π* to *Nπ*. This transition behaviour will be verified numerically in §6 (figure 2).

Turning next to the second possibility of bound states in equation (4.10), equation (4.11b) furnishes two values of *δ*,
4.16
for each *ϵ*, as long as *W*(*ϵ*,*N*)≤1. Hence, *ϵ* must exceed a certain threshold, , where . In view of equation (4.12), this bifurcation point is determined from the equation
4.17
Here, the two values of *δ* in equation (4.16) actually correspond to equivalent bound-state configurations (under mirror reflection of the *x*-axis), as can be easily verified from equation (4.6b) with *d*=(1/4)*π*+*δ*. Therefore, equation (4.11b) describes only one solution branch. For , in particular, *δ*→±*π*/4, and this solution branch approaches a limiting bound-state configuration that, according to equation (4.6b), involves an on-site and an off-site solitary wave.

We remark that the double root of equation (4.11b) for is also a solution of equation (4.11a) for the same value of *ϵ*. Hence, the solution branch (4.11b) bifurcates from the point *δ*=0 on the lower branch of the solutions (4.11a), which themselves bifurcate at , 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, and , 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 equations (4.11), (4.12) and (4.17), that . 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 and in terms of *N* is consistent with the condition *N*≫1/*ϵ* imposed earlier for the purpose of matching the tails (4.2) and (4.3) of non-local solitary waves forming a bound state.

We now consider the case when the two non-local solitary waves participating in a bound state have envelopes with opposite polarity, and the lower sign in equation (4.4) applies. Hence, . This condition can be satisfied in two ways:
4.18a
and
4.18b
for some integers *m*,*n* and 0<*d*<*π*/2. As a result, the expressions (4.7) for the solitary-wave separation are also valid here, so the matching condition (4.4) gives rise to the same equations (4.11) 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 and, in the limit , 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 equation (4.18a), the separation distance *S* of the two non-local solitary waves increases from *Nπ* to (*N*+1)*π*. (This transition behaviour is also brought out by the numerical results presented in §6; figure 3.) This transition contrasts with the case of same envelope polarity examined earlier, where the separation distance *S* of the two non-local solitary waves decreases from (*N*+1)*π* to *Nπ* when marching along the solution curve from the on-site branch to the off-site branch.

Finally, according to equations (4.11b), (4.16) and (4.17), the third solution branch bifurcates at and, for , represents bound states with one solitary wave on-site and the other off-site. Moreover, it follows from equation (4.18) that, in the course of continuation of the bound-state solution from the on-site branch to this mixed-site branch, one of the two on-site solitary waves remains on-site and does not move, while the other on-site solitary wave moves away and becomes off-site.

## 5. Power curves

In applications, gap-soliton solution branches are often described through
5.1
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 (4.2) and (4.3) in the overlap region , one may write the following uniformly valid approximation to bound states involving two non-local solitary waves (to the leading order): 5.2 here again the upper (lower) sign corresponds to the case where the envelopes of the solitary waves have the same (opposite) polarity. Inserting equation (5.2) into equation (5.1), the power associated with a bound state may then be approximated as 5.3 where 5.4 and 5.5

The integrals above can be evaluated by substituting the Fourier series
5.6
for the even *π*-periodic function *p*^{2}(*x*), into equations (5.4) and (5.5) and then integrating term-by-term. Specifically, in the limit *ϵ*≪1, *ϵS*≫1, we find
5.7
and
5.8
From equation (4.8), however, it is clear that *ϵ*^{2}*S*e^{−ϵS/β}≫e^{−πβ/ϵ}; hence, the second term in equation (5.7) is sub-dominant in comparison with equation (5.8), and equation (5.3) yields the following asymptotic expression for *P*,
5.9
with
5.10

Consider first the case of bound states comprising two solitary waves with envelopes of the same polarity, where the upper sign in equation (5.9) applies. Corresponding to the first possibility for *S* in equation (4.10), there are two bound-state solution branches, *δ*=*δ*(*ϵ*), governed by equation (4.11a), and the associated power according to equation (5.9) is
5.11
where *S*_{0}=(*N*+1/2)*π* as before. On the other hand, corresponding to the second possibility for *S* in equation (4.10), there is only one solution branch, given by equation (4.16), and its power is
5.12
In both equations (5.11) and (5.12), the second term in the brackets derives from the interaction of the tails of the non-local 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 (5.11) and (5.12), 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 (4.11a) is traversed from the on-site to the off-site branch. As, for given *ϵ*, the power (5.11) 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 (5.12) of the mixed-site solution branch (4.11b) always lies between the powers of the on-site and off-site branches. 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 low-power branch, the mixed-site branch lies in the middle and the off-site branch is the high-power branch. Also, in the (*ϵ*,*δ*) plane, the middle solution branch bifurcates from the lower (off-site) branch (figure 1); thus, on the power curve, the intermediate-power (mixed-site) branch bifurcates from the high-power (off-site) branch. These general features of the power curves will be verified numerically in §6 (figures 2 and 5).

Next, in the case of bound states with non-local solitary waves having envelopes of opposite polarity, where the lower sign in equation (5.9) applies, a similar calculation based on equation (4.18), along with equation (4.10), yields 5.13 for the two solution branches bifurcating at , and 5.14 for the single solution branch bifurcating at . Based on these power formulae, it is again possible to show that, on the associated power curve, the on-site branch is the low-power branch, the mixed-site branch is the immediate-power branch, and the off-site branch is the high-power branch, just as in the case of same envelope polarity above. However, the intermediate-power branch now bifurcates from the low-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 in §6 (figure 3).

## 6. 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 (2.4) is solved numerically by the Newton-conjugate-gradient method (Yang 2010), using the sinusoidal potential (2.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 (2.12) of the envelope soliton (2.11) take the values
6.1
Moreover, the constant *C* in the tail expressions (3.3b) and (3.4b) is found by solving a certain recurrence relation as discussed in Hwang *et al.* (2011): *C*=1.307. On the other hand, when *σ*=−1 (self-defocusing nonlinearity), gap solitons bifurcate from the upper edge *μ*_{0}=2.266735 of the first Bloch band, where the diffraction coefficient is negative, and the soliton parameter values are
6.2
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 2*a*, 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 low-power branch, the bound state at point c comprises two on-site gap solitons which are separated by (*N*+1)*π*=11*π* (figure 2*c*); on the high-power branch, the bound state at point e comprises two off-site gap solitons which are separated by *Nπ*=10*π* (figure 2*e*). 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*π* (figure 2*d*), 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 2*f*. 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 2*b*, it is clear that the middle branch bifurcates from the high-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 3*a* displays the numerical power diagram near the bifurcation point *μ*_{c}=2.04442, and figure 3*b* illustrates the analytical power curves given by equations (5.13) and (5.14) near the analytical bifurcation point , where as obtained from equation (4.15). (To compute the analytical power curves, we first solve equations (4.11) for *δ* as a function of *ϵ* for , and then use .) Both the numerical and analytical power curves in figure 3*a*,*b* exhibit three branches, as expected. In addition, the middle (mixed-site) branch now bifurcates from the low-power (on-site) branch, again in agreement with the analysis in the end of §5. Quantitatively, the numerical and analytical power curves are also in reasonable agreement, especially considering that the analytical bifurcation point here is , which is not all that small; moreover, for this value of *ϵ* and *N*=10, *ϵS*∼4, which is not all that large. Figure 3*c,d* shows quantitative comparison of numerics and analysis for the bound-state profiles on the low- and high-power branches at *μ*=2.0438 (near the bifurcation point). The corresponding analytical bound-state profiles (red dashed curves) are given by equation (5.2). Here, the locations of the two solitary waves involved in the bound state, namely, , are determined by equation (4.18); these values are dependent on *δ*(*ϵ*), which is found by solving equations (4.11) with . The analytical solution profiles show good agreement with the numerical ones.

As another quantitative comparison, we examine the dependence of the bifurcation point on the parameter *N* that controls the separation distance of the two solitary waves in a bound state. The analytical and numerical values of this bifurcation point against *N* are displayed in figure 4 for solitons with the same as well as opposite envelope polarities and self-focusing nonlinearity (*σ*=1). Here, the analytical bifurcation point is obtained from solving equation (4.15), and this value is the same for both cases of envelope polarity (with the same *N*). For both envelope polarities, the numerical values for approach the analytical ones when , 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 5*a*), and the middle branch bifurcates from the upper branch (figure 5*b*), 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*π* (figure 5*c*), 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 5*c* 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 (see Yang 2010, §6.1.1). As the two on-site gap solitons in figure 5*c* 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 signs, as in figure 5*c*.

## 7. Concluding remarks

In this paper, an asymptotic theory was developed for stationary gap-soliton bound states consisting of two fundamental gap solitons in one-dimensional 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 low-power branch, the mixed-site branch the middle branch and the off-site branch is the high-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 (2.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) used two-dimensional high-order gap solitons of various shapes, and those gap solitons are closely related to the one-dimensional bound states constructed in the present paper. For the purpose of assessing the potential usefulness of these one-dimensional 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 two-dimensional periodic media, will be left to future studies.

Our final remark is that, as the bound states constructed in this paper bifurcate near Bloch-band edges at finite amplitude, their power curves always terminate before band edges and never reach them (figures 2*a* and 5*a*). 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.

## Acknowledgements

The work of T.R.A. was supported in part by the National Science Foundation (grant DMS-098122), and the work of G.H. and J.Y. was supported in part by the Air Force Office of Scientific Research (grant USAF 9550-09-1-0228) and the National Science Foundation (grant DMS-0908167).

- Received May 25, 2011.
- Accepted August 1, 2011.

- This journal is © 2011 The Royal Society