Dynamical Barrier for the Formation of Solitary Waves in Discrete Lattices
P.G. Kevrekidis, J.A. Espinola-Rocha, Y. Drossinos, A. Stefanov
aa r X i v : . [ n li n . PS ] N ov Dynamical Barrier for the Formation of Solitary Waves in Discrete Lattices
P.G. Kevrekidis and J.A. Espinola-Rocha
Department of Mathematics and Statistics, University of Massachusetts, Amherst, MA 01003
Y. Drossinos
European Commission, Joint Research Centre,I-21020 Ispra (Va), ItalyandSchool of Mechanical and Systems Engineering,University of Newcastle upon Tyne,Newcastle upon Tyne NE1 7RU, United Kingdom
A. Stefanov
Department of Mathematics, University of Kansas, 1460 Jayhawk Blvd., Lawrence, KS 66045–7523 (Dated: October 31, 2018)We consider the problem of the existence of a dynamical barrier of “mass” that needs to beexcited on a lattice site to lead to the formation and subsequent persistence of localized modes for anonlinear Schr¨odinger lattice. We contrast the existence of a dynamical barrier with its absence inthe static theory of localized modes in one spatial dimension. We suggest an energetic criterion thatprovides a sufficient, but not necessary, condition on the amplitude of a single-site initial conditionrequired to form a solitary wave. We show that this effect is not one-dimensional by consideringits two-dimensional analog. The existence of a sufficient condition for the excitation of localizedmodes in the non-integrable, discrete, nonlinear Schr¨odinger equation is compared to the dynamicsof excitations in the integrable, both discrete and continuum, version of the nonlinear Schr¨odingerequation.
I. INTRODUCTION
In the past few years, there has been an explosion of interest in discrete models that has been summarized innumerous recent reviews [1]. This growth has been spurted by numerous applications of dynamical lattice nonlinearmodels in areas as diverse as the nonlinear optics of waveguide arrays [2], the dynamics of Bose-Einstein condensatesin periodic potentials [3], micro-mechanical models of cantilever arrays [4], or even simple models of the complexdynamics of the DNA double strand [5]. Perhaps the most prototypical model among those that emerge in thesesettings is the, so-called, discrete nonlinear Schr¨odinger equation (DNLS) [6]. DNLS may arise as a direct model, asa tight binding approximation, or even as an envelope wave expansion: it is, arguably, one of the most ubiquitousmodels in the nonlinear physics of dispersive, discrete systems.In at least one of these settings (namely, in the nonlinear optics of waveguide arrays with the focusing nonlinearity),the feature that will be of interest to the present work has been observed experimentally. In particular, it has beennoted, to the best of our knowledge firstly in Ref. [7], that when an injected beam of light into one waveguide hadlow intensity, then the beam dispersed through quasi-linear propagation. On the other hand, in the same work,experiments with high intensity of the input beam led to the first example of formation of discrete solitary waves inwaveguide arrays. A very similar “crossover” from linear to nonlinear behavior was also observed very recently in arraysof waveguides with the defocusing nonlinearity [8]. The common feature of both works is that they used the DNLSequation as the supporting model to illustrate this behavior at a theoretical/numerical level. However, this crossoverphenomenon is certainly not purely discrete in nature. Perhaps the most famous example of a nonlinear wave equationthat possesses such a threshold is the integrable continuum nonlinear Schr¨odinger equation [9]. Specifically, it is well-known that, e.g., in the case of a square barrier of initial conditions of amplitude V and width L , the product V L determines the nature of the resulting soliton, and if it is sufficiently small the initial condition disperses without theformation of a solitonic structure [10]. On the other hand, the existence of the threshold is not a purely one-dimensionalfeature either. For instance, experiments on the formation of solitary waves in two-dimensional photorefractive crystalsshow that low intensities lead to diffraction, whereas higher intensities induce localization [11, 12]. Moreover, similarphenomena were observed even in the formation of higher-order excited structures such as vortices (as can be inferredby carefully inspecting the results of Refs. [13, 14]). It should be mentioned that the latter field of light propagationin photorefractive crystals is another major direction of current research in nonlinear optics; see, e.g., Ref. [15] for arecent review.This crossover behavior between linear and nonlinear dynamics may be understood qualitatively rather simply. Inthe case of power law nonlinearities of order p , which are relevant in these settings, a small intensity O( ǫ ), where ǫ ≪
1, yields a nonlinear contribution O( ǫ p ) that is negligible with respect to the linear terms of the equation. Onthe other hand, if ǫ ≫
1, the opposite will be true and the nonlinear terms will dominate the linear ones, yieldingessentially nonlinear behavior. A key question regards the details of this crossover and what determines its moreprecise location for an appropriately parametrized initial condition. This is the question we address herein. We arguethat the problem related to the experiments described above reduces, at the mathematical level, to a DNLS equationwith a Kronecker- δ initial condition parametrized by its amplitude. Then, a well defined value of the initial-stateamplitude exists such that initial states with higher amplitude always give rise to localized modes. The conditionmay be determined by comparing the energy of the initial state with the energy of the localized excitations that themodel supports. This sufficient, but not necessary, condition for the formation of localized solitary waves provides anintuitively and physically appealing interpretation of the dynamics that is in very good agreement with our numericalobservations. We also consider variants of this process in different settings: for reasons of completeness, we present italso in the continuum NLS equation, noting the significant differences that the latter case has from the present one.As yet another example of very different (from both its non-integrable sibling and its continuum limit) dynamicalbehavior, we also present the case of the integrable discrete NLS (so-called Ablowitz-Ladik [16, 17, 18]) model. Inaddition to the one-dimensional DNLS lattice, we also consider the two-dimensional case where the role of bothenergy and beam power (mathematically the squared l -norm) become apparent. We should note here that our toolof choice for visualizing the “relaxational process” (albeit in a Hamiltonian system) of the initial condition will beenergy-power diagrams. Such diagrams have proven very helpful in visualizing the dynamics of initial conditions in adiverse host of nonlinear wave equations. In particular, they have been used in the nonlinear homogeneous systemssuch as birefringent media and nonlinear couplers as is discussed in Chapters 7 and 8 of Ref. [19]. They have alsobeen used in a form closely related to the present work (but in the continuum case; see also the discussion below)for general nonlinearities in dispersive wave equations in [20], while they have been used to examine the migration oflocalized excitations in DNLS equations in [21].Our presentation is structured as follows. In section II, we present the analytically tractable theory of the integrable“relatives” of the present model: we review the known theory for the NLS model and develop its analog for theintegrable discrete NLS case. Then, in section III, we present our analytical and numerical results in the one- andtwo-dimensional DNLS equation. In the last section, we summarize our findings and present our conclusions, as wellas highlight some important questions for future studies. II. THRESHOLD CONDITIONS FOR THE INTEGRABLE NLS MODELSA. The Continuum NLS Model
For reasons of completeness of the presentation and to compare and contrast the results of the non-integrable casethat is at the focus of the present work, we start by summarizing the threshold conditions for the continuum NLSmodel [10]. For the focusing NLS equation iu t = − u xx − | u | u, (1)with squared barrier initial data u ( x,
0) = (cid:26) V , − L ≤ x ≤ L , otherwise , (2)(the inverse of) the transmission coefficient, S ( E ), which is the first entry of the scattering matrix, is given by S ( E ) = ν ( E ) cos(2 ν ( E ) L ) − iE sin(2 ν ( E ) L ) , (3)with ν ( E ) = p E + V where E is the spectral parameter and V the amplitude of the barrier. It is well-knownthat the number of zeros of this coefficient represents the number of solitons produced by the square barrier initialcondition [10].It can be proved that the roots of this equation are purely imaginary. (This initial condition satisfies the single-lobeconditions of Klaus-Shaw potentials, from which it follows that the eigenvalues are purely imaginary [22]). Let usdefine η ≥ E = iη. Then, Eq. (3) becomes p − η cos (cid:16) V p − η L (cid:17) + η sin (cid:16) V p − η L (cid:17) = 0 . (4)We can verify that Eq. (4) does not have any roots (i.e., leads to no solitons in Eq. (1)) if V < π . (5)Furthermore, the condition to generate n solitons, i.e., so that Eq. (4) has n roots is(2 n − π < V L < (2 n + 1) π , (6)or, equivalently, the count of eigenvalues is given by2 π V L − < n < π V L + 12 . (7)The limit V → ∞ together with L → V L = const . In this instance, the number ofeigenvalues stays the same.
B. The Ablowitz-Ladik Model
We now turn to the integrable discretization of Eq. (1) and examine its dynamics. The one-dimensional integrable,discrete nonlinear Schr¨odinger model (so-called Ablowitz-Ladik (AL) model [16, 17, 18]) reads: i ˙ u n = −
12 ( u n +1 + u n − − u n ) − | u n | ( u n +1 + u n − ) . (8)In this case, there exists a Lax pair of linear operators [16, 17, 18] L n = Z + M n , (9) B n = (cid:18) z − z − (cid:19) D + 12 (cid:0) ZM n − Z − M n − (cid:1) − DM n M n − , (10)with the definitions for the matrices Z = (cid:18) z z − (cid:19) , D = (cid:18) − (cid:19) , and M n = (cid:18) U n − U ⋆n (cid:19) , (11)where z is the spectral parameter, and U n = U n ( t ) is a solution of the equation. These two operators (9) and (10)define the system of differential-difference equationsΨ n +1 = L n Ψ n , (12) i ddτ Ψ n = B n Ψ n , (13)for a complex matrix function Ψ n . Then, the compatibility condition of Eqs. (12) and (13) i ddτ Ψ n +1 = i (cid:18) ddτ Ψ m (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) m = n +1 (i.e., the Lax equation) becomes the AL model. U n = U n ( t ) is referred to as the potential of the AL eigenvalueproblem.For U n decaying rapidly at ±∞ , and for n → ±∞ , from Eq. (12) we have:Ψ n +1 ∼ Z Ψ n . We normalize this type of solutions as follows: let Ψ n denote the solution of Eq. (12) such thatΨ n ∼ Z n as n → + ∞ , and let Φ n be the solution of Eq. (12) such thatΦ n ∼ Z n as n → −∞ . Ψ n and Φ n are known as the Jost functions. Each of these forms a system of linearly independent solutions of theAL eigenvalue problem (12). These sets of solutions are inter-related by the scattering matrix S ( z ),Φ n = Ψ n S ( z ) . (14)The first column of this equation is given by(Φ ) n = S ( z )(Ψ ) n + S ( z )(Ψ ) n , (15)where (Φ ) n denotes the first column of Φ n . Similar definitions apply to (Ψ ) n and (Ψ ) n . Since Φ n ∼ Z n as n → −∞ ,then (Φ ) n ∼ z n (cid:18) (cid:19) . To obtain decay, (Φ ) n → n → −∞ , we demand | z | > . (16)On the other hand, Ψ n = ((Ψ ) n , (Ψ ) n ) ∼ (cid:18) z n z − n (cid:19) as n → ∞ . Therefore, if | z | > ) n → ∞ and (Ψ ) n →
0, as n → ∞ . Now, from Eq. (15) it follows that (Φ ) n → ∞ as n → ∞ , unless S ( z ) = 0 . We therefore seek solutions z , z , . . . , z N of the equation S ( z k ) = 0 , k = 1 , , . . . , N, (17)such that | z k | >
1. Then, (Φ ) n ( z k ) = S ( z k )(Ψ ) n ( z k ) , k = 1 , , . . . , N. From this, it follows that (Φ ) n ( z k ) decays at ±∞ :(Φ ) n ( z k ) → n → ±∞ . It is then said that (Φ ) n ( z k ) is an eigenfunction ( k = 1 , , . . . , N ), with corresponding eigenvalue z k .In the case of U n ( t = 0) = U δ n,n , the Jost function isΦ n = Z n − ( Z + M )Φ , for n ≥ , (18)with M M = (cid:18) U − U (cid:19) . (19)Furthermore, Ψ n = Z n for n ≥ is the identity matrix: Φ = I . Hence, Eq. (18) reads:Φ n = Ψ n Z − ( Z + M ) , for n ≥ . Its comparison with Eq. (14) leads to a scattering matrix S ( z ) = Z − ( Z + M ) . We thus obtain that the transmission coefficient S ( z ) = 1 , which never vanishes. This means that the one-site potential (i.e., a single-site initial condition) does not admitsolitonic solutions, independently of the amplitude U of initial excitation. This theoretical result has also beenconfirmed by numerical simulations for different values of U , always leading to dispersion of the solution (not shownhere). III. THRESHOLD CONDITIONS FOR THE NON-INTEGRABLE DNLS MODEL
We now turn to the non-integrable DNLS lattice that has the general form: i ˙ u n = − ǫ ∆ u n − | u n | σ u n , (20)where u n is a complex field (corresponding to the envelope of the electric field in optics, or the mean field wavefunctionin optical lattice wells in the BECs), ∆ u n = ( u n +1 + u n − − u n ) is the discrete Laplacian, and ǫ is the ratio of thetunneling strength to the nonlinearity strength. Using a scaling invariance of the equation, we can scale ǫ = 1, throughrescaling t → ǫt and u n → u n / √ ǫ . Importantly for our considerations, this model has an underlying Hamiltonian(which is also a conserved quantity in the dynamics) of the form: H = X n | u n +1 − u n | − σ + 1 | u n | σ +2 . (21)The only other known conserved quantity in the dynamics of Eq. (20) is the beam power (in Bose-Einstein condensates,the normalized number of atoms) P = P n | u n | .The specific problem that we examine here addresses the following question. Consider a “compactum” of mass u n = Aδ n, ; what is the critical value of A that is necessary for this single-site initial condition to excite a localizedmode? That such a threshold definitely exists is illustrated in Fig. 1. The case of the leftmost panel is subcritical,leading to the discrete dispersion of the initial datum. This follows the well-known t − / amplitude decay which isimplied by the solution of the problem in the absence of the nonlinearity: u n ( t ) = Ai n J n (2 t ) , (22)where J n is a Bessel function of order n . Notice that this can also be shown to be consistent with the findings of[23]. On the other hand, the rightmost panel shows a nonlinearity-dominated regime with the rapid formation of asolitary wave strongly localized around n = 0. In the intermediate case of the middle panel, the system exhibits along oscillatory behavior reminiscent of a separatrix between the basins of attraction of the two different regimes. n t0 5 10 15 20 25 30 35−50050 00.20.40.60.80 5 10 15 20 25 30 3500.51 | u , , | t −2 t | u , , | | u , , | n t0 50 100 150−10010 0123 | u , , | n t0 50 100 150−10010 024 FIG. 1: Subcritical case (left), critical (middle), and supercritical (right) initial single-site excitations on the lattice. In allcases, the initial condition is u n = Aδ n, , with A = 1 in the left, 2 in the middle, and 2 . | u n | forsites n = 0 , , t − decay for comparison. Clearly, the damped oscillation of the field modulus has an envelope of t − / . We argue herein that the essence of this separatrix lies in the examination of the stationary (localized) states ofthe model. Such standing wave solutions of the form u n = exp( i Λ t ) v n , which are exponentially decaying for v n as afunction of n , can be found for arbitrary frequency Λ (and arbitrary power P in one spatial dimension). This is awell-known result in one dimension; see e.g., Ref. [24]. On the other hand, the single-site initial condition discussedabove has an energy of: H ss = 2 A − σ + 1 A σ +2 , (23)where the subscript denotes single site. Figure 2 summarizes succinctly the power dependence of the energy for thesetwo cases for σ = 1. Both the energy of the stationary solutions as a function of their power and the single-siteenergy as a function of single-site power ( P ss = A ), are shown. Note that the two curves H ss ( P ss ) and H ( P ) do notintersect (except at the trivial point H = P = 0) since for ǫ = 0 single-site states are not stationary ones.The examination of Fig. 2 provides information on the existence of a sufficient condition for the formation of alocalized mode and on the dynamics of the single-site initial condition. Figure 2 shows that localized solutions existfor arbitrarily small values of the input power P , and that the energy of the localized states is negative. This impliesthat the crucial quantity to determine the fate of the process is the energy H and not the power P . The role of thepower will become more evident in the two-dimensional setting mentioned below. Moreover, if the system starts ata given point on the curve defined by H ss = 2 P ss − P σ +1 ss / ( σ + 1), due to conservation of total H and total P , itcan only end up in a stationary state in the quadrant H < H ss and P < P ss . That is to say, some of the initialenergy and power are typically “shed off” in the form of radiation (i.e., converted to other degrees of freedom whichis the only way that “effective dissipation” can arise in a purely Hamiltonian system), so that the initial conditioncan “relax” to the pertinent final configuration. As mentioned above, a localized solution with the same power asthat of the initial condition exists for arbitrary A . However , emergence of a localized mode occurs only for thoseinitial conditions whose core energy (i.e., the energy of a region around the initially excited state) is negative, afterthe profile is “reshaped” by radiating away both energy and power.Therefore, if H ss <
0, then the compactum of initial data will always yield a localized excitation: this inequalityprovides the sufficient condition for the excitation of solitary waves. The condition on the energy, in turn, providesa condition on the single-site amplitude that leads to the formation of solitary waves, namely, solitary waves alwaysform if
A > A ∗ with A ∗ = [2 ( σ + 1)] σ . (24)For the case of σ = 1 considered in Figs. 1 and 2, this amplitude value is A ∗ = 2 in agreement with our numericalobservations of Fig. 1. Whether an initial state with H ss > H FIG. 2: The blue solid line shows the energy H versus the power P of the discrete solitary wave solutions. Above it, the greensolid line shows the energy versus power of the initial condition, obtained from Eq. (23) and P ss = A . The horizontal linedenotes H = 0 and its intersection with the initial condition curve defines the single-site, initial amplitude A ∗ . The dynamicalevolution of three different supercritical initial states with A = 2 . A = 2 .
3, and A = 2 . One should make a few important observations here. Firstly, we note the stark contrast of the non-integrablediscrete model and both of its integrable (continuum and discrete) counterparts. In the (singular) continuum limit, itis possible to excite a single soliton or a multi-soliton depending on the barrier height and width. On the other hand,in the integrable discretization one-site excitation never leads to solitary wave formation, contrary to what is the casehere where either none or one solitary wave may arise, depending on the amplitude of the initial one-site excitation.Restricting our consideration to the non-integrable model, we observe that even though a localized solution withthe same power as that of the initial condition exists for arbitrary A , formation of a localized mode will always occuronly for A > A ∗ , as determined by the previous energetic criterion. Secondly, the answer that the formation alwaysoccurs for A > A ∗ generates two interconnected questions: what is the threshold for the formation of the localizedmode, and given an initial A > A ∗ , which one among the mono-parametric family of solutions will the dynamics of themodel select as the end state of the system ? [It should be noted in connection to the latter question that single-siteinitial conditions have always been found in our numerical simulations to give rise to at most a single-site-centeredsolution i.e., multipulses cannot be produced by this process.] Some examples of this dynamical process are illustratedin Fig. 2, where the energy and power of a few sites (typically 20-40) around the originally excited one are measuredas a function of time and are parametrically plotted in the H − P plane. As it should, the relevant curve startsfrom the H ss − P ss curve, and asymptotically approaches, as a result of the dynamical evolution, the H − P curveof the stationary states of the system. However, the relaxation process happens neither at fixed energy, nor at fixedpower . Instead, it proceeds through a more complex, dynamically selected pathway of loss of both H and P to relaxeventually to one of the relevant stationary states. This is shown for three different values of super-critical amplitudein Fig. 2 ( A = 2 .
1, 2 . . dH/dP ≈ const. However, the later stages ofrelaxation no longer preserve this constant slope. It is worthwhile to highlight here that the fact that the process offormation of a localized mode is neither equienergetic, nor does it occur at fixed power is something that has beenpreviously observed in the continuum version of the system in [20] (cf. with the discussion in p. 6095 therein).The dynamics of the system in the H − P space can be qualitatively understood by considering the frequenciesΛ associated with the instantaneous H , P values along the system trajectories. At each instantaneous energy H a frequency Λ H may be defined as the frequency of a single-site breather stationary state with energy H . Sucha frequency is unique as the stationary-state energy is a monotonically decreasing function of the frequency [24].Similarly, a frequency Λ P may be uniquely associated with the instantaneous power P , i.e., Λ P is the frequency of astationary state with power P . Note, however, that the stationary-state power is a monotonically increasing functionof Λ. At the final stationary state where the system relaxes Λ H = Λ P = − dH/dP . Hence, the system trajectoriesare such that Λ H increases (consequently the energy decreases) and Λ P decreases (the power decreases). The finalstationary state is reached when the two frequencies become equal, the point where they meet depending on their rateof change along the trajectory, i.e., on their corresponding “speeds” along the trajectory. Nevertheless, the precisemechanism of selection of the particular end state (i.e., of the particular “equilibrium Λ”) that a given initial statewill result in remains a formidable outstanding question that would be especially interesting to address in the future.We now turn to the two-dimensional variant of the above one-dimensional non-integrable lattice. Equation (20)remains the same, but for the two-dimensional field u n,m , and the discrete Laplacian becomes the five-point stencil∆ u n,m = ( u n +1 ,m + u n − ,m + u n,m +1 + u n,m − − u n,m ). For the initial condition u n,m = Aδ n, δ m, , the H ss − P ss curve is given by: H ss = 4 P ss − P σ +1 ss / ( σ + 1) . (25)In d -dimensions the first term would be 2 dP ss . We once again find the branch of standing waves of the equation.While this branch of solutions is also well-known [25], there are some important differences with the one-dimensionalcase. Firstly, a stable and an unstable branch of solutions exists (per the well-known Vakhitov-Kolokolov criterion); inthe wedge-like curve indicating the standing waves only the lower energy branch is stable. Furthermore, the maximalenergy of the solutions is no longer H = 0, but finite and positive (in fact, for σ = 1, it is H ≈ . H ss with H we find two solutions: one with A ≈ . A ≈ .
73. However,for the lower one there are no standing wave excitations with the corresponding power (this reveals the role of thepower in the higher dimensional problem). Hence, the relevant amplitude that determines the sufficient conditionin the two-dimensional case is the latter. Indeed, the case with A = 2 .
65 shown in the left panel of Fig. 3 showsspatio-temporal diffraction, while that of A = 2 .
75 in the middle panel illustrates robust localization, indicating thesuper-critical nature of the latter case. Notice that the pathway of the dynamics is presented in the right panel forthree super-critical values of A = 2 . A = 2 . A = 3 .
1. Once again, the dynamics commences on the H ss − P ss curve, as it should, eventually relaxing on the stable standing wave curve. As before, the initial dynamics follows aroughly constant dH/dP , but the relaxation becomes more complex at later phases of the evolution. H FIG. 3: The left panel shows the spatio-temporal evolution of | u n,m ( t ) | for the sub-critical case with amplitude A = 2 . A = 2 .
75. Finally the right panel is analogous to Fig. 2, showing the H − P diagram for the solution branch (wedge-like blue solid line), the H ss − P ss graph of the initial conditions (green solidline) and the trajectories of three super-critical cases in H − P space for A = 2 .
75, 2 . . H ≈ . IV. CONCLUSIONS AND FUTURE CHALLENGES
The above study has examined the presence of a sharp crossover between the linear and nonlinear dynamics of aprototypical dynamical lattice model such as the discrete nonlinear Schr¨odinger equation. This crossover has alreadybeen observed in media with the focusing nonlinearity [7] (as considered here). It has also been observed very recentlyin media with the defocusing nonlinearity [8]. The latter can be transformed into the former under the so-calledstaggering transformation u n = ( − n w n , where w n is the field in the defocusing case. As a result, the solitary wavesof the defocusing problem discussed in [8] will be “staggered” (i.e., of alternating phase between neighboring sites), yetthe phenomenology discussed above will persist. The crossover was quantified on the basis of an energetic comparisonof the initial-state energy with the branch of corresponding stable localized solutions “available” in the model. Asufficient, but not necessary, condition for the excitation of a localized mode based on the initial-state, single-siteamplitude was discussed, and it was successfully tested in numerical simulations. Similar findings were obtained inthe two-dimensional analog of the problem: a crossover behaviour dictated by the energy was found, but the crossoverwas also affected by the power and its excitation thresholds. Furthermore, these results were contrasted with the caseof the continuum version of the model, where depending on the strength of the excitation, also multi-solitons can beobtained and with the integrable discrete Ablowitz-Ladik model where no single-site excitation can produce a solitarywave, independently of the excitation amplitude.However, a number of interesting questions emerge from these findings that are pertinent to future studies. Perhapsthe foremost among them concerns how the dynamics “selects” among the available steady-state excitations withenergy and power below that of the initial condition the one to which the dynamical evolution leads. According toour results, this evolution is neither equienergetic, nor power-preserving (see also [20]), hence it would be extremelyinteresting to identify the leading physical principle which dictates it. Another question generalizing the more principalone asked herein concerns the excitation of multiple sites (possibly at the same amplitude), starting with two, andinquiring the nature of the resulting state of the system. On a related note, perhaps this last question is most suitableto be addressed first in the context of the integrable model where the formulation presented herein can be generalizedand yield definitive answers for the potential excitation of solitons. Studies along these directions are currently inprogress and will be reported in future publications. V. ACKNOWLEDGEMENTS
YD gratefully acknowledges discussions with L. Isella. PGK gratefully acknowledges support from grants NSF-DMS-0505663, NSF-DMS-0619492 and NSF-CAREER. [1] S. Aubry, Physica D , 201, (1997); S. Flach and C.R. Willis, Phys. Rep. , 181 (1998); D. Hennig and G. Tsironis,Phys. Rep. , 333 (1999);[2] D. N. Christodoulides, F. Lederer and Y. Silberberg, Nature , 817 (2003); Yu. S. Kivshar and G. P. Agrawal,
OpticalSolitons: From Fibers to Photonic Crystals , Academic Press (San Diego, 2003).[3] V.V. Konotop and V.A. Brazhnyi, Mod. Phys. Lett. B , 627 (2004); O. Morsch and M. Oberthaler, Rev. Mod. Phys. , 179 (2006); P.G. Kevrekidis and D.J. Frantzeskakis, Mod. Phys. Lett. B , 173 (2004).[4] M. Sato, B. E. Hubbard, and A. J. Sievers, Rev. Mod. Phys. , 137 (2006).[5] M. Peyrard, Nonlinearity , R1 (2004).[6] P.G. Kevrekidis, K.O. Rasmussen, and A.R. Bishop, Int. J. Mod. Phys. B , 2833 (2001).[7] H.S. Eisenberg, Y. Silberberg, R. Morandotti, A.R. Boyd and J.S. Aitchison, Phys. Rev. Lett. , 3383 (1998).[8] M. Matuszewski, C.R. Rosberg, D.N. Neshev, A.A. Sukhorukov, A. Mitchell, M. Trippenbach, M.W. Austin, W. Kro-likowski and Yu.S. Kivshar, Opt. Express , 254 (2006).[9] C. Sulem and P.L. Sulem, The Nonlinear Schr¨odinger Equation , Springer-Verlag (New York, 1999).[10] V.E. Zakharov and A.B. Shabat, Soviet Physics-JETP , 62 (1972).[11] J.W. Fleischer, T. Carmon, M. Segev, N.K. Efremidis and D.N. Christodoulides, Phys. Rev. Lett. , 023902 (2003); J.W.Fleischer, M. Segev, N.K. Efremidis and D.N. Christodoulides, Nature , 147 (2003).[12] H. Martin, E.D. Eugenieva, Z. Chen and D.N. Christodoulides, Phys. Rev. Lett. , 123902 (2004).[13] D.N. Neshev, T.J. Alexander, E.A. Ostrovskaya, Yu.S. Kivshar, H. Martin, I. Makasyuk and Z. Chen, Phys. Rev. Lett. , 123903 (2004).[14] J.W. Fleischer, G. Bartal, O. Cohen, O. Manela, M. Segev, J. Hudock and D.N. Christodoulides, Phys. Rev. Lett. ,123904 (2004).[15] J.W. Fleischer, G. Bartal, O. Cohen, T. Schwartz, O. Manela, B. Freeman, M. Segev, H. Buljan, and N.K. Efremidis,Optics Express , 1780 (2005).[16] M.J. Ablowitz and J.F. Ladik, J. Math. Phys. , 598 (1975).[17] M.J. Ablowitz and J.F. Ladik, J. Math. Phys. , 1011 (1976).[18] M.J. Ablowitz, B. Prinari and A.D. Trubatch, Discrete and Continuous Nonlinear Schr¨odinger Systems , Cambridge Uni-versity Press (Cambridge, 2004).[19] N.N. Akhmediev and A. Ankiewicz,
Solitons: Nonlinear pulses and beams , Chapman & Hall (London 1997).[20] N.N. Akhmediev, A. Ankiewicz, and R. Grimshaw, Phys. Rev. E , 6088 (1999).[21] B. Rumpf, Phys. Rev. E , 016609 (2004).[22] M. Klaus and J. K. Shaw, Phys. Rev. E , 036607 (2002).[23] A. Stefanov and P.G. Kevrekidis, Nonlinearity , 1841 (2005).[24] M. Johansson and S. Aubry, Phys. Rev. E , 5864 (2000).[25] P. G. Kevrekidis, K. Ø. Rasmussen, and A. R. Bishop, Phys. Rev. E , 2006-2009 (2000).[26] S. Flach, K. Kladko and R.S. MacKay, Phys. Rev. Lett. , 1207 (1997).[27] M.I. Weinstein, Nonlinearity , 673 (1999).[28] P.G. Kevrekidis, K.Ø. Rasmussen and A.R. Bishop, Phys. Rev. E61