Demonstration of the stability or instability of multibreathers at low coupling

Demonstration of the stability or instability of multibreathers at low coupling JFR Archilla 1, J Cuevas, B S´anchez–Rey and A Alvarez* Nonlinear Physics Group of the University of Sevilla Dep. F´ısica Aplicada I, ETSI Inform´atica and *Facultad de F´ısica, p. 5 Avda Reina Mercedes s/n, 41012 Sevilla, Spain Abstract Whereas there exists a mathematical proof for one–site breathers stability, and an unpublished one for two–sites breathers, the methods for determining the stability properties of multibreathers rely in numerical computation of the Floquet multipliers or in the weak nonlinearity approximation leading to discrete non–linear Schr¨odinger equations. Here we present a set of multibreather stability theorems (MST) that provides with a simple method to determine multibreathers stability in Klein–Gordon systems. These theorems are based in the application of degenerate perturbation theory to Aubry’s band theory. We illustrate them with several examples. Key words: Discrete breathers; Multibreathers; Intrinsic localized modes; PACS: 63.20.Pw, 63.20.Ry, 66.90.+r. 1 Introduction Discrete breathers are time–periodic, localized oscillations in discrete systems due to a combination of nonlinearity and discreteness. They have become a well understood phenomenon since the publication of the proof of existence in Ref. [1]. This proof is the origin of exact and powerful numerical methods to calculate them and to determine their stability [2]. A deeper insight has been achieved since the introduction of Aubry’s band theory [3]. The last reference, together with Ref. [4] can be considered as reviews, although a new one is badly needed after the huge development of the subject in the last years. 1Corresponding author. E-mail: [email protected] Preprint submitted to Physica D 9 February 2003 The stability of one–site breathers is proofed in Refs. [2,3] under rather general conditions, the stability of the possible two–sites breathers in analyzed in an unpublished theorem in Ref. [5], page 69. Other methods use the hypothesis of small amplitude oscillations or the rotating wave approximation [6,7] leading to the Nonlinear Schr¨odinger equation (DNLS), which is the actual equation analyzed, and, there are, of course, the efficient but slow numerical methods mentioned above. Here we propose a method based in the properties of the bands in Aubry’s bands theory and perturbation theory, to obtain the stability properties of any multibreather at low coupling. In many cases it involves only some modification of the coupling matrix and the knowledge of the hardness/softness of the on–site potential. We will briefly summarize the band theory in the next section, while we introduce the notation and some basic concepts, on which our theory is based. However, we refer to Ref. [3,1,5] for detailed explanations, as it will be long and repetitive to expose it in detail. In section 3 we develop the method for symmetric on–site potentials or in– phase multibreathers, which is synthesized in a theorem, commented in section 4. Its scope is enlarged to non–symmetric potentials in section 3, and to generalized Klein–Gordon systems in section 6. The method is applied to several interesting examples in sections 7 and 8. In Appendix A we calculate the value of a magnitude γused in our theory, and in Appendix B we relate the curvature of the bands with the characteristics of the on site potential. 2 Band theory and notation 2.1 The Newton operator We consider Klein–Gordon systems with linear coupling described by dynamical equations of the form: ¨un+V′(un) + ε N ∑ m=1 Cn mun= 0 n= 1, . . . , N (1) where the variables unare functions of time t,V(un) is an homogeneous on– site potential, V′its derivative, ˙undenotes derivation with respect to time, Nis the number of oscillators, Cis a coupling constant matrix, which can describe nearest neighbor or long–range interaction, and includes the boundary conditions, and εis the coupling parameter. We use a notation similar to quantum mechanics, i.e., |u⟩ ≡ [u1(t), . . . , uN(t)]†(†meaning the transpose matrix). Defining V(u) = [V(u1), . . . , V (uN)]†and analogously its derivatives, 2 equation (1) can be written as: |¨u⟩+|V′(u)⟩+ε C |u⟩= 0 .(2) We suppose that the functions un(t) are time–periodic (with period Tand frequency ωb), time–reversible solutions and, therefore, they can be written as cosine Fourier series with real coefficients (some of these assumptions will be relaxed later). The (linear) stability properties of a given solution |u⟩depend on the properties of the characteristic equation for the Newton operator Nε given by Nε(u)|ξ⟩≡|¨ ξ⟩+V′′(u)·|ξ⟩+ε C |ξ⟩=E|ξ⟩,(3) where ·product is the list product, i.e., f(u)·|ξ⟩is the column matrix with elements f(un(t)) ξn(t). If E= 0, this equation describes the evolution of small perturbations |ξ⟩of |u⟩. 2.2 The Floquet matrix Any solutions of Eq. (3) can be determined by the column matrix of the initial conditions for positions and momenta Ω(0) = [ξ1(0), . . . ξN(0), π1(0), . . . πN(0)]†, with πn(t) = ˙ ξn(t), and Ω(t) describes its evolution in the space of coordinates and momenta. A base of solutions is given by the 2 Nfunctions with initial conditions Ων(0), ν= 1, . . . , 2N, with Ων l(0) = δν l. We will identify often a given solution ξ(t) of the Newton equation with the corresponding matrix of initial conditions Ω(0). As the Newton operator depends on the T–periodic solution u(t), it is also periodic, and the evolution of the solutions of Eq. (3) can be studied by means of the Floquet operator, which maps ξ(t) into ξ(t+T). In a finite system this is equivalent to the Floquet matrix FEbeing given by:   {ξn(T)} {πn(T)}  =FE  {ξn(0)} {πn(0)}  (4) These matrices can be easily constructed numerically, by integrating Eq. (3) 2Ntimes from t= 0 to t=T, with initial conditions Ων(0). Then, the ν column of FEis given by Ων(T). 2.3 Stability and bifurcations The eigenvalues {λl}2N l=1 of F0, called the Floquet multipliers, determine the linear stability of the solution u(t). If there is any eigenvalue with |λl|>1, the corresponding eigenfunction ξl(t) grows with time and u(t) is unstable; if λl= 1, ξl(t) is T-periodic; if λl=−1, ξl(t) is 2 T–periodic. However, F0(and 3 any FE) is a symplectic matrix, as it is derived from the symplectic system Eq. (3), which implies that if λlis a non–zero eigenvalue, so it is 1/λ.FEis also real, which in turn implies that λ∗ land 1/λ∗are also multipliers. Therefore, the eigenvalues of FEcome in groups of four, if they are not real, or in pairs, if they are real or if |λl|= 1. We can write the Floquet multipliers as {exp(iθl)}, with θl, in general, complex numbers. The complex number iθlare called the Floquet exponents, and θlthe Floquet arguments. The only possibility for the stability of u(t) is that all the eigenvalues of F0 have moduli 1, i.e., they are at the unit circle. Therefore, the condition for linear stability of u(t) can be described as all the Floquet arguments for F0 being real. Moreover, if u(t) is stable and, then, every |λl|= 1, the Floquet multipliers come in complex conjugate pairs with arguments (±θl) or as double 1 or −1 (θ= 0 or θ=±π). If a parameter like the coupling εis changed the multipliers of F0change continuously. Therefore, a bifurcation to an unstable solution can only take place in three different forms: a) Harmonic instability: two complex eigenvalues moving along the unit circle collide at λ= 1 (θ= 0); b) Subharmonic instability: two complex eigenvalues moving along the unit circle collide at λ=−1 (θ=±π); c) Oscillatory or Hopf instability: two pairs of complex eigenvalues moving along the unit circle collide at ±θ= 0 and abandon the unit circle as a quadruplet (λ, 1/λ, λ∗,1/λ∗). A Floquet multiplier of F0is always known. Calculating the derivative of Eq. (2) with respect to time we obtain that N˙u= 0. ˙uis also T–periodic, therefore, it is an eigenfunction of F0with multiplier 1, and as they come in pairs, there is always a double 1 multiplier. ˙uis called the phase mode because its meaning is that if ˙u(t) is a solution of Eq. (2), ˙u(t+dt) is also a solution. While the solution u(t) exists, this double eigenvalue is always there. Due to the possible forms of the bifurcations, if all the multipliers are isolated except the double 1, the system is structurally stable, i.e., there is a neighborhood in the space of the parameters where u(t) is stable. 2.4 Aubry’s band theory However, it turns out that much information can be obtained by studying also the Floquet arguments for E= 0, with is known as Aubry’s band theory [3]. The set of points (θ, E), with θa real Floquet argument of FEhave a band structure. The fact that the Floquet multipliers come in pairs of complex conjugate pairs brings about that if (θ, E) belongs to a band, (−θ, E) does it too, i.e., the bands are symmetric with respect to θand dE/ dθ(0) = 0. As ˙u has eigenvalue 1 (θ= 0), there is always a band tangent to the axis E= 0 at θ= 0. There are at most 2 Npoints for a given value of Eand, therefore, there are at most 2 Nbands crossing any horizontal axes in the space of coordinates 4 (θ, E). The condition for linear stability of u(t) is equivalent to the existence of 2 Nbands crossing the axis E= 0 (including tangent points with their multiplicity). If a parameter like the coupling εchanges, the bands evolve continuously, and they can loss crossing points with E= 0, bringing about the instability. We know more about the eigenfunctions of the Newton operator. The fact that it is periodic, or, in other words, that it commutes with the operator Pof translation in time a period, Pf(t) = f(t+T), implies that its eigenfunctions can be chosen simultaneously as eigenfunctions of P, which by Bloch theorem, are given by ξ(θ, t) = χ(θ, t) exp(i θ t/T), χ(θ, t) being a column matrix of T–periodic functions. It is straightforward to check that exp(iθ) is the corresponding Floquet multiplier, and θits Floquet argument (FEis simply the representation of Pin the base Ωl). 2.5 Bands at the anticontinuous limit The key concept on the demonstration of breather existence [1], single breather stability [3] and the present paper is the anticontinuous limit, i.e., the system with all the oscillators uncoupled, ε= 0, in Eq. (1–2). At the anticontinuous limit, Eq. (1) reduces to Nidentical equations: ¨un+V′(un) = 0 .(5) Supposing that we consider time–reversible solutions around a single minimum of V, there are only three different solutions: a) oscillators at rest un= 0; b) excited oscillators with identical un(t), hereafter denoted u0(t); c) excited oscillators with a phase difference of πwith the previous ones, given by un(t) = u0(t+T/2). Each site index is given a code σn, which takes elements in {0,1,−1}, where σn= 0 represents an oscillator at rest, un= 0; σn= 1, an oscillator with solution u0(t); and σn=−1, the solution u0(t+T/2). The matrix of codes σ= [σ1, . . . , σN]†represents the state of the system at the anticontinuous limit. We suppose that there are poscillators at rest and N−pexcited oscillators. Equation (3) at ε= 0 becomes N0(u)|ξ⟩≡|¨ ξ⟩+V′′(u)·|ξ⟩=E|ξ⟩,(6) or, equivalently, N−pidentical equations: ¨ ξn+V′′(un)ξn=E ξn,(7) with only a periodic eigenfunction for E= 0, the phase mode ˙un(t) for the isolated, excited oscillators. The other one is the growth mode, (see Appendix B) 5 which is not bounded and will not be used in this work. They give rise to N−p bands tangent to the axis E= 0, shown in Fig. 1. They have positive curvature at (θ, E) = (0,0) if the on–site potential is soft ( dH/ dωb<0) and negative if it is hard ( dH/ dωb>0), as it is demonstrated in Appendix B. The possibility of dH/ dωb= 0 is excluded by the conditions of the breather existence theorem [1], which means that the on-site potential is truly nonlinear for the isolated oscillators with the breather frequency ωb. If, as we are supposing here, the on–site potentials V(un) are identical, the solutions un(t) are also identical except for a change of phase and the bands are superposed, if not, we might have different bands but with the same general shape. The remaining pequations corresponding to the oscillators at rest are of the form: ¨ ξn+ (ω0)2ξn=E ξn,(8) with ω0=√V′′(0). They have only a T-periodic solution for E= 0 , the null solution, due to the non–resonance condition for breather existence nωb=ω0 (p∈N) [1], i.e., none of the harmonics of the breather resonates with the oscillators at rest. They provide pidentical bands which are easily calculated solving the equation above. Its solutions are ξ= exp(±i√ω2 0−E t), with Floquet multipliers exp(±i√ω2 0−E T) and Floquet arguments θ=±√ω2 0−E T. That is, the bands are given by E=ω2 0−ω2 b(θ/T)2, where θcan be reduced to the first Brillouin zone [−π, π] by the addition of 2πp,p∈Z. If the oscillators have different rest frequencies, the rest bands are not superposed but they have the same characteristics. The rest bands are also shown in Fig. 1. Note that this figure is only a sketch, for clarity, as very often some excited oscillators bands are very flat and difficult to appreciate at the same scale. This sketch reproduces, however, the basics facts of the band structure. When the coupling is switched on, the degeneracy of the bands is generically raised. All the tangent bands at (θ, E) = (0,0) except one, than continues there, move upwards or downwards. If the on–site potential is soft, Fig. 1 (left), and a band moves upwards, a double tangent point with the axis E= 0 is lost, or, in other terms, a pair of Floquet arguments of F0becomes complex and the solution u(t) is unstable. If the on site potential is hard, Fig. 1 (right), the same occurs when a band moves downwards. 3 Multibreathers stability In this work use degenerate perturbation theory [8] to demonstrate the stability or instability of the breathers of any code. Degenerate perturbation theory 6 establishes that if N0is a linear operator with a degenerate eigenvalue E0, with eigenvectors {|vn⟩}, which are ortonormal with respect to a scalar product, i.e., ⟨vn|vm⟩=δn m, and if ε˜ Nis a perturbation of N0, with εsmall; then, to first order in ε, the eigenvalues of N0+ε˜ Nare E0+ε λi, with λibeing the eigenvalues of the perturbation matrix Qwith elements Qn m =⟨vn|˜ N|vm⟩. Note that perturbation theory as described in the reference cited, is time independent perturbation theory, but our time variable tis their spatial coordinate x. Let us consider again Eq. (3) with zero coupling, ε= 0. As explained in the previous section, if there are N−pexcited oscillators, there N−pzero eigenvalues corresponding to the T–periodic eigenfunction ˙un(the phase modes of the isolated oscillators), or in other words, a N−ptimes degenerate eigenvalue E0= 0 if we restrict the domain of Nto periodic functions. What we need to know is the sign of this degenerate eigenvalues when the coupling is switched on. To apply degenerate theory to multibreather stability we need to identify the perturbation operator, a suitable scalar product and a ortonormal basis of the eigenspace with eigenvalue E0= 0. The scalar product is defined as: ⟨ξ1|ξ2⟩= N ∑ n=1 ∫T/2 −T/2 ξ∗ 1(t)ξ2(t).dt . (9) Let us suppose initially that all the excited oscillator are identical and vibrate in phase. We will denote by u0(t) these identical solutions. In this case, all the phase modes are also identical and will be denoted as ˙u0. The N−pelements of the basis are |n⟩=1 µ[0, . . . , 0,˙u0,0,...,0]†,(10) with the non-zero element at the position n,nbeing the index of the excited oscillators, and µ=√∫T/2 −T/2( ˙u0)2dt. It is straightforward to check that they are ortonormal. It is enough to know the eigenvalues corresponding to periodic and real solutions, i.e., with Floquet argument θ= 0, as the intersections of the bands with the axis θ= 0 correspond to periodic solutions of (3). Therefore the N−p|n⟩’s form the basis of the degenerate eigenvalue E0= 0 needed to apply perturbation theory. To obtain the perturbation operator, we expand in Taylor series Eq. (3) at ε= 0 and obtain to first order in ε: Nε(u)|ξ⟩=|¨ ξ⟩+V′′(u)·|ξ⟩+ε(V′′′(u)·uε·|ξ⟩+C|ξ⟩)= = (E0+ε λi)|ξ⟩,(11) 7 −2 0 2 −1 −0.5 0 0.5 1 1.5 θ E −2 0 2 −1 −0.5 0 0.5 1 1.5 θ Fig. 1. Sketch of the band structure for a soft on–site potential (left) and a hard on–site potential (right). The continuous lines correspond to excited oscillators and the dashed ones to oscillators at rest. with uε=(∂u ∂ε )ε=0 and uis also the ε= 0 solution. Therefore the perturbation operator of N0is ˜ N |ξ⟩=V′′′(u)) ·uε·|ξ⟩+C|ξ⟩,(12) where uand uεare calculated at ε= 0. We do not know uε, but deriving with respect to ε, at ε= 0, the dynamical equations (2) we obtain: |¨uε⟩+V′′(u)·|uε⟩+C|u⟩= 0 or N0|uε⟩=−C|u⟩.(13) The perturbed eigenvalues Eiare ε λi,λibeing the eigenvalues of the perturbation matrix Q, with elements ⟨n|˜ N|m⟩. The matrix ˜ Cof elements ⟨n|C|m⟩ is simply the matrix Cwithout the columns and rows corresponding to the oscillators at rest. The other terms are: ⟨n|V′′′(u)·uε|m⟩= 1 µ2∫T/2 −T/2 [. . . , 0,˙u0 n,0, . . . ] [. . . , 0, V ′′′(u0)um,ε ˙u0 m,0, . . . ]†dt= δn m µ2∫T/2 −T/2 ˙u0V′′′(u0)un,ε ˙u0dt , (14) with un,ε =(∂un ∂ε )ε=0. Thus, only the diagonal elements ⟨n|V′′′(u)·uε· |n⟩ 8 are non–zero. To calculate the last integral in (14) we will integrate by parts and use that the integral in a period of the derivative of a periodic function is zero. Besides, the functions un,ε are periodic as the coefficients of their Fourier series are given by the derivatives with respect to εof the Fourier coefficients of un. In the deduction below, all the integral limits are −T/2 and T/2, and the terms between brackets from integration by parts will be zero. The last integral in Eq. (14) becomes: [˙u0un,ε V′′(u0)]T/2 −T/2−∫V′′(u0) ˙u0˙un,ε dt−∫V′′(u0) ¨u0un,ε dt= −[V′(u0) ˙un,ε,]T/2 −T/2+∫V′(u0) ¨uε,n,dt−∫V′′(u0) ¨u0un,ε dt= −∫¨u0(¨un,ε +V′′(u0)un,ε) dt(15) The term between parentheses, is the ncomponent of the lhs of Eq. (13), i.e., it becomes −∑mCn m u0 m, where u0 m=u0, if the oscillator mis excited, and zero otherwise, i.e., it is −∑m˜ Cn m u0=−(∑m˜ Cn m)u0. Equation (15) becomes: (∑ m ˜ Cn m)∫¨u0u0dt = (∑ m ˜ Cn m)([˙u0u0]T/2 −T/2−∫( ˙u0)2dt)=−(∑ m ˜ Cn m)µ2.(16) That is, equation (14), leads to: ⟨n|V′′′(u)·uε·|n⟩=−∑ m ˜ Cn m at ε= 0 .(17) Therefore the diagonal elements of the perturbation matrix Qare Qn n =−∑ m ˜ Cn m +˜ Cn n =−∑ ∀m=n ˜ Cn m =∑ ∀m=n Qn m .(18) To summarize, the perturbation matrix Qis given by: Qn m =˜ Cn m, n =m , Qn n =−∑ ∀m=n Qn m ,(19) ˜ Cbeing the coupling matrix without the prows and columns corresponding to oscillators at rest. This result can be very easily extended, to the case when there are oscillators out of phase, at least for an even potential V(u). In this case u0(t+T/2) = −u0(t) and ˙u0(t+T/2) = −˙u0(t). Let us define a new code ˜σ, which is equal to σ, but with +1 instead of zeros, and perform the ansatz un= ˜σn˜un. By 9 2δn m −δ|n−m|,1and CN1=C1N=−1, that is: C=           2−1 0 . . . 0−1 −1 2 −1 0 . . . 0 . . . . . . . . . . . . . . . . . . −1 0 . . . 0−1 2           .(34) If the two oscillator are in phase, the perturbation matrix is Q=   1−1 −1 1   , with eigenvalues (λ1, λ2) = (0,2), the first one corresponds to the phase mode, the second to the antisymmetric eigenmode. Thus, one of the two bands tangent to the axis E= 0 at zero coupling will remain there, and the other will move upwards. If the on–site potential is soft, i.e., with positive curvature, a intersection/tangent point will be lost, and the two site breather is unstable. For the hard on–site potential, the curvature of the bands is negative, and when the band moved upwards, the intersection points will be kept, and the system is be stable. If ε < 0 the signs of the eigenvalues ε Eiare reversed and so are the conclusions. For the out–of–phase, two-site breather, with code σ= [1,−1], the perturbation matrix is: Q=  −γ γ γ−γ  , with eigenvalues (λ1, λ2) = (0,−2γ). The conclusions above are reversed, being the instability mode, when appropriate, the symmetric one. Note that the boundary conditions in Chave no effect. Theorem 4 (Two–site breathers) Consider time–reversible, two–site breathers in an homogeneous Klein-Gordon system with attractive linear coupling, then: a) The two-site breathers with codes ±[1,1] and soft on–site potential are unstable and with ±[−1,1] are stable. b) The two-site breathers with code ±[1,1] and hard on–site potential are stable and with ±[−1,1] are unstable. c) If the coupling is repulsive, the conclusions are reversed. 16 7.2 Three–site breathers Let us renumber the sites 1,2,3 and suppose that the code of the site in the middle is σ2,= +1 as a reference, i.e., σ= [σ1,1, σ3]. We define αi= 1 if σi= 1 and αi=−γif σi=−1, for i= 1,2. Using Cfrom Eq. (34) the perturbation matrix is: Q=       α1−α10 −α1α1+α3−α3 0−α3α3        .(35) Its nonzero eigenvalues are λ±=α1+α3±√(α1+α3)2−3α1α3. In order to have both eigenvalues with the same sign, we need that α1α3>0, i.e., either α3=α1= 1, or α3=α1=−γ. That is, if σ1=σ3= 1, λ±= 1,3>0, if σ1=s3=−γ,λ±=−γ, −3γ < 0, and if σ1σ2=−1, then λ±have different signs. Therefore: Theorem 5 (Three–site breathers) Consider time–reversible, three–site breathers in an homogeneous K-G system with nearest–neighbor, linear coupling, attractive for ε > 0. Then: a) The three–site breathers with codes σ=±[−1,1,−1] are stable if ε > 0 (ε < 0)and the on–site potential is soft (hard). b) The three–site breathers with code σ=±[1,1,1] are stable if ε > 0 (ε < 0) and the on–site potential is hard (soft). c) All other three-site breathers are unstable for any sign of εand type of on–site potential. We can compare the predictions of the values given by the theory with the numerically calculated values of the intersection points of the bands with the axis θ= 0. As an example, the numerically calculated eigenvalues and the predictions given by theorems 1 and 2 for the three–site breather with code [−1,1,−1] are plotted in Figure 3. The very good accuracy is evident. 7.3 A model with long–range interaction An example of system with long range interaction is the twist model used in Ref. [9]. The dynamical equations are: ¨un+V′(un) + ε n+N/2 ∑ m=n−N/2 cos[θtw(n−m)] |n−m|3um= 0 (36) 17 0.01 0.02 0.03 0.04 0.05 0.06 0.07 −0.25 −0.2 −0.15 −0.1 −0.05 0 ε E Fig. 3. Comparison between the eigenvalues Eof the Newton operator obtained numerically (points), predicted by the symmetric MST, theorem 1, (dashed lines) and by the non–symmetric MST, theorem 2, (continuous line), which is the right one to be used. Morse potential, code [−1,1,−1], frequency ωb=γ= 0.8. Vbeing the Morse potential, and θtw the angle between two dipole moments corresponding to neighboring base pairs in a simplified model of DNA. The coupling matrix elements are: Cn m =cos[θtw(n−m)] |n−m|3,(37) Some of the results, also checked numerically and with good accuracy with respect to the eigenvalues E, given by theorem 1, in spite of Vbeing non– symmetric are: Code θtw = 0 θtw =π 11 Stable Unstable 1 -1 Unstable Stable 101 Stable Stable 111 Stable Unstable 1-1 1 Unstable Stable 11-1 Unstable Stable As an example of the effect of the symmetry coefficient γin theorem 2, we can consider the code [1,−1,1] and plot the Newton eigenvalues Eas a function 18 0 0.5 1 1.5 2 2.5 3 −4 −3 −2 −1 0 1 2 3 4 θtw λ Fig. 4. Eigenvalues λof the perturbation matrix Qfor a three-site breather with code [−1,1,−1] on the twist model. Continuous line: γ=ωb= 1; dashed line: γ=ωb= 0.8; and dotted line: γ=ωb= 0.1. Only values of γvery far from 1 lead to a different stability prediction. of θtw. The perturbation matrix is: Q=       γa −b−γa b −γa 2γa −γa b−γa γa −b        ,(38) with a= cos(θtw) and b= cos(2θtw)/8. Figure 4 shows the dependence of the eigenvalues λon the twist angle θtw for three values of the symmetry coefficient γ. It can be seen that only values of γvery far from 1 lead to a change on the stability prediction. 8 Multibreathers, phonobreathers and dark breathers. Parity instabilities. For larger, time–reversible multibreathers in homogeneous Klein–Gordon systems with nearest–neighbor, linear coupling, it is trivial to obtain the stability properties of the two different solutions with wave number q= 0, all the oscillators in phase, or q=π, neighboring oscillators out of phase. Consider any number N′of contiguous oscillators in phase in a system coupled by the 19 matrix Cin Eq. (34). The N′×N′perturbation matrix is given by: Q=               1−1 0 . . . 0 0 −1 2 −1 0 . . . 0 . . . . . . . . . . . . . . . . . . 0 0 . . . −1 2 −1 0 0 . . . 0−1 1               .(39) The characteristic equation for Q x =λ x corresponds to the following equations: x1−x2=λ x1, 2xn−xn+1 −xn−1=λ xn,∀n= 1,2, xN′−xN′−1=λ xN′.(40) Substitution of the trial solution xn= cos(q n −ϕq), leads to ϕq=q/2, q= m π/N′, with m= 0, . . . , N′−1, and λm= 4 sin2(m π/(2 N′)), i.e., positive eigenvalues without a degenerate zero. In case all the oscillators of the group are out of phase, Q→ −γQ and λm→ −γλm. Therefore: Theorem 6 If ε > 0, the time–reversible, in-phase, multibreathers are unstable (stable) with soft (hard) on–site potential. For any out–of–phase multibreathers, and ε < 0the conclusions are reversed. Although, we have no mathematical proof, according to numerical calculations of the eigenvalues of the perturbation matrices corresponding to groups with different codes, the numbers of negative and positive eigenvalues are equal to the numbers of −1 and +1 in {σnσn+1}N′−1 n=1 . Therefore being unstable for any on–site potential. Note that if there are several non–contiguous multibreathers, the perturbation matrix if composed of independent blocks of perturbation submatrices of the same type, each one with a zero eigenvalue by construction as seen in subsection 4.1. The stability of such a system cannot be determined by the MST. The only possible conclusion is that if there are eigenvalues of different sign this system would be unstable. Phonobreathers are multibreathers without oscillators at rest. For time–reversible phonobreathers in a system either with free–ends, or fixed-ends boundary conditions, the perturbation matrix Qis the same as in Eq. (39). If the boundary conditions are periodic, the only change is that Q1 1 =QN,N = 2, QN1=Q1,N =−1, and the eigenvalues are the same. That is, the mode q= 0 with attractive coupling and hard potential is stable. Changing the hardness, 20 the type of coupling, or the mode, changes each time the stability. In this way, we reproduce the results in Refs. [7,6] as a consequence of modulational instability and using the DNLS approximation, respectively, and in theorem 9 in Ref. [3], based in the properties of the action. Parity instabilities are a interesting effect for phonobreathers with periodic boundary conditions. It is a consequence of the parity of the total number of sites N. If, for example, Vis soft, the coupling is attractive, and Nis even, all the eigenvalues of the perturbation matrix are negative, i.e., the system is stable. But if Nis odd, there is an isolated positive eigenvalue, bringing about the instability. There is a clear physical origin of this phenomenon, as the boundaries are equivalent to an inhomogeneity. According to our theory it is also obvious: we have, in fact two oscillators at the ends coupled and in phase, which gives rise to a positive eigenvalue. Dark breathers are multibreathers with only one or a few oscillators at rest. Their stability depends on the characteristics of the contiguous groups of excited oscillators, as a consequence of theorem 6. We have to be careful about the number of particles, for example the πmode with code around the dark site [. . . , 1,−1,0,1,−1, . . . ], needs an odd number of sites to avoid having two in–phase oscillators at the ends, and a parity instability. The one with [. . . , −1,1,0,1,−1, . . . ] needs an even number of sites. If the system has free of fixed periodic conditions, the perturbation matrix is decoupled in two and the stability is undefined as commented above. We have compared the predictions of the instabilities and the values of the intersection points of the bands with the axis θ= 0, in many cases, as for example in Ref. [10], with very good results. The question of up to which values of the coupling parameter εthe theory is valid, is however unsolved. Generically speaking values up to ε= 0.1 (compared to the rest frequency ω0= 1), are safe, but there are some exceptions. In the same reference, for the dark breather hard on–site potential and attractive coupling, a mode not predicted by the MST gives rise to an instability as soon as ε= 0.022. As the width of the phonon band can be calculated, and the width of the background is given by the eigenvalues of the perturbation matrices, it is also possible to predict the values of the ε, for which oscillatory and subharmonic instabilities occur, but we do not extend further here. 9 Summary and conclusions We have developed using degenerate perturbation theory a method for obtaining the stability properties of multibreathers of any code at low coupling. 21 The method is synthesized in three different versions of a theorem, referred to as the multibreather stability theorem or MST for short. The two simplest versions correspond to time–reversible multibreathers, with linear coupling and homogeneous on–site potentials. If the on–site potential is symmetric or we consider only in–phase multibreathers, it involves only some simple modifications of the coupling matrix. If none of those conditions are fulfilled then it involves the analytical or numerical calculations of a magnitude γ, called the symmetry coefficient, as its value is 1 for symmetric on–site potentials. For soft potentials and values of the frequency until about one half of the rest frequency, a good approximation is the multibreather frequency, result which we demonstrate is exact for the Morse potential. For hard potentials its value its close to unity. The generalized version of the MST is of not so straightforward applicability, although its complexity depends on the characteristics of the system and multibreather considered. We give some examples of application of the method and compare it with numerical results. The systems considered are two and three–site breathers, with nearest neighbor or long–range interactions, multibreathers, phonobreathers and dark breathers. A parity instability can appear in finite systems depending on the parity of the number of oscillators. All these examples are interesting in themselves, but also illustrate the use of the MST, and show that it provides a powerful method of easy applicability to determine multibreathers stability. Some other applications under study are multibreathers in 2D and 3D systems, prediction of the oscillatory and subharmonic instabilities and the peculiarities of some disordered systems Acknowledgements This work has been supported by the European Union under the RTN project, LOCNET, HPRN–CT–1999–00163. J Cuevas acknowledges an FPDI grant from ‘La Junta de Andaluc´ıa’. A Calculation of the symmetry parameter γfor a Morse potential In this appendix, we calculate the value of the symmetry parameter γin theorem 2 for the particular case of a Morse potential. The choice of this kind of potential relies in the fact that the expressions of the orbits are easy to manage. We will proof that, for this potential, γ=ωb. Let us suppose an isolated oscillator submitted to a Morse potential. The 22 energy of the system is: E=1 2˙x2+1 2(exp(−x)−1)2.(A.1) We look for solutions with period T= 2π/ωbthat are time–reversible, i.e, with ˙x(0) = 0, or, in other words, at t= 0 the system is at a turning point. There are two of them, obtained from the equation above: x1=−log(1 + √2E)<0 and x2=−log(1−√2E)>0. If we consider the solution with x(0) = x1, t(x) is given by: t(x) = ∫x x1 dx √2(E −(exp(−x)−1)2/2) ,(A.2) which leads to: t(x) = 1 √1−2E(π 2−arcsin [1 √2E(1 + 2E −1 exp(−x))]).(A.3) By substitution of t=T/2 = π/ωband x=x2we obtain ωb=√1−2E. The inversion of Eq. (A.2) leads to: x(t) = log 1−√1−ω2 bcos ωbt ω2 b .(A.4) To calculate γ, we need the Fourier coefficients of x(t). If the solution x(t) is expressed as: x(t) = z0+ 2 ∞ ∑ k=1 zkcos(k ωbt),(A.5) the coefficients {zk}∞ k=0 are given by: zk=1 T∫T/2 −T/2 x(t) cos(k ωbt) dt . (A.6) The calculation of these integrals is straightforward and leads to: z0= log 1 + ωb 2ω2 b ;zk=−1 n(1−ωb √1−ω2 b)k .(A.7) The parameter γis given by: γ=−∑k≥1(−1)kk2z2 k ∑k≥1k2z2 k .(A.8) By substitution of the Fourier coefficients in Eq. (A.7), we obtain γ=−∑k≥1(−r)k ∑k≥1rkwith r=1−ωb 1 + ωb .(A.9) 23 As |r|<1, γcan be easily calculated using the expressions for the sum of an infinite geometric progression: ∑ k≥1 rk=r 1−r,∑ k≥1 (−r)k=−r 1 + r(A.10) leading to γ=ωb, as we wanted to proof. B Curvature of the bands In this section we demonstrate a key point used in our theory: the fact that the band corresponding to an isolated excited oscillators has negative curvature if the on site potential is hard and positive if it is soft. Unfortunately, the demonstration is quite long if we want it to be self–consistent. A similar one can be found in Ref. [11] and different ones in Ref. [5,12]. Let us consider the dynamical equation of an isolated oscillator with Hamiltonian H= 1/2p2+V(u), with ua single–value, real function of time and p= ˙u: ¨u+V′(u) = 0 (B.1) The characteristic equation for the Newton operator Ncorresponding to a given periodic solution of Eq. (B.1)is Nξ≡¨ ξ+V′′(u)ξ=Eξ , (B.2) where ξis a C2function of t. For each eigenvalue Ethere are only 2 independent solutions, determined, for example, by the values of ξ(0) and ˙ ξ(0). We can obtain them for E= 0, by deriving Eq. (B.1) with respect to tand with respect to ωb: ... u+V′′(u) ˙u= 0 ⇔ N ˙u= 0 ; ¨uw+V′′(u)uw= 0 ⇔ Nuw= 0 .(B.3) As we have a single oscillator we can suppose that u(t) is time symmetric with a suitable origin of time, i.e., ˙uis a time-antisymmetric function of twith the same period T= 2 π/ωbas u(t): ˙u(t+T) = ˙u(t), ˙u(0) = ˙u(T) = 0. The properties of uw=du/dωb, can be also obtained. We can write uas a cosine Fourier series u(t) = ∑∞ k=0 ckcos(ωbt), with the coefficients ckdepending on ωb. Therefore, du/dωb=∑∞ k=0(dck/dωb) cos(ωbt) + ∑∞ k=0 −k t sin(ωbt), or: uw≡du dωb =γ(t) + t ωb ˙u(t) (B.4) γ(t) being a time symmetric, periodic function of time with the same period Tas u. This functions constitute a base of the eigenspace corresponding to E= 0. 24 On the other hand, N(B.2) is a periodic operator because it depends on a periodic function u. Its eigenfunctions are given by the Bloch theorem: ξ(θ, t) = χ(iθ, t) exp(θt/T),(B.5) with χa T–periodic function of time and θ=θ(E). If we consider only values of Efor which θis real, the set of points (θ, E(θ) constitutes a band. If ξis a solution of Nξ=Eξ, the fact that Nis a real operator implies that ξ∗is another solution, i.e., (−θ, E) also belongs to the band and this is symmetric with respect to θand, therefore, dE/ dθ(0) = 0. As Nis time–symmetric ξ(θ, −t) is another, and has to be proportional to ξ∗as there can be only two independent solutions. The Floquet matrix FEis defined as:    ξ(T) ˙ ξ(T)  =FE   ξ(0) ˙ ξ(0)   .(B.6) The eigenvalues {λl}of FEare the Floquet multipliers, if we write them as λl= exp iθl, with θreal or complex, iθlare the Floquet exponents, and θl, the Floquet arguments. An eigenfunction with θlreal is bounded, if θl= 1 the eigenfunction has period T. The condition for linear stability of the solution uis that all the Floquet arguments of F0are real. It is easy to check some properties. The Bloch functions Eq. (B.5) are eigenfunctions of FEwith Floquet exponent θ, therefore FEis diagonalizable (over C) with Floquet multipliers exp(±iθ) and Floquet arguments ±θ. However, at E= 0, there is only one eigenfunction of F0, ˙u, the other independent function of the subspace Nξ= 0, uw, it is not an F0–eigenfunction. That is F0has a degenerate Floquet multiplier 1, or Floquet argument 0 and can be transformed into a Jordan block. Another consequence of the latter is that (θ, E) = (0,0) belongs to the band. We need to calculate the curvature of the band at this point. Let us come back to the Bloch functions ξ(θ, t) in Eq. (B.5) and ξ∗(θ, t) corresponding to a pair of symmetric points (±θ, E) of the band, and suppose that θ→0 (and then E→0). At θ= 0 the two functions collide in ξ(0, t) = χ(0, t) = χ∗(0, t), therefore ξ(0, t) is a real, T–periodic function. But we already now that at E= 0, there is a T–periodic function solution of Nξ= 0, ˙u. Therefore, ξ(0, t) = ˙u(Nis lineal, so we can adjust the norm of ξin order that the latter equation is fulfilled). We can reobtain the missing function that spans ker(N). If we derive Eq. (B.2) with respect to θ, we obtain Nξθ=Eθξ+Eξθand at E= 0 (Eθ(0) = 0) we 25