Exact numerical solutions for dark waves on the discrete nonlinear Schrödinger equation

Exact numerical solutions for dark waves on the discrete nonlinear Schrödinger equation Bernardo Sánchez-Rey Departamento de Física Aplicada I, E. U. P., Universidad de Sevilla, Virgen de África 7, 41011 Sevilla, Spain Magnus Johansson Department of Physics and Measurement Technology, Linköping University, SE-581 83 Linköping, Sweden 共Received 12 November 2004; published 30 March 2005兲 In this paper we study numerically existence and stability of exact dark waves on the 共nonintegrable兲 discrete nonlinear Schrödinger equation for a finite one-dimensional lattice. These are solutions that bifurcate from stationary dark modes with constant background intensity and zero intensity at a site, and whose initial state translates exactly one site each period of the internal oscillations. We show that exact dark waves are characterized by an oscillatory background whose wavelength is closely related with the velocity. Faster dark waves require smaller wavelengths. For slow enough velocity dark waves are linearly stable, but when trying to continue numerically a solution towards higher velocities bifurcations appear, due to rearrangements in the oscillatory tail in order to make possible a decreasing of the wavelength. However, in principle, one might control the stability of an exact dark wave adjusting a phase factor which plays the role of a discreteness parameter. In addition, we also study the regimes of existence and stability for stationary discrete gray modes, which are exact solutions with phase-twisted constant-amplitude background and nonzero minimum intensity. Also such solutions develop envelope oscillations on top of the homogeneous background when continued into moving phase-twisted solutions. DOI: 10.1103/PhysRevE.71.036627 PACS number共s兲: 05.45.Yv, 42.65.Wi, 03.75.Lm, 63.20.Ry I. INTRODUCTION There is a large current interest in effects arising from competition between nonlinearity and dispersion in spatially periodic systems, which in many cases can be well captured by simple nonlinear lattice models 关1兴. In particular, much recent experimental activity has been devoted to studies of arrays of nonlinear optical waveguides 关2兴, as well as to studies of Bose-Einstein condensates 共BECs兲in optical lattices 关3兴. Under certain conditions 共notably strong periodic potentials and, for BECs, a large number of atoms per well兲, the dynamics can in both cases 共Refs. 关4,5兴, respectively兲approximately be described by the discrete nonlinear Schrödinger 共DNLS兲model 关Eq. 共1兲below, see, e.g., Ref. 关6兴 for a review of its properties and applications兴. As predicted by the DNLS model, the existence of self-localized discrete solitons has been experimentally confirmed in both contexts 共Refs. 关7,8兴, respectively兲. These are particular examples of a generic class of time-periodic, spatially localized modes in nonlinear lattices, the so-called “intrinsic localized modes” or “discrete breathers” 关9兴. An important consequence of the discreteness 共or spatial periodicity兲in these systems is, that self-localization is possible also when the nonlinearity in itself is defocusing 共repulsive interactions兲. This is a result of the band gap structure of the dispersion relation for linear waves, and the anomalous dispersion 共or diffraction, in the case of spatial solitons兲at the Brillouin zone boundary. In particular, in the DNLS description the attractive and repulsive cases are mathematically equivalent through a “staggering” transformation 共−1兲j, corresponding to reversing the phases in every second potential well. Thus, in such a situation discrete solitons 共breathers兲have exponentially decaying staggered tails, and have been observed experimentally in optically induced nonlinear photonic lattices 关10兴, waveguide arrays 关11兴,as well as BECs in periodic potentials 关8兴. On the other hand, it is well known 关12兴that without the periodic potential 共continuum system兲, a defocusing nonlinearity generally leads to the existence of dark solitons, consisting of a localized dip in a homogeneous background intensity 关like the dark-soliton solution of the defocusing integrable nonlinear Schrödinger 共NLS兲equation 关13兴兴. Such modes are ubiquitous in nonlinear optics 关12兴, and also observed in BECs with repulsive interatomic interaction 共e.g., Ref. 关14兴兲. When adding a periodic potential, it is possible to generate stationary discrete dark solitons 共breathers兲not only in the standard defocusing case, but also for focusing 共attractive兲nonlinearity where staggered dark solitons may exist, analogously to the staggered bright solitons in the defocusing case. Such staggered discrete 共immobile兲dark solitons, predicted from the DNLS model 关16兴, have indeed been experimentally observed in waveguide arrays 关15,11兴. For BECs in optical lattices, dark solitons have not yet been experimentally reported, but theoretically analyzed in Ref. 关17兴. However, as was found first numerically for the DNLS model in Ref. 关16兴and later confirmed in Ref. 关18兴, the discreteness may induce oscillatory instabilities for the stationary dark modes. It was found 关18兴, that for infinite chains the strength of these instabilities decay in an exponential-like way approaching the continuum NLS limit, while for finitesize systems stabilization due to boundary effects may occur. Typically, these instabilities lead to a spontaneous motion of the dark mode, together with some radiation 关16,18兴.Itis therefore an interesting question, whether such moving dark solitons 共breathers兲may exist as exact solutions in the DNLS model, and if so, whether they may be stable in regimes where the stationary dark modes are unstable. Addressing these issues is the main purpose of the present paper. In integrable models, such as the Ablowitz-Ladik discretization of the NLS equation, traveling dark solitons exist and can be obtained exactly analytically 关19兴. The evidence that PHYSICAL REVIEW E 71, 036627 共2005兲 1539-3755/2005/71共3兲/036627共11兲/$23.00 ©2005 The American Physical Society036627-1 moving dark breathers may exist as exact solutions in a nonintegrable lattice was given by Feddersen, who showed numerically 关20兴, for the DNLS as well as for the discrete Davydov equations, that traveling dark waves bifurcated from stationary solutions. However, the solutions he found all had a wide envelope and were close to continuum NLS dark solitons. Attempting a continuation towards narrower modes, he found that it always stopped at some point, but he did not analyze this scenario in detail. On the other hand, much recent attention has been drawn to the study of mobile discrete 共bright兲breathers in general oscillator chains, and it has been found, numerically 关21,22兴 as well as mathematically 关23兴, that such solutions may exist but typically as nonlocalized solutions with a smallamplitude oscillatory tail. Apreliminary calculation of one of the present authors 共M.J.兲关24兴showed, that generally also the tails of numerically exact moving dark DNLS solitons did not have a constant amplitude, but exhibited smallamplitude envelope oscillations 共in Ref. 关25兴, such solutions were termed “dark nanopterons”兲. In addition, it has also been preliminary reported 关26兴, that also stationary gray modes, with nonzero minimum intensity and phase twisted background 共cf. gray solitons in continuum NLS models 关12兴兲, may exist in the DNLS model. In this paper, we will provide extensive numerical confirmation of these preliminary results. By performing numerical continuations from the stationary dark DNLS soliton, we observe the following main effects of discreteness: 共i兲An exact moving discrete dark wave generically has envelope oscillations around the constant-amplitude background and 共ii兲gray solitons with phase-twisted constant-amplitude background exist as stationary solutions trapped by the lattice potential. It should also be emphasized that, due to the generic nature of the DNLS equation, similar results should be expected also for dark breathers in other types of anharmonic lattice models. Stationary dark breathers have been discussed as well for Klein-Gordon 关25,27,28兴, as for Fermi-PastaUlam 共FPU兲models 关29,30兴, and the validity of the DNLS approximation for weak coupling and small-amplitude oscillations explicitly confirmed numerically. The outline of this paper is as follows. In Sec. II we introduce the DNLS model, and recapitulate briefly some properties of its stationary dark modes from 关18兴. In Sec. III we describe our numerical results for existence and stability of exact traveling dark waves; Sec. III A discusses the continuation versus velocity, and Sec. III B the continuation versus internal frequency, which essentially plays the role of discreteness parameter. Section IV discusses properties of gray solitons with a phase gradient, which may be either stationary or moving. Finally, in Sec. V some concluding remarks are made. II. STATIONARY DARK MODES Our starting point is the 共attractive兲DNLS equation i ␺ ˙j+兩 ␺ j兩2 ␺ j+C共 ␺ j+1 + ␺ j−1 −2 ␺ j兲=0, 共1兲 where ␺ j共t兲is the complex amplitude of the oscillator at site jon a lattice, and Cis a coupling constant. Without loss of generality, we assume C⬎0. We consider a finite lattice 共j =1,2,...,N兲and adopt, as usual, periodic boundary conditions: ␺ j+N共t兲= ␺ j共t兲. Let us assume that ␺ j共t兲is a solution of the DNLS Eq. 共1兲 and define ␺ j ⬘共t兲= ␺ j共t兲ei ␻ 0t.共2兲 By direct substitution into 共1兲it is easy to check that ␺ j ⬘共t兲 verifies i ␺ ˙j ⬘+兩 ␺ j ⬘兩2 ␺ j ⬘+C共 ␺ j+1 ⬘+ ␺ j−1 ⬘−2 ␺ j ⬘兲+ ␻ 0 ␺ j ⬘=0. 共3兲 As a first step we are going to study stationary dark solutions to 共3兲of the form ␺ j ⬘共t兲= ␾ jei⍀⬘t,共4兲 with time independent ␾ jand ⍀⬘. These solutions correspond to stationary dark solutions of the DNLS Eq. 共1兲with frequency ⍀=⍀⬘− ␻ 0. We obtain them numerically 共as an example see Fig. 1兲by continuation of the code ...+1,−1, +1,−1,0,+1,−1,+1,−1... at the anticontinuous limit C =0. As it is well known 关18兴, they show a constant background intensity 共兩 ␺ j兩2→⍀+4C;兩j兩→⬁兲, with zero amplitude at a lattice site. In particular, we are interested in studying the linear stability of these solutions versus the period T=2 ␲ /⍀⬘. Figure 2 shows the moduli of the Floquet eigenvalues for the choice ␻ 0=4. For low enough values of the period 共high enough frequencies, equivalent to anticontinuum limit兲the solutions are stable, but increasing Tabove a critical value T⯝0.48 we observe oscillatory instabilities. These oscillatory instabilities are size dependent and their computation for increasing system sizes suggests that they decay to zero for large enough Tfor any finite chain, but only asymptotically for T→⬁共continuum limit兲in the infinite chain. This result is equivalent to the one observed in Ref. 关18兴for increasing coupling 关with the used parameter values, Fig. 2 appears as a rescaling of Fig. 2 in Ref. 关18兴, with T↔2 ␲ Cand Moduli of Floquet eigenvalues ↔e2 ␲ Re共␭兲兴. FIG. 1. Stationary dark breather for ⍀⬘= ␲ /3, ␻ 0=4 and C=1. B. SÁNCHEZ-REY AND M. JOHANSSON PHYSICAL REVIEW E 71, 036627 共2005兲 036627-2 We remark that, for direct physical applications, the form 共3兲is generally the most relevant form of the DNLS equation, with ␻ 0being the frequency of an uncoupled oscillator in the small-amplitude limit, and ⍀⬘=2 ␲ /Ta physical oscillation frequency. Thus, the frequency ⍀of a solution to 共1兲 normally does not in itself describe a physical oscillation frequency but rather a frequency shift relative to ␻ 0, and therefore ⍀does not need to be sign definite. Evidently, the sign of a frequency also depends on the sign-convention used in 共2兲and 共4兲. III. EXACT DARK WAVES Traveling dark waves were found to appear, together with some radiation, from unstable stationary dark breathers in numerical simulations 关16,18,25,27兴. Typically these moving dark breathers continuously emit radiation and decay 共its minimum intensity increases兲, especially if the dark mode is narrow. However, analogously to the search for moving ’bright’ discrete breathers 关21,22兴, one can search for traveling dark breathers which after minternal oscillations, each of period T, return to a state identical to its initial state but translated ksites. In other words, we can look for solutions to 共3兲that fulfill ␺ j ⬘共mT兲= ␺ j−k ⬘共0兲ei ␣ ⬘,共5兲 where ␣ ⬘is some arbitrary phase. The dark solitary wave numerically found in Ref. 关20兴is a special case of such a solution. The reason for working with Eq. 共3兲instead of Eq. 共1兲is that ␺ j共t兲= ␺ j ⬘共t兲e−i ␻ 0twill be a solution of Eq. 共1兲fulfilling the condition ␺ j共mT兲= ␺ j−k共0兲ei ␣ with a phase ␣ = ␣ ⬘− ␻ 0mT.共6兲 Thus, for fixed parameters C,T,m, and k, we are able to investigate the dependence of an exact mobile solution of the DNLS system on the phase ␣ , just by computing a solution to Eq. 共3兲of the form 共5兲, and doing the continuation versus ␻ 0. In other words, for the DNLS model there are no effects of commensurability, in contrast to the case for general oscillator chains 关21兴. Note that for a stationary dark mode, we must have ␻ 0⬍共2 ␲ /T兲+4Cto have a nonzero background intensity, and that continuation in the direction of larger ␻ 0 共at fixed Tand C兲towards this limit value is equivalent to approaching the continuum dark-soliton limit, i.e., the “dip” in the envelope becomes wider. To obtain numerically an exact traveling dark wave, we have imposed in the Newton iteration the condition 共5兲and, as a guess, we have taken the corresponding stationary solution and applied a sinusoidal perturbation to its imaginary part, which plays a role similar to the velocity in stationary breathers on Klein-Gordon systems. With this method, for m=1 and k=−1, we have found exact dark solutions moving to the left with a phase ␣ ⬘= ␲ , which means that ␺ j−2 ⬘共2T兲= ␺ j−1 ⬘共T兲ei ␲ = ␺ j ⬘共0兲.共7兲 Therefore the envelope velocity of these dark waves is c =2/2T=1/T, since the lattice spacing is 1. In Fig. 3 we show the dynamics of an exact dark wave for the case T=4.5. Note that Eq. 共7兲implies that 兩 ␺ j−1 ⬘共T兲兩2=兩 ␺ j ⬘共0兲兩2. In Fig. 4, we show snapshots of dark waves at times when their envelopes are site centered, which is when their dip intensity is minimal. We observe that the minimum intensity of these moving solutions is small but nonzero. In that sense, they are “gray” modes rather than “black” modes. Another important property is that the background intensity is not constant. These waves show a characteristic oscillating tail which seems to be necessary in order to avoid phonon radiation. We can also appreciate that for fixed Cand ␻ 0, the intensity of the background decreases with the velocity. The same behavior is observed for the minimum intensity. Note that, for a stationary constant-amplitude traveling background wave of the form 共4兲, with ␾ j= ␾ 共0兲eiQj,tobe compatible with the condition 共7兲, its amplitude must fulfill 兩 ␾ 共0兲兩2=Q+共2n+1兲 ␲ T+2C关1 − cos共Q兲兴 − ␻ 0共8兲 for some integer n. For large T, the family of solutions in Fig. 4 approaches the stationary dark mode, corresponding to a constant-amplitude background with n=0 and Q= ␲ . When continuing these solutions towards smaller T, the observed FIG. 2. Oscillatory instabilities of the stationary dark breathers vs the period for different lattice sizes and parameter values C=1, ␻ 0=4. FIG. 3. Time evolution of an exact dark wave. The initial state translates one site each period Tof the internal oscillations. EXACT NUMERICAL SOLUTIONS FOR DARK WAVES ON…PHYSICAL REVIEW E 71, 036627 共2005兲 036627-3 oscillations develop against a background which still, as illustrated in Fig. 4共b兲, approximately can be described by 共8兲 with n=0 and Q= ␲ .共In the smooth continuation, the average phase twist is fixed by using periodic boundary conditions in a finite-size system.兲Thus, these dark modes move through the interaction with the tail envelope oscillations, and not through an average tail phase gradient eiQj,Q ⫽0, ␲ . However, solutions with a phase gradient also exist, and we will return to this point in Sec. IV. One can study the linear stability of these solutions linearizing the map ⑀ j共0兲→ ⑀ j−1共T兲around the exact solution ␺ j ⬘. The linearized equations for ⑀ jare i ⑀ ˙j+2兩 ␺ j ⬘兩2共t兲 ⑀ j+ ␺ j ⬘2共t兲 ⑀ j *+C共 ⑀ j+1 + ⑀ j−1 −2 ⑀ j兲+ ␻ 0 ⑀ j=0. 共9兲 This linearized map defines an extended Floquet matrix with properties analogous to the ordinary Floquet matrix for stationary solutions 共cf. Ref. 关21兴兲. A. Continuation versus T As in the case of stationary dark breathers in finite-size systems, the mobile solutions are stable for high enough values of T共slow enough velocities兲if ␻ 0is not too small. 共In this subsection, we put ␻ 0=4C, so that T→⬁becomes a continuum limit.兲If we perform a numerical continuation of this family of stable solutions towards lower values of T 共higher velocities兲we find bifurcations as Fig. 5 shows. These bifurcations are due to harmonic instabilities 共a pair of eigenvalues collides at +1 and leaves the unit circle along the real axis兲. After emerging the instability the Newton method stops to converge and the continuation path is interrupted. However, we can jump into another family of solutions able to move faster increasing the continuation step. In Fig. 6 we compare a solution just at the bifurcation point T1=8.456 共pluses兲whose continuation stops at 8.385 with a solution found increasing the continuation step up to 8.3 共circles兲. What happens if we continue the latter solution towards the bifurcation point? We observe a fast growth and rearrangements in the oscillatory structure of the tails 共see full circles for T=8.450 in Fig. 6兲until the continuation finally stops at 8.468. The signature of these huge changes in the oscillating tail is an avoided collision at +1 between two conjugated Floquet eigenvalues and the appearance of oscillatory instabilities, as Figs. 7共a兲and 7共b兲show, respectively. The following harmonic bifurcations we found in Fig. 5 occur at T2⯝8.023, T3⯝6.872, T4⯝6.234, T5⯝6.029, T6 ⯝5.034, and T7⯝4.604. For T⬍4.1 more bifurcations appear at T8⯝4.091, T9⯝4.025, T10⯝3.907, and T11⯝3.780 but most of them are not visible because the continuation path is interrupted immediately after a pair of eigenvalues collide at +1. The structure of all these bifurcations is similar to the one described earlier. Figure 8 shows the exact mobile solution at the bifurcation point T5=6.029 共circles兲, which can be continued until T=5.797, and the family of solutions obtained increasing the continuation step up to T=5.6 共squares兲. A three-dimensional plot in Fig. 9 shows the continuation path of this new solution towards the bifurcation point T5. An avoided collision at +1 occurs at 5.6704 关see Fig. 10共a兲兴. After the avoided collision the order of the oscillating tail is destroyed. A projection of Fig. 9 on to the 共j,兩 ␺ j兩2兲plane 共see Fig. 11兲shows that this is due to the FIG. 4. 共a兲Snapshots of exact moving dark breathers corresponding to T=4,T=4.5, and T=6 from top to bottom. Note how the intensity of the background decreases when Tincreases. Other parameters values: C=1, ␻ 0=4. 共b兲Re关 ␺ n共0兲兴 共full circles兲,Im关 ␺ n共0兲兴 共open circles兲, and 兩 ␺ n共0兲兩2共stars兲for an exact moving solution with T=4. FIG. 5. Instabilities found along a numerical continuation towards lower values of T. Shortly after the harmonic bifurcations the continuation path stops. At this point we have to jump into another family of solutions increasing the continuation step. For T⬍7we begin to observe oscillatory instabilities which become stronger. Other parameter values: C=1, ␻ 0=4, N=61. B. SÁNCHEZ-REY AND M. JOHANSSON PHYSICAL REVIEW E 71, 036627 共2005兲 036627-4 activation of an extended mode with wave number close to ␲ . The phenomenon resembles a lot resonances due to higher harmonics observed in Klein-Gordon lattices, e.g., resonances between standing-wave phonons observed in Ref. 关25兴, and the breather-phonon resonance described in Ref. 关31兴along a breather-phonobreather transition. In Fig. 10共b兲 we also show the oscillatory instabilities associated with this process. Small bubbles of oscillatory instabilities are present in the system for T⬍7. They are not visible in Fig. 5 due to the scale on the Yaxis. The size of these bubbles increases as T decreases, and in fact we can appreciate them clearly in Fig. 5 for T⬍5. All the earlier descriptions show that the nature of the oscillating background changes along the continuation path towards smaller T. Clearly, as the velocity increases the number of oscillations increases and, therefore, the number of sites per oscillation decreases. For example, in Fig. 6 we have roughly five oscillations for T=8.3 with about 12 sites per oscillation. In Fig. 8 we can appreciate six oscillations with ten sites per oscillation for T=5.6. For T=4.5 共see Fig. 4兲the number of oscillations increases up to seven with approximately nine sites per oscillation. We again can observe seven oscillations for T=4 but now with eight sites per oscillation 共see Fig. 4 anew兲. We conclude that the harmonic bifurcations along the continuation towards smaller Tare associated with a qualitative change of tail nature, corresponding either to an increase of the number of oscillations, a decrease of the number of sites per oscillation, or both. FIG. 6. Exact mobile dark breathers just at the first bifurcation T1=8.456 共pluses兲and for T=8.3 共circles兲after jumping over the bifurcation. The full circles correspond to the continuation of this last solution towards the bifurcation point. Other parameter values: C=1, ␻ 0=4. FIG. 7. 共a兲When we reverse the continuation path after jumping over a bifurcation, an avoided collision between a pair of conjugated Floquet eigenvalues gives raise to huge changes in the oscillatory tail of the solutions. 共b兲These rearrangements enclose oscillatory instabilities. FIG. 8. Exact mobile dark breathers just at the bifurcation T5 =6.029 共circles兲and for T=5.6 共squares兲after jumping over the bifurcation. Other parameter values: C=1, ␻ 0=4. FIG. 9. Inversion of the continuation path after jumping over the fifth bifurcation. When we approach the bifurcation point the smooth oscillating tails are destroyed. EXACT NUMERICAL SOLUTIONS FOR DARK WAVES ON…PHYSICAL REVIEW E 71, 036627 共2005兲 036627-5 We can obtain an analytical estimate for the dependence of the oscillating tail wave number on the velocity, by considering small-amplitude tail oscillations in the linearized equations. With the ansatz for the tail ␺ j ⬘共t兲=关 ␾ 共0兲+u0ei共qj+ ␻ t兲+v0e−i共qj+ ␻ t兲兴ei共 ␲ j+⍀⬘t兲,共10兲 with ⍀⬘=2 ␲ /Tand ␾ 共0兲=冑⍀⬘− ␻ 0+4C, condition 共7兲is fulfilled if ␻ T−q=2 ␲ n,ninteger. Moreover, inserting the ansatz into 共3兲and linearizing yields the dispersion relation for small-amplitude oscillations around the stationary solution 共cf., e.g., Ref. 关32兴兲, ␻ 2=8Csin2共q/2兲关2Csin2共q/2兲 +共 ␾ 共0兲兲2兴. Combining these two conditions thus determines the possible wave numbers qfor small-amplitude tail oscillations at given T,C, ␻ 0,as q+2 ␲ n±4CT冏sin q 2冏冑sin2q 2+2 ␲ − ␻ 0T 2CT +2=0. 共11兲 The variation of the spatial oscillation period 2 ␲ /qwith T for fixed ␻ 0,C, corresponding to the relevant solutions of 共11兲, is illustrated in Fig. 12. As can be seen, the curve corresponding to a positive sign and n=−1 in 共11兲agrees qualitatively well with the numerical results reported earlier. Note that this solution corresponds to opposite signs of qand ␻ , and thus a wave traveling to the right, i.e., in the opposite direction to the dark breather itself. Thus, the traveling dark breather emits backward radiation behind at the same time as it absorbs radiation in front, which makes possible its existence as an exact solution. B. Continuation versus ␻ 0 For fixed period 共fixed velocity兲we can also continue a given solution versus the phase parameter ␻ 0. We find that the background and the stability of the dark waves are very sensitive to this parameter ␻ 0. When ␻ 0decreases, the dark mode narrows and becomes more discrete, and the background intensity and the oscillation amplitudes increase very quickly. As an example, a mobile dark breather for T=12 and ␻ 0=4 共circles in Fig. 13兲can be continued until ␻ 0=3.64. For that value, we can see in Fig. 13 共full circles兲that the oscillations are very large. If we increase ␻ 0, the dark mode widens and becomes more continuumlike, and we observe the opposite effect. For ␻ 0=4.25 共pluses兲the background intensity has already been reduced to its half, and the backFIG. 10. 共a兲An avoided collision between a pair of conjugated Floquet eigenvalues points out the excitation of the extended mode shown in Fig. 11. 共b兲This excitation gives raise to numerous oscillatory instabilities. FIG. 11. Projection of Fig. 9 on to the 共j,兩 ␺ j兩2兲plane. Notice the excitation of an extended mode with wave number close to ␲ . FIG. 12. Spatial oscillation period 2 ␲ /qvs Tfor smallamplitude tail oscillations of the form 共10兲fulfilling 共7兲, as obtained from 共11兲. Solid 共dashed兲lines correspond to a positive 共negative兲 sign in 共11兲, with, from top to bottom, nranging from −1 to −7 共1 to 6兲. As in previous figures, C=1 and ␻ 0=4. B. SÁNCHEZ-REY AND M. JOHANSSON PHYSICAL REVIEW E 71, 036627 共2005兲 036627-6 ground oscillations are invisible on the scale of the figure, but still present 共they are of order 10−6兲. This continuation path ends, for the considered system with N=61, at ␻ 0 =4.516 共close to the limit value ␻ 0=2 ␲ /T+4Cfor the stationary mode in the infinite system兲with a background intensity around 0.007. This result is very interesting for various reasons. On the one hand, varying ␻ 0one can find stable mobile solutions with huge oscillating tails if the velocity is slow enough. An example is the solution shown in Fig. 13 with full circles. On the other hand, for a given unstable mobile solution one can adjust the phase ␻ 0in order to stabilize it. In Fig. 14 we show how we can avoid the harmonic instability at T =8.4⬍T1decreasing ␻ 0. Another alternative would be to jump into another stable family of dark mobile solutions increasing the continuation step from ␻ 0=4.013 共where the continuation path stops兲to ␻ 0=4.02. For higher velocities the scenario becomes more complex due to the appearance of more and stronger instabilities. But again we can find windows of stability managing ␻ 0as Fig. 15 shows for the case T=5. The bifurcations observed in the continuations versus ␻ 0 at fixed Tare essentially the same types of bifurcations as those discussed in Sec. III A in the continuation versus Tat fixed ␻ 0. Generally, since we have two-parametric families of solutions, each bifurcation should define a bifurcation curve in the 共 ␻ 0,T兲plane 关cf., e.g., a similar scenario for second-harmonic resonances of standing waves in KleinGordon chains, illustrated in Ref. 关25兴, Fig. 16共a兲兴. Performing the continuations at fixed Tand fixed ␻ 0, respectively, we will generally hit these curves from different directions. Note also that expression 共11兲for the linear resonances depends not only on Tbut also on ␻ 0, since the background amplitude depends on ␻ 0. Thus, the linear spatial oscillation period changes also when varying ␻ 0at fixed T, so that a similar set of resonances should be expected 共and is also numerically observed兲also in this continuation. 共Another way of expressing this is, that the bifurcation curves in the 共 ␻ 0,T兲plane generally are not horizontal or vertical lines.兲 IV. STATIONARYAND MOVING GRAY SOLITONS WITH PHASE GRADIENT From a “quasicontinuum” NLS approximation of the DNLS Eq. 共3兲, taking into account the discrete dispersion for the background wave ␾ 共0兲ei共Qj+⍀⬘t兲but neglecting other effects of discreteness, one finds 关16兴that the envelope velocity cand the minimum intensity 兩 ␺ min兩2for wide, continuumlike discrete dark solitons should be related by c=2Csin共Q兲±兩 ␺ min兩.共12兲 Thus, the minimum intensity determines the envelope velocity relative to the group velocity of the background wave. It also determines the total phase shift across such a solution 共relative to the background wave兲as ␦ =2 arccos兩 ␺ min/ ␾ 共0兲兩. In particular, the quasicontinuum NLS approximation predicts the existence of a stationary 共c=0兲gray soliton with 兩 ␺ min兩=2Csin共Q兲. The condition 兩 ␺ min兩⬍兩 ␾ 共0兲兩would thus limit the existence of such solutions to 兩2Csin共Q兲/ ␾ 共0兲兩⬍1. However, as was discussed in Ref. 关16兴, the quasicontinuum NLS approximation does in general not accurately describe wide small-amplitude 关 ␺ min close to ␾ 共0兲兴gray solitons in the discrete system, since it does not take into acFIG. 13. Exact dark waves for T=12 and three different values of ␻ 0: 4.0 共circles兲, 3.64 共full circles兲, and 4.25 共pluses兲. C=1. FIG. 14. Instabilities vs the phase ␻ 0for a fixed value of the period T=8.4. FIG. 15. Instabilities vs the phase ␻ 0for a fixed value of the period T=5. EXACT NUMERICAL SOLUTIONS FOR DARK WAVES ON…PHYSICAL REVIEW E 71, 036627 共2005兲 036627-7 count effects of higher-order dispersion. Instead, by taking into account also third and fourth-order dispersion in the continuum approximation and performing a multiscale expansion, it is possible to derive a Korteweg-de Vries 共KdV兲 equation for the soliton amplitude in the limit of small amplitude 共see also Ref. 关33兴兲. The envelope of the corresponding family of soliton solutions, parametrized by the small parameter ␮ , can be written as 关16兴: 兩 ␺ j共t兲兩 = ␾ 共0兲−12 ␮ 2G2 G1sech2关 ␮ 共n−共Vg+c0+W兲t兲兴, 共13兲 where Vg=2Csin共Q兲,c0 2=−2关 ␾ 共0兲兴2Ccos共Q兲,G1 =4 ␾ 共0兲Ccos共Q兲关−3+Vgc0/4C2cos2共Q兲兴,G2=C2cos2共Q兲关1 +共 ␾ 共0兲兲2/6Ccos共Q兲兴−Vgc0/3, and W=−2 ␮ 2G2/c0.共In Ref. 关16兴only the case Vg=0 was considered, leading to a vanishing third-order dispersion. However, the third-order dispersion present when Vg⫽0 only leads to a renormalization of the coefficients in the KdV equation; see also Ref. 关34兴.兲 Thus, the condition for a stationary small-amplitude gray soliton becomes Vg+c0+W=0. The limit of vanishing soliton amplitude ␮ 2G2→0 then gives the limit for existence of a stationary gray solution as Vg=−c0, i.e., it may exist when 冏2Csin2共Q兲 关 ␾ 共0兲兴2cos共Q兲冏⬍1, 共14兲 which notably agrees with the expression obtained from the quasicontinuum NLS limit only when 兩tan共Q兲兩=1.共Note that, as was discussed in Ref. 关16兴, the family 共13兲of solutions also exists for G2/G1⬍0, leading to “antidark” pulses.兲 To search numerically for stationary gray solutions as exact solutions to the DNLS Eq. 共3兲, we impose the boundary conditions ␾ N+1=eiQ ␾ N, ␾ 0=e−iQ ␾ 1and perform a continuation versus Q, using the dark 共black兲solutions with Q= ␲ and ␾ 共0兲=1 previously obtained 关18兴for C⬎0 as trial solutions in the Newton scheme, with an additional phase torsion eiQj imposed on the background. 共Our method is analogous to that used in Ref. 关35兴for Klein-Gordon lattices.兲To keep the background intensity to unity, we vary the frequency ⍀ as ⍀=1+2C关cos共Q兲−1兴. An example with Q=0.98 ␲ is shown in Fig. 16. This solution can be continued towards larger C, and generally becomes smoother and “grayer” 共i.e., 兩 ␺ min兩increases兲as Cis increased for fixed Q. The scenario with oscillatory 共Krein兲instabilities is similar to the scenario for the black solitons 关18兴as long as Qis close to ␲ . However, some qualitative differences in the continuation scenario are observed: 共i兲For each Q, ␲ /2⬍Q⬍ ␲ , there is a maximum value of Cfor existence of a stationary gray mode. To numerical accuracy, this maximum value is given by Cmax共Q兲 =兩cos共Q兲兩/2 sin2共Q兲as predicted by 共14兲, and corresponds to the point where 兩 ␺ min兩= ␾ 共0兲, and the smooth continuumlike gray soliton bifurcates with the plane wave with wave number Q. Note that Cmax is a monotonously increasing function of Qin the interval 关 ␲ /2, ␲ 关, with Cmax共 ␲ /2兲=0 and Cmax共 ␲ 兲→+⬁. FIG. 16. Stationary discrete gray breather in DNLS, corresponding to a background wave with Q=0.98 ␲ .共a兲shows, for C =0.075, intensity 兩 ␺ j兩2共minimum intensity is 兩 ␺ 31兩2⬇0.0052兲;共b兲 shows nearest-neighbor phase shifts ␣ j+1− ␣ j共with notation ␺ j =兩 ␺ j兩ei ␣ j兲共note that the total phase shift across the central site ␣ 32 − ␣ 30 is not ␲ as it would be for a black breather兲. Lower figures with stability eigenvalues show that the solution is close to the threshold of Krein instability: C =0.075 is stable 共c兲while C =0.076 is unstable 共d兲. B. SÁNCHEZ-REY AND M. JOHANSSON PHYSICAL REVIEW E 71, 036627 共2005兲 036627-8 共ii兲Attempting a continuation towards the anticontinuous limit C=0 at fixed Q, there are two possible qualitatively different scenarios depending on Q. For ␲ /2⬍Qⱗ0.92 ␲ the continuation is monotonic, and ends at the constantamplitude phase-twisted solution 共...,e−2iQ,e−iQ,+1,eiQ, −1,e−iQ,+1,eiQ,e2iQ,...兲,atC=0, where the “−1” denotes the central site. The total phase twist with respect to the background wave at the anticontinuous limit is thus ␦ =2共2Q− ␲ 兲. The scenario is illustrated in Fig. 17 for Q =0.6 ␲ . Thus, for each Qthere is an “optimal” value of C where the stationary gray soliton contrast is largest. Note also that the instability scenario is different than for Q= ␲ : the unstable complex eigenvalues behave similarly close to the continuum limit, but when decreasing Cthey do not return to the imaginary axis but collide with each other on the real axis, so that close to the anticontinuous limit there are two unstable real eigenvalues. Thus, strictly speaking, the gray solitons are always unstable for infinite systems for these Q. For 0.92 ␲ ⱗQⱗ ␲ , a monotonic continuation of the continuumlike gray soliton to the anticontinuous limit C=0 at fixed Qis not possible due to a bifurcation, associated with a collision of eigenvalues at ␭=0. The location of this bifurcation line in the Q-Cplane is illustrated in Fig. 18共a兲. Similarly, the continuation of the anticontinuous solution discussed earlier to larger Cis now interrupted by another bifurcation close to C=0.14 共depends only weakly on Q兲. Thus, there is a region in the Q-Cplane where these two solutions exist simultaneously, but have different properties 关see Fig. 18共b兲兴. The solution continued from the anticontinuous limit is always unstable, while the solution continued from the continuous limit is stable only in the area in the lower right corner in Fig. 18共a兲, above which the oscillatory instability sets in. Thus, stable stationary discrete gray solitons only exist 共for the infinite chain兲for 0艋Cⱗ0.076 and 0.975ⱗQ/ ␲ 艋1. In a similar way as described in Sec. III, it is possible to find also exact moving gray discrete solitons with an average phase torsion Q⫽ ␲ , by using perturbations of the stationary FIG. 17. 共a兲Continuation of stationary discrete gray soliton from continuous limit C=Cmax 关from 共14兲兴 to anticontinuous limit C=0 at fixed Q=0.6 ␲ . The smallest minimum intensity 兩 ␺ min兩2 ⯝0.9398 is obtained for C⯝0.09. Only the central part of a larger chain is shown. 共b兲shows the real part of the stability eigenvalues during the continuation. The unstable eigenvalues are complex for Cⲏ0.06 and real for Cⱗ0.06. FIG. 18. 共a兲Solid line: Location of the bifurcation interrupting a monotonic continuation at fixed Qof a continuumlike gray soliton towards the anticontinuous limit C=0. Dashed line: Location of the Krein collision causing oscillatory instability. A stable stationary discrete gray soliton exists only in the area between the lines 共lower right part of the figure兲.共b兲Patterns for the central parts of two simultaneously existing solutions at Q=0.95 ␲ ,C=0.13. 共⫹兲: Solution continued from continuous limit; 共⫻兲: solution continued from anticontinuous limit. EXACT NUMERICAL SOLUTIONS FOR DARK WAVES ON…PHYSICAL REVIEW E 71, 036627 共2005兲 036627-9