Radiationless traveling waves in saturable nonlinear Schrödinger lattices

Radiationless Traveling Waves in Saturable Nonlinear Schro ¨dinger Lattices T. R. O. Melvin, 1 A. R. Champneys, 1 P. G. Kevrekidis, 2 and J. Cuevas 3 1 Department of Engineering Mathematics, University of Bristol, BS8 1TR, United Kingdom 2 Department of Mathematics and Statistics, University of Massachusetts, Amherst, Massachusetts 01003-4515, USA 3 Grupo de Fı ´sica No Lineal, Departamento de Fı ´sica Aplicada I, Escuela Universitaria Polite ´cnica, Virgen de A ´frica, 7, 41011 Sevilla, Spain (Received 30 March 2006; published 18 September 2006) The long-standing problem of moving discrete solitary waves in nonlinear Schro ¨dinger lattices is revisited. The context is photorefractive crystal lattices with saturable nonlinearity whose grand-canonical energy barrier vanishes for isolated coupling strength values. Genuinely localized traveling waves are computed as a function of the system parameters for the first time. The relevant solutions exist only for finite velocities. DOI: 10.1103/PhysRevLett.97.124101 PACS numbers: 42.65.Tg, 05.45.Yv, 63.20.Pw Recently, the topic of discrete solitons (intrinsic localized modes) in photorefractive materials has received significant attention; see [1] for a review. This interest was initialized by the experimental realization of twodimensional periodic lattices in photorefractive crystals [2] in which solitons were observed [3,4]. Further work has revealed a wealth of additional coherent structures such as dipoles, quadrupoles, soliton trains, vector, necklace, and ring solitons; see, e.g., [5–8]. Since photorefractive materials feature the so-called saturable nonlinearity, these results for periodic lattices have spawned a parallel interest in genuinely discrete saturable nonlinear lattices [9–11]. One particularly intriguing result of these studies is that the stability properties of the localized modes are substantially different from their regular discrete nonlinear Schro ¨dinger (DNLS) analogs. In DNLS, it is well known [12–14] that site-centered localized modes are always stable, while intersite-centered modes are unstable (and are stabilized only in the continuum limit where these two branches degenerate into the well-known continuum sech-soliton of the integrable cubic NLS). For the photorefractive nonlinearity (i.e., the so-called Vinetskii-Kukhtarev model originally proposed in Ref. [15] and revisited in Refs. [9–11]), depending on the coupling strength, the intersite-centered modes may have lower energy than their on site counterparts. Hence, one should expect that the ensuing sign reversal of the so-called Peierls-Nabarro (PN) energy barrier EEIS EOS (where the subscripts denote intersite and on site, respectively) should cause an exchange between the linear stability properties of the two modes. A related, even more fundamental, question in nonlinear lattice models of DNLS type is whether exponentially localized self-supporting excitations that move with a constant wave speed can exist, the so-called moving discrete solitons. For continuum models, this question is in some sense trivial since the equations posed in a moving frame remain of the same fundamental type. Yet for lattices, passing to a moving frame leads to so-called advance-delay equations that are notoriously hard to analyze. One recent (negative) result in this direction [16] for the so-called Salerno model shows that, starting from the integrable Ablowitz-Ladik equation, mobile discrete solitary waves acquire nonvanishing tails as soon as parameters deviate from the integrable limit. Hence, exponentially localized fundamental (single-humped) moving discrete solitons cannot be constructed. While these results settle a longstanding controversy—see, e.g., [17–19]—they are also somewhat unsatisfactory since they do not give conditions under which moving discrete solitons might exist for generic, nonintegrable lattices. In the arguably simpler problem of kinks in FrenkelKontorova lattices (and, more generally, in discrete KleinGordon models), genuinely traveling fundamental (‘‘charge one’’) topological solitons are known not to occur unless there is a ‘‘competing’’ nonlinearity that causes vanishings of the PN barrier; see, e.g., [20,21]. Drawing an analogy from this, for (the more complicated, complex field) DNLS models with pure cubic nonlinearity, the PN barrier never disappears, so one might not expect genuinely localized moving solitons (in some sense anticipating the above-mentioned negative result). However, the present saturable model has a nonlinear term in its denominator, which, under Taylor expansion, is effectively cubic for low intensities but exhibits saturation for large intensities. This effective ‘‘competition’’ between nonlinear behaviors leads to the possibility of the vanishing of the PN barrier and, hence, suggests this model as a good candidate for finding moving discrete solitons. Such a search is the theme of this Letter. We note that such a connection was also partially discussed in Ref. [22], based on the properties of the stationary localized modes. We start by examining the existence of stationary localized modes as the strength of intersite coupling is varied, initializing our computations from the anticontinuum limit, where explicit solutions are available. A linear stability analysis of the solutions is carried out, and we find that the regime of stability-instability alternation does not coincide with the vanishings of the PN barrier but rather with the vanishings of the grand-canonical free-energy of the sysPRL 97, 124101 (2006) PHYSICAL REVIEW LETTERS week ending 22 SEPTEMBER 2006 0031-9007=06=97(12)=124101(4) 124101-1 ©2006 The American Physical Society tem GEP, where is the frequency (chemical potential) of the solution and Pis its l2norm. After identifying this sequence of transparent points, we proceed to obtain, for the first time in discrete nonlinear Schro ¨dinger type models, genuinely localized, single-humped, uniformly traveling solutions in their vicinity. This is done by examining traveling solutions as embedded solitons [23], such that, while resonant with the linear spectrum in the traveling frame, the tail amplitude can vanish for isolated, definitely nonzero values of the speed, for a given set of system parameters. This is reminiscent of recent work in the Klein-Gordon context [24], where the so-called Stokes constant is shown to vanish for isolated values of the speed, indicating the presence of traveling localized solutions. The model of interest, as discussed above, will be a discrete lattice with photorefractive nonlinearity, namely, i_ un"2un 1junj2un:(1) Here "denotes the coupling strength between adjacent sites and is the coefficient of the nonlinearity (proportional to the voltage in photorefractive crystals) and can be scaled out for our purposes; 2unun1un12un denotes the discrete Laplacian. We will therefore set 1 in what follows. We seek stationary solutions of the form uneitvn, starting from the so-called anticontinuum (AC) limit "0. In that limit, the principal solutions of interest are the on site solution with un0  1=1 p(all other sites bearing vanishing amplitudes) and the intersitecentered solution with un0un01  1=1 pand zero amplitudes elsewhere. We initialize these exact profiles from the AC limit and subsequently use continuation in "to solve the stationary equations for the wave profile via Newton’s method for finite ". We also examine the linear stability of the obtained modes using linearization around a solution profile v0 nof the form uneitv0 n anet bne?t and resolving the ensuing matrix eigenvalue problem to Ofor the eigenvalue and the eigenvector pair (an,b? n) (where ?stands for complex conjugation). In what follows, we will also use the (Hamiltonian) energy of solutions EX n "jun1unj2log1junj2;(2) and the l2norm, PPnjunj2. Our stationary calculations as a function of the intersite coupling (rather than of the frequency as in Ref. [9]) are shown in Fig. 1. Figures 1(a) and 1(b) show the key imaginary eigenvalues of the problem, while Fig. 1(c) shows the maximal real eigenvalues. A solid line is used in the figures to represent the intersite-centered (IS) mode and a dashed line the on site-centered (OS) mode. Observe that for small values of "the OS mode is stable while the intersite mode is unstable. Consider the OS mode. At " 0:1, an eigenvalue emerges from the band edge of the continuous spectrum at iand approaches the origin of the spectral plane. This feature also occurs in the pure cubic DNLS but the eigenvalue only arrives at the origin of the spectral plane as "!1(the continuum limit). Here the eigenvalue arrives at the origin for "0:25, becomes real for 0:25 <"<0:445, then becomes imaginary again for 0:445 <"<0:7, then real again, etc. Furthermore, this alternation of stability occurs hand in hand with the alternation of stability of the IS mode, which starts out unstable, becomes stable, then unstable again, with the transitions occurring at precisely the same "values. Note from the figure that stability alternation occurs, not at the points of vanishing of the PN energy barrier, but when the grandcanonical energy barrier of the model GEP between the IS and OS modes vanishes. This provides an alternative viewpoint to the results of Ref. [9], which is also rather natural from the point of view of statistical mechanics for a system with a conserved ‘‘particle number’’ (l2 norm); see, e.g., [25]. At these very points, the relevant eigenvalue (pertaining to translation) crosses zero, instantaneously restoring an effective translational invariance in the model. Intuitively, one might expect that the transparent points represent parameter values at which genuinely traveling localized solutions of the photorefractive lattice might bifurcate with zero wave speed. Hence, we shall look for solutions to Eq. (1) near such points using a traveling wave substitution unt ncteiknt. This leads to the advance-delay equation in the traveling coordinate z nct ic 0z2" z" z1eik  z1eik 1 1j zj2 z;(3) where 0denotes differentiation with respect to z. Equation (3) is rotationally invariant; therefore, the transformation z ~ zeikz can be used to obtain the same 0 0.5 1 1.5 2 100 λi ε 0 0.5 1 1.5 2 0 0.2 0.4 ε λi A B 0 0.5 1 1.5 2 10−5 100 ε Max(λr) 0 0.5 1 1.5 2 −20 −10 0 log(|∆ E|), log(|∆ G|) ε C D FIG. 1 (color online). (a) and (b) show the key imaginary eigenvalues for the on site (dashed line) and the intersite (solid line) mode as a function of "[(a) on a linear scale, (b) on an exponential scale]. The band edge of the continuous spectrum is at 0:5. From there bifurcates an eigenmode for ">0:1(for the on site case) which arrives at the origin of the spectral plane for "0:25 and becomes real. On the contrary, for the intersite solution a previously real eigenmode exits as imaginary for "> 0:25. These modes alternate again for "0:445,"0:7, etc. The real parts of the corresponding eigenvalues are shown in (c). (d) shows logjEj between on site and intersite modes and by a dashed line the quantity logjGj, where GEP(see text for details). PRL 97, 124101 (2006) PHYSICAL REVIEW LETTERS week ending 22 SEPTEMBER 2006 124101-2 equation as above for ~ and ~ kc . Therefore, k 0will be used hereafter in (3) without loss of generality. The linear spectrum is obtained by substitution of z eiz into the linear part of Eq. (3) c 14"sin2 2:(4) It is interesting to note that the asymptotic behavior of the two sides of Eq. (4) for small is such that a resonance with the linear spectral bands is unavoidable. However, in the hope of obtaining a localized solution, such resonances need to be minimized, and, hence, we need to investigate parameter regions where Eq. (4) has only one root, >0. This can be efficiently done [21] by considering the conditions for double roots, namely, c2"sin, which allow us to distinguish among regions with different root multiplicity. The results are plotted in Fig. 2, from which it can be seen that, for each ">0, the number of spectral bands increases to infinity as c!0. This result serves to counter the expectation that, near the transparency points, traveling solutions of small speed may exist. The best hope for finding genuinely localized solutions is then within parameter regions of Fig. 2in which there is a single resonance with linear waves. Such solutions would represent embedded soliton structures [21,23] for which the radiation mode component exactly vanishes in the tail of the solitary wave. At best, this could happen in a codimension-one set in parameter space, that is, along discrete curves lying within the single resonance bands of Fig. 2. Solutions of this type are sought using a pseudospectral method; see [18,20]. To this end, we set zX N j1 ajcos!jzibjsin!jz;(5) with !jj L, to transform (3) into a system of algebraic equations in the long finite interval L=2;L=2.aj;b j2 Rare the coefficient of the Fourier series, for which the ensuing algebraic equations are solved at the collocation points zmLm 2N1. The solution is obtained by means of the Powell hybrid method [26], with an error tolerance of 1013. From the resulting aj;b j’s, zcan be reconstructed and numerically continued using AUTO [27] with an error tolerance of 1010 to investigate the effect of varying parameters on the solution shape and its tail amplitude. Generally, the thus-obtained solutions will be weakly nonlocal ones with nonzero oscillatory tails. However, appending an additional tail condition for the solution, yielding a signed measure of the amplitude of the tail, we can obtain truly localized solutions as isolated zeros of this function . Given our choice of phase in the substitution (5) that the real component of is odd about zL=2(and even for the imaginary component), a good choice for such a tail amplitude measure is Im L 2. Continuation of weakly nonlocal solitary wave solutions as "varies is shown as a function of 1=" in Fig. 3(top left). Note the topological distinction between the S-shaped branches that contain true zeros of and the U-or n-shaped branches that do not [28]. Once the isolated points where the tail disappears are found, these can be continued in any pair of desired parameters, e.g., "; c,by appending the condition 0to the numerical problem. Thus, we can construct the parametric diagram of the existence of such genuinely traveling (with constant prescribed speed), yet truly localized (i.e., radiationless due to 0) solutions. A regular sequence of branches so obtained are found as "increases, the first three of which are shown in the top right panel in Fig. 3. Notice that each branch terminates (as it should) at the edges of the multispectral bands. Nevertheless, at the low-cend the branches are close to straight lines, which, if continued to zero wave speed, would hit the transparent points. Thus, while there is no actual bifurcation as such of intrinsic localized modes from these transparent points, their presence apparently has a strong influence. The two rightmost upper panels in Fig. 3show solution profiles on the first (lowest ") branch of localized solutions. Note how the profile delocalizes when the edge of the multiresonance band is reached. Any attempt to continue zero-tail solutions into the band resulted in a lack of convergence for large enough Land N. The lower two panels show the result of direct integration of a representative solution on branch I, clearly illustrating their genuinely traveling nature and indicating that they are stable against the perturbation introduced by numerical discretization. In conclusion, to our knowledge, this is the first time genuinely localized traveling waves have been found in a discrete lattice of nonlinear Schro ¨dinger type. While it is true that our construction has relied on numerical continuation, we have demonstrated the robustness of the solutions we have obtained and have provided a new rationale as to why they should exist. This rationale is based on the 0 1 2 3 4 0 1 2 3 4 ε c 010 −1 0 4 λ0 6 0 15 λ 050 − 10 0 30 λ 00.1 0.2 0.3 0.4 0 0.05 0.1 0.15 0.2 ε c 13 3 1 57 9 (i) (ii) (iv) (iii) FIG. 2. Intersections of the single and double root conditions (see text) for 0:5and varying cand ". The shaded area shows where there is more than one branch of linear waves. Values indicate the number of roots of (4). Subplots show leftand right-hand sides of (4) for "; c1;3, one root (i); 3;3, three roots (ii); and 3:5;0:6, seven roots (iii). (iv) displays the overlapping of the bands for small "and c; only the first six bands have been shown. PRL 97, 124101 (2006) PHYSICAL REVIEW LETTERS week ending 22 SEPTEMBER 2006 124101-3 existence of transparent points for stationary localized modes when a generalized Peierls-Nabarro barrier vanishes. However, we have also shown that this is not enough; one also needs the wave speed to be sufficiently large to overcome the multispectral bands. Within each band, exponentially localized solutions occur as embedded solitons which correspond to regular zeros of an appropriately defined measure of the tail amplitude. We anticipate that the principal features of the model problem studied (competing nonlinear terms for the existence of transparent points and sufficiently high wave speeds to overcome the multispectral bands) should be more general features of a variety of similar Hamiltonian lattices. Indeed, we have found similar results for a pure cubic-quintic nonlinearity. More realistic models of photorefractive crystals should be expected to behave similarly, giving the possibility of experimental observation of the waves we have constructed. There, the robust, radiationless nature of the new moving pulses may be particularly appealing. Another relevant question is the extension of the present ideas to two-dimensional settings [29], where deciding which directions of wave propagation within the lattice can lead to radiationless pulses will also be of interest. [1] J. W. Fleischer, Opt. Express 13, 1780 (2005). [2] N. K. Efremidis et al., Phys. Rev. E 66, 046602 (2002). [3] J. W. Fleischer et al., Phys. Rev. Lett. 90, 023902 (2003). [4] H. Martin et al., Phys. Rev. Lett. 92, 123902 (2004). [5] J. Yang et al., Phys. Rev. Lett. 94, 113902 (2005). [6] X. Wang, Z. Chen, and P. G. Kevrekidis, Phys. Rev. Lett. 96, 083904 (2006). [7] D. N. Neshev et al., Phys. Rev. Lett. 92, 123903 (2004). [8] J. W. Fleischer et al., Phys. Rev. Lett. 92, 123904 (2004). [9] L. Hadz ˇievski et al., Phys. Rev. Lett. 93, 033901 (2004). [10] M. Stepı ´cet al., Phys. Rev. E 69, 066618 (2004). [11] J. Cuevas and J. C. Eilbeck, nlin.PS/0501050. [12] P. G. Kevrekidis, K. Ø. Rasmussen, and A. R. Bishop, Int. J. Mod. Phys. B 15, 2833 (2001). [13] D. N. Christodoulides, F. Lederer, and Y. Silberberg, Nature (London) 424, 817 (2003). [14] T.Kapitula and P. Kevrekidis, Nonlinearity 14, 533 (2001). [15] V. L. Vinetskii and N. V. Kukhtarev, Fiz. Tverd. Tela 16, 3714 (1974) [Sov. Phys. Solid State 16, 2414 (1975)]. [16] J. Go ´mez-Garden ˜es, L. M. Florı ´a, M. Peyrard, and A. R. Bishop, Chaos 14, 1130 (2004). [17] Nonlinear Coherent Structures in Physics and Biology, edited by M. Remoissenet and M. Peyrard (SpringerVerlag, Berlin, 1991). [18] D. B. Duncan et al., Physica (Amsterdam) 68D, 1 (1993). [19] M. J. Ablowitz, Z. H. Musslimani, and G. Biondini, Phys. Rev. E 65, 026602 (2002). [20] A. V. Savin, Y. Zolotaryuk, and J. C. Eilbeck, Physica (Amsterdam) 138D, 267 (2000). [21] A. Aigner, A. R. Champneys, and V. M. Rothos, Physica (Amsterdam) 186D, 148 (2003). [22] A. Maluckov, Lj. Hadzievski, and M. Stepic, Physica (Amsterdam) 216D, 95 (2006). [23] J. Yang, B. A. Malomed, and D. J. Kaup, Phys. Rev. Lett. 83, 1958 (1999). [24] O. F. Oxtoby, D. E. Pelinovsky, and I. V. Barashenkov, Nonlinearity 19, 217 (2006). [25] J. L. Lebowitz, H. A. Rose, and E. R. Speer, J. Stat. Phys. 50, 657 (1988); K. O ¨. Rasmussen et al., Phys. Rev. Lett. 84, 3740 (2000). [26] M. J. D. Powell, Numerical Methods for Nonlinear Algebraic Equations (Gordon and Breach, New York, 1970). [27] E. J. Doedel et al.,AUTO97 continuation and bifurcation software for ordinary differential equations, available via: ftp://ftp.es.concordia.ca/directory/doedel/auto. [28] A. R. Champneys, J.-M. Vanden-Broeck, and G. J. Lord, J. Fluid Mech. 454, 403 (2002). [29] R. Vicencio and M. Johansson, Phys. Rev. E 73, 046602 (2006). III II I 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 1.2 ε c III III 0.2 0.4 0.6 −40 −20 0 20 40 0 1 2 c z Re(ψ) 0.2 0.4 0.6 −40 −20 0 20 40 −1 0 1 c z Im(ψ) n t 20 40 60 80 0 20 40 60 80 100 0 20406080 0 0.2 0.4 0.6 0.8 1 n |u| FIG. 3 (color online). (Top left) Continuation of weakly nonlocal solitary waves for various values of 1=" showing three zeros in  for c0:7,0:5,L60. Zeros of at "0:76, 1.02, 1.36. (Top second from left) Continuation of the 3 zeros of , varying " and cwith 0:5. Circles on the c0axis indicate the transparency points. The top rightmost two panels show the continuation of the (real and imaginary parts of the) solution of branch I for different values of the speed. A typical example of direct integration of the solution of branch I with "0:911 396 and c0:894 153 (bottom panels). The space time contour plot of the solution modulus (bottom left panel) and the modulus before (solid line) and after (dashed line) 100 time steps (bottom right panel) are shown. PRL 97, 124101 (2006) PHYSICAL REVIEW LETTERS week ending 22 SEPTEMBER 2006 124101-4