Breathers in inhomogeneous nonlinear lattices : an analysis via centre manifold reduction

Breathers in inhomogeneous nonlinear lattices : an analysis via centre manifold reduction. Guillaume Jamesa)∗ , Bernardo S´anchez-Reyb)and Jesus Cuevasb) a)Institut National Polytechnique de Grenoble and CNRS, Laboratoire Jean Kuntzmann (UMR 5224), tour IRMA, BP 53, 38041 Grenoble Cedex 9, France. e-mail : [email protected] b)Grupo de F´ısica No Lineal, Universidad de Sevilla, Departamento de F´ısica Aplicada I, Escuela Universitaria Polit´ecnica, c/. Virgen de ´ Africa 7, 41011-Sevilla, Espa˜na. e-mail : b[email protected], jcuev[email protected] October 7, 2008 Abstract We consider an infinite chain of particles linearly coupled to their nearest neighbours and subject to an anharmonic local potential. The chain is assumed weakly inhomogeneous, i.e. coupling constants, particle masses and on-site potentials can have small variations along the chain. We look for small amplitude and time-periodic solutions, and in particular spatially localized ones (discrete breathers). The problem is reformulated as a nonautonomous recurrence in a space of time-periodic functions, where the dynamics is considered along the discrete spatial coordinate. Generalizing to nonautonomous maps a centre manifold theorem previously obtained for infinite-dimensional autonomous maps [Jam03], we show that small amplitude oscillations are determined by finite-dimensional nonautonomous mappings, whose dimension ∗Corresponding author 1 depends on the solutions frequency. We consider the case of two-dimensional reduced mappings, which occurs for frequencies close to the edges of the phonon band (computed for the unperturbed homogeneous chain). For an homogeneous chain, the reduced map is autonomous and reversible, and bifurcations of reversible homoclinic orbits or heteroclinic solutions are found for appropriate parameter values. These orbits correspond respectively to discrete breathers for the infinite chain, or “dark” breathers superposed on a spatially extended standing wave. Breather existence is shown in some cases for any value of the coupling constant, which generalizes (for small amplitude solutions) an existence result obtained by MacKay and Aubry at small coupling [MA94]. For an inhomogeneous chain the study of the nonautonomous reduced map is in general far more involved. Here this problem is considered when the chain presents a finite number of defects. For the principal part of the reduced recurrence, using the assumption of weak inhomogeneity, we show that homoclinics to 0 exist when the image of the unstable manifold under a linear transformation (depending on the defect sequence) intersects the stable manifold. This provides a geometrical understanding of tangent bifurcations of discrete breathers commonly observed in classes of systems with impurities as defect strengths are varied. The case of a mass impurity is studied in detail, and our geometrical analysis is successfully compared with direct numerical simulations. In addition, a class of homoclinic orbits is shown to persist for the full reduced mapping and yields a family of discrete breathers with maximal amplitude at the impurity site. 1 Introduction It is now well established that many nonlinear networks of interacting particles sustain time-periodic and spatially localized oscillations commonly denoted as discrete breathers. In spatially periodic systems, breathers are also called intrinsically localized modes [ST88] in distinction to Anderson modes of disordered linear systems [And58]. The properties of discrete breathers have been analyzed in an important number of numerical works (see the reviews [FW98, VMZ03, DLHMS04]) and their existence in periodic systems has been proved analytically in different contexts, see [MA94, Aub98, SM97, AS98, AKK01, Pan05, Fla95, Jam03] and references therein. In the context of numerical simulations or experiments discrete breathers often denote a larger class of spatially localized oscillations, such as metastable states, oscillations with a certain degree of periodicity, or even chaotic oscillations interacting with a noisy extended background [IKSF04, GLC05]. Nonlinear waves of this type are now actually detected in real materials [SS04, Swa99, SES99, EH02, Man06] and also generated in artificial systems such as Josephson junction arrays, micromechanical cantilever arrays and coupled optical waveguides (see references in [CFK04]). 2 They are thought to play a role in various physical processes such as the formation of local fluctuational openings in the DNA molecule [PB89, Pey04], which occurs in particular during thermal denaturation experiments. Beyond spatially periodic systems, it is a fundamental and challenging problem to understand breather properties in nonlinear and inhomogeneous media, such as non-periodic or disordered crystals, amorphous solids and biological macromolecules. For example the interplay between nonlinearity and disorder can provide an alternative interpretation for slow relaxation processes in glasses [KA99-00]. In quasione-dimensional media, moving localized waves interacting with impurities [KBS94, CPAR02, FPM94], extended defects [TP96] or local bends of the lattice (see [CK04] and its references) can remain trapped and release vibrational energy at specific sites. The modelling of thermal denaturation of DNA and the analysis of its local fluctuational openings, also known as denaturation bubbles, represents another problem where heterogeneity is important. In order to describe these phenomena, a nonlinear model at the scale of the DNA base pair has been introduced by Peyrard and Bishop [PB89] and further improved by Dauxois et al [DPB93]. The model describes the stretching xn(t) of the H-bonds between two bases, in the nth base pair along a DNA molecule (a large value of xncorresponding to a local opening). Each bond fluctuates in an effective anharmonic potential Vand interacts with its nearestneighbours. The model is described by an Hamiltonian system, and can be coupled with a thermostat to study the effect of thermal noise in denaturation experiments. This model accurately describes the thermal denaturation of certain real DNA segments provided their heterogeneity is taken into account [CG98]. In its simplest form, the model incorporates different dissociation energies for the adenine-thymine (AT) and guanine-cytosine (GC) base pairs. The Hamiltonian of the system reads H= +∞ X n=−∞ m 2˙x2 n+Vsn(xn) + k 2(1 + ρ e−β(xn+1+xn)) (xn+1 −xn)2,(1) where Vsn(x) = Dsn(1 −e−asnx)2is a Morse potential depending on the base pairs sequence sn∈ {AT, GC}. The case ρ= 0 yields a particular case of a Klein-Gordon lattice, i.e. the model consists in a chain of anharmonic oscillators with harmonic nearest-neighbours coupling. For parameters corresponding to real DNA sequences, Langevin molecular dynamics of (1) have shown that some locations of discrete breathers heavily depend on the sequence and seem to coincide with functional sites in DNA [KRBCU04], but at the present time this topic remains controversial [vECLP06]. ¿From a mathematical point of view, Albanese and Fr¨ohlich have proved the existence of breathers for a class of random Hamiltonian systems describing an infinite array of coupled anharmonic oscillators [AF91] (see also the earlier work [FSW86] of 3 Fr¨ohlich et al concerning quasiperiodic localized oscillations). These breather families can be parametrized by the solutions frequencies, which belong to fat Cantor sets (i.e. with nonzero Lebesgue measure) of asymptotically full relative measure in the limit of zero amplitude. These solutions are nonlinear “continuations” of a given Anderson mode from the limit of zero amplitude, and the gaps in their frequency values originate from a dense set of resonances present in the system. For disordered Klein-Gordon lattices, complementary numerical results on the continuation of breathers with respect to frequency or the transition between breathers to Anderson modes are available in [KA99-00, AMM99]. In addition, the existence of breathers in inhomogeneous Klein-Gordon lattices (with disordered on-site potentials) has been proved by Sepulchre and MacKay [SM98, SM97] for small coupling k. The proof is based on the continuation method previously introduced by MacKay and Aubry [MA94] for an homogeneous chain (method considerably generalized in [SM97]). For k= 0 the system reduces to an array of uncoupled non-identical anharmonic oscillators, and the simplest type of discrete breather consists of a single particle oscillating while the others are at rest. Under a nonresonance condition [SM98, SM97], this solution can be continued to small values of k(in most cases at fixed frequency) using the implicit function theorem, yielding a spatially localized solution. In this paper we provide complementary mathematical tools for studying timeperiodic oscillations (not necessarily spatially localized) in inhomogeneous infinite lattices. The theory is developed in a very general framework, and applied to breather bifurcations in inhomogeneous Klein-Gordon lattices as lattice parameters and breather frequencies are varied. We start from a general Klein-Gordon lattice with Hamiltonian H= +∞ X n=−∞ Mn 2˙x2 n+DnV(Anxn) + Kn 2(xn+1 −xn)2(2) (case ρ= 0 of (1) with more general inhomogeneities). The potential Vis assumed sufficiently smooth in a neighbourhood of 0 with V0(0) = 0, V00(0) = 1. The general theory is a priori valid for small inhomogeneities and small amplitude oscillations. In particular, in our application to system (2) we assume Mn,Dn,An,Knto be close (uniformly in n) to positive constants. However, considering an example of Klein-Gordon lattice with a mass defect, we check using numerical computations that our tools remain applicable up to strongly nonlinear regimes, and sometimes for a large inhomogeneity. Our analysis is based on a centre manifold reduction and the concept of spatial dynamics. This concept was introduced by K. Kirchg¨assner [Kir82] for nonlinear elliptic PDE in infinite strips, considered as (ill-posed) evolution problems in the unbounded space coordinate, and locally reduced to a finite-dimensional ODE on 4 an invariant centre manifold. This idea was transposed to the context of travelling waves in homogeneous infinite oscillator chains by Iooss and Kirchg¨assner [IK00], and centre manifold reduction has been subsequently applied to the analysis of travelling waves and pulsating travelling waves in different one-dimensional homogeneous lattices [Ioo00, JS05, IJ05, Sir05, PR05, IP06]. Indeed, looking for travelling waves in an oscillator chain yields an advance-delay differential equation (a system of such equations in the case of pulsating travelling waves), which can be reformulated as an infinite-dimensional evolution problem in the moving frame coordinate, and locally reduced to a finite-dimensional ODE under appropriate spectral conditions. In [Jam01], one of us has proved the existence of breathers in Fermi-Pasta- Ulam (FPU) lattices using a similar technique in a discrete context. The dynamical equations for time-periodic solutions were reformulated as an infinite-dimensional recurrence relation in a space of time-periodic functions, and then locally reduced to a finite-dimensional mapping on a centre manifold, where breathers corresponded to homoclinic orbits to 0. A general centre manifold theorem for infinite-dimensional maps with unbounded linearized operator has been proved subsequently [Jam03] and has been used to analyze breather bifurcations in diatomic FPU lattices [JN04, JK07] and spin lattices [Nob04]. More generally, the dynamical equations of many onedimensional lattices can be reformulated as infinite-dimensional maps in loop spaces as one looks for small amplitude time-periodic oscillations ([Jam03], section 6.1). As shown in the present paper, the centre manifold reduction theorem readily applies to homogeneous Klein-Gordon lattices, where Mn=m,Dn=d,An=a, Kn=kin (2) and m, d, a, k > 0. This reduction result rigorously justifies (in the weakly nonlinear regime) a formal one-Fourier mode approximation previously introduced in [BCKRBW00]. The equations of motion read md2xn dt2+da V 0(a xn) = k(xn+1 −2xn+xn−1), n ∈Z.(3) Looking for time-periodic solutions (with frequency ω) and setting xn(t) = ˜xn(ω t), (3) can be formulated as an (ill-posed) recurrence relation (˜xn+1,˜xn) = F(˜xn,˜xn−1) in a space of 2π-periodic functions. Using the theorem of [Jam03], one can locally reduce the problem to a finite-dimensional mapping on a centre manifold whose dimension depends on the frequency ω. More precisely, equation (3) linearized at xn= 0 admits solutions in the form of linear waves (phonons) with xn(t) = Acos (qn −ωqt), whose frequency satisfies the dispersion relation mω2 q=a2d+ 2k(1 −cos q).(4) The frequencies ωqlie in a band [ωmin, ωmax] with ωmin >0. In the nonlinear case, the dimension of the centre manifold depends on how many multiples of ω belong to (or are close to) the phonon band. When ω≈ωmax or ω≈ωmin (with 5 no additional resonance), the centre manifold is two-dimensional if solutions are searched even in time, which reduces (3) locally to a two-dimensional reversible mapping on the centre manifold. For appropriate parameter values, this map admits small amplitude homoclinic solutions to 0 corresponding to breather solutions of (3). Breather solutions in this system have been proved to exist by MacKay and Aubry [MA94] for small values of the coupling parameter k. Known regions of breather existence are considerably extended here, since we prove the existence of small amplitude breathers for arbitrary values of kin some cases and for frequencies close to the phonon band edges (see theorem 7 p.32). In addition we prove the existence of “dark breather” solutions, which converge towards a nonlinear standing wave as n→ ±∞ and have a much smaller amplitude at the centre of the chain. These solutions correspond to heteroclinic orbits of the reduced two-dimensional map. Furthermore, we extend this analysis to the case when small lattice inhomogeneities are present. The dynamical equations of the inhomogeneous system (2) take the form Mn d2xn dt2+DnAnV0(Anxn) = Kn(xn+1 −xn)−Kn−1(xn−xn−1), n ∈Z,(5) and time-periodic solutions can be obtained as orbits of a nonautonomous map (xn+1, xn) = F(λn, xn, xn−1), where the nonconstant lattice parameters are embedded in a multicomponent parameter λn. Fixing Mn=m+mn,Dn=d+dn, An=a+an,Kn=k+kn, we consider the case when constant lattice parameters m, d, a, k > 0 are perturbed by uniformly small sequences (mn)n∈Z, (dn)n∈Z, (an)n∈Z, (kn)n∈Z. We prove (see theorem 3 p.19) that small amplitude time-periodic solutions with frequencies close to ωmin or ωmax are determined by a two-dimensional nonautonomous map. Moreover, we generalize this reduction result in the case when several multiples of ωare close to the band [ωmin, ωmax], which yields a higherdimensional reduced problem (see theorem 4 p.21). In fact we prove this type of reduction result in a very general framework, for infinite-dimensional mappings with small nonautonomous perturbations, considered in a neighbourhood of a non-hyperbolic fixed point, or close to a bifurcation. The linear autonomous part of the map must satisfy a property of spectral separation (see theorem 1 p.15), but a large number of one-dimensional lattices with finite-range coupling fall within this category. We obtain a direct proof of the reduction result by observing that any nonautonomous mapping un+1 =F(λn, un) can be seen as a projection of an extended autonomous mapping, to which the centre manifold theorem of [Jam03] can be applied under appropriate assumptions. The centre manifold of the extended map is infinite-dimensional, but this case is also covered in [Jam03]. The reduced nonautonomous mapping for the original system can be interpreted as a projection on a finite-dimensional subspace of the extended autonomous mapping 6 restricted to the invariant centre manifold. We use this reduction result to analyse the case when equation (5) presents a mass defect at a single site, all other lattice parameters being independent of n. In that case, the linearized problem admits a spatially localized mode (usually denoted as an impurity mode or defect mode), and a nonlinear continuation of this mode can be computed [FPM94], corresponding to a Lyapunov family of periodic orbits. Klein-Gordon systems with a coupling defect or a harmonic impurity in the on-site potential share similar characteristics [CPAR02], as well as nonlinear lattices with a different type of nonlinearity [KZK97]. In addition to this simple localization phenomenon, single impurities can have more complex effects in a nonlinear system. Indeed it is a common feature to observe a complex sequence of tangent bifurcations between (deformations of) site-centered and bond-centered breathers in some neighbourhood of the defect as the strength of an impurity is varied [CK04, KZK97]. Using numerical computations we show some examples of such bifurcations in the present paper, as one varies the strength of a mass defect in system (5). From a physical point of view it is quite important to understand how a local change in the lattice parameters modifies the set of spatially localized solutions. For example, this could contribute to explain how a mutation at a specific location of an homogeneous sequence of (artificial) DNA would modify the structure of fluctuational openings [KRBCU04]. This paper provides a qualitative explanation of such tangent bifurcations, which reveals also very precise quantitatively when compared with numerical simulations of the Klein-Gordon model. According to the previously described reduction theorem, for a small mass defect of size ², small amplitude breather solutions of (5) with frequencies below (and close to) ωmin are described by a two-dimensional nonautonomous mapping vn+1 =f(vn, ω) + ² g(n, vn, ω, ²). Here we consider the principal part of the reduced mapping as (vn, ω, ²)≈(0, ωmin,0). We show that this truncated reduced map admits an homoclinic orbit to 0 (corresponding to an approximate breather solution for the oscillator chain) if, for ²= 0, the image of the unstable manifold of 0 under a certain linear shear intersects its stable manifold. The linear shear is O(²)-close to the identity. When the on-site potential is soft (i.e. the period of small oscillations in this potential increases with amplitude), these manifolds have very complicated windings characteristic of homoclinic chaos, hence the set of their intersections changes in a complex way as the linear shear varies, or equivalently as one varies the mass defect. This phenomenon explains the existence of the above mentioned tangent bifurcations, at least for small defect sizes, and for small amplitude breathers with frequencies close to the phonon band. In addition, we show (by comparison with direct numerical simulations of the Klein-Gordon model) that this picture remains valid quite far from the weakly nonlinear regime. Let us note that, to obtain an exact solution of (5) from an orbit of the truncated map, it would be necessary to control the effect of higher order terms (with respect 7 to vn,ω−ωmin,²) present in the full reduced mapping and prove the persistence of this solution. This result is obtained for Lyapunov families of periodic orbits, which correspond here to discrete breathers with maximal amplitude at the impurity site (see theorem 8 in section 4.1.5). In that case, the corresponding orbits of the reduced mapping appear through a pitchfork bifurcation when ωreaches the linear defect mode frequency. The persistence of the above mentioned tangent bifurcations of discrete breathers is a much more complex problem, which study would require asymptotical techniques beyond all algebraic orders (more details in this respect are given in section 4.1.4). This problem is not examined here from the analytical side, but we compare instead numerically computed solutions of (5) with approximate solutions deduced from the truncated map. The very good agreement leads us to conjecture that most of the tangent bifurcations existing for the truncated problem persist for the full reduced system. Lastly we consider the more general case when system (5) admits a finite number of defects, i.e. perturbations mn,dn,an,knhave a compact support (as above these perturbations are assumed to be small, of order ²). We show that the approach developped for a single impurity can be extended to this case (see lemma 7 p.49), where the linear shear is replaced by a more general linear near-identity transformation A². The linear transformation A²provides a useful tool for studying breather bifurcations in Klein-Gordon lattices with a finite number of impurities, as for the single impurity case that we have analyzed in detail. By computing the principal part of A²as ²is small and frequencies are close to ωmin, we show that the effect of the parameter sequence (λn)n∈Zon the set of small amplitude breather solutions should mainly depend on weighted averages of the defect values. The outline of the paper is as follows. Section 2 presents the centre manifold reduction theory for time-periodic oscillations in weakly inhomogeneous nonlinear lattices. We treat the case of Klein-Gordon lattices in detail in sections 2.1 and 2.3, and formulate the reduction theory in a much more general setting in section 2.2. Section 3 concerns spatially homogeneous Klein-Gordon lattices. Existence theorems for small amplitude breather and dark breather solutions are deduced from the dynamics of two-dimensional reversible maps on invariant centre manifolds. The case of weakly inhomogeneous Klein-Gordon chains is considered in section 4, where the truncated reduced map is analyzed for a finite number of defects. A geometrical condition for the existence of homoclinic orbits to 0 is derived in section 4.2, and some homoclinic bifurcations are studied in detail in section 4.1 for a single mass defect. In the latter case, breather solutions are numerically computed in section 5 and the results are successfully compared with our analytical findings. 8 2 Reduction result for small inhomogeneities In this section we consider system (5) in the limit of small inhomogeneities. We show that all small amplitude time-periodic solutions are determined by a finitedimensional nonautonomous map, whose dimension depends on the frequency domain under consideration. For this purpose we reformulate (5) as a map in a loop space, perturbed by a small nonautonomous term (section 2.1). Then we prove in section 2.2 a general centre manifold reduction theorem for infinite-dimensional maps with small nonautonomous perturbations. This result is based on the centre manifold theorem proved in [Jam03] for autonomous systems. Our general result is applied to the inhomogeneous Klein-Gordon lattice, which yields the above mentioned reduction result (section 2.3). 2.1 The Klein-Gordon system as a map in a loop space We set xn(t) = yn(ω(k/m)1/2t) in equation (5), where ynis 2π-periodic in t(hence xnis time-periodic with frequency ω(k/m)1/2). The constant a > 0 being fixed, we also define ˜ V(x) = a−2V(ax). Equation (5) becomes ω2(1+²n)d2yn dt2+Ω2(1+ηn)˜ V0((1+γn)yn) = yn+1 −yn−(1+κn) (yn−yn−1), n ∈Z (6) where Ω2=a2d/k and 1+²n= (1+ mn m)/(1+ kn k), 1+ηn= (1+ dn d)(1+ an a)/(1+ kn k), γn=an a, 1 + κn= (1 + kn−1 k)/(1 + kn k). The sequences (²n)n∈Z, (ηn)n∈Z, (γn)n∈Z (κn)n∈Zwill be assumed sufficiently small in `∞(Z), where `∞(Z) is the classical Banach space of bounded sequences on Z, equiped with the supremum norm. To simplify the notations, we shall drop the tilde in the sequel when referring to the renormalized potential ˜ V. Moreover we shall use the shorter notations {²}when referring to sequences (²n)n∈Z. To analyze system (6) we use the same approach as in [Jam03] for spatially homogeneous systems. We reformulate (6) as a (nonautonomous) recurrence relation in a space of 2π-periodic functions of t, and locally reduce the (spatial) dynamics to one on a finite-dimensional centre manifold. We restrict our attention to the case when ynis even in tin order to deal with lower-dimensional problems. More precisely, we assume yn∈H2 #for all n∈Z, where Hn #={y∈ Hn per(0,2π), y is even }and Hn per(0,2π) denotes the classical Sobolev space of 2π- periodic functions (H0 per(0,2π) = L2 per(0,2π)). Since our analysis concerns small amplitude solutions and small inhomogeneities, the first step consists in studying the linearized system at yn= 0 when ²n, ηn, γn, κn 9 where fn∈Ck((Xc∩Ω) ×Λ, Xc)is defined by fn(., {λ}) = πc(L+N(., λn)) ◦(I+ψ(., τn{λ})). ii) Conversely, if {Yc}is a solution of (14) such that Yc n∈Ωfor all n∈Z, then Yn=Yc n+ψ(Yc n, τn{λ})satisfies (13). iii) If L+N(., λ)commutes with a linear isometry T∈ L(X)∩L(D)then Tψ(., {λ}) = ψ(., {λ})◦Tand Tfn(., {λ}) = fn(., {λ})◦T. Properties i), ii) reduce the local study of (13) to that of the nonautonomous recurrence relation (14) in the subspace Xc. Note that the dependency of ψand the reduced map fnwith respect to sequences {λ}is nonlocal. In what follows we give a simple proof of theorem 1 which relies on the fact that the nonautonomous mapping (13) can be seen as a projection of an extended autonomous mapping, to which the centre manifold theorem proved in [Jam03] can be applied. This procedure will explain why the result of theorem 1 can be seen as a centre manifold reduction, since the reduction function ψwill appear as one component of the function having the centre manifold as its graph. The reduced nonautonomous mapping (14) will be interpreted as a projection of the extended autonomous mapping restricted to the invariant centre manifold. Theorem 1 has been proved in [Jam03] in the case of an autonomous mapping, when the sequence {λ}is absent or replaced by a simple parameter λ∈Rp. To recover this autonomous case we introduce the additional variable Sn= τn{λ} ∈ E. Note that for any fixed n∈Z,Sndenotes a bounded sequence in Rp(to simplify the notations we use the symbol Sninstead of {Sn}). Given a sequence {λ} ∈ Ewe also note δ0{λ}=λ0. Equation (13) can be rewritten Yn+1 =L Yn+N(Yn, δ0Sn), Sn+1 =τ1Sn,(15) which consists in an autonomous mapping in X×E. In what follows we apply the theory of [Jam03] to system (15). As we shall see the corresponding centre manifold will be infinite-dimensional due to the second component of (15). The case of infinite-dimensional centre manifolds has been treated in [Jam03], with the counterpart that theory is restricted to maps in Hilbert spaces. Consequently the first step is to search for Snin a suitable Hilbert space instead of the Banach space E. For this purpose we consider the space of sequences h−1={{u}/ uk∈Cp,k{u}k−1<+∞}, where k{u}k2 −1=Pk∈Z(1 + k2)−1kukk2. The space h−1defines a Hilbert space equiped with the scalar product h{u},{v}i =Pk∈Z(1 + k2)−1uk·vk, where ·denotes the usual scalar product on Cpand kk the associated norm. For all n∈Zwe 16 now search for Snin the space H=h−1∩(Rp)Zconsisting of real sequences in h−1. Note that E⊂H, the embedding being continuous. Since sequences in Hmay be unbounded and N(Y, .) is defined on a neighbourhood Vof λ= 0 in Rp, we replace (15) by a locally equivalent problem (Yn+1, Sn+1) = F(Yn, Sn) (16) where F(Y, S) = ¡L Y +N(Y, γ(δ0S) ) , τ1S¢, γ:Rp→ V is a C∞cut-off function satisfying kγ(x)k ≤ kxk,γ(x) = xfor kxk< r, γ(x) = 0 for kxk>2r,rbeing chosen small enough (with B(0,2r)⊂ V). Problem (16) consists in an autonomous mapping in X×H. In order to apply the centre manifold theorem of [Jam03] we need to study the spectrum of DF(0) = L×τ1. One has clearly σ(DF(0)) = σ(L)∪σ(τ1), where σ(τ1) is determined in the following lemma. Lemma 1 The spectrum σ(τ1)of τ1:H→Hconsists of the unit circle. Proof. Consider the complexification h−1of H. Given a sequence {f} ∈ h−1and z∈C, we look for {u} ∈ h−1satisfying (zI −τ1){u}={f}.(17) Equation (17) can be solved in a simple manner using Fourier series. Recall that the periodic Sobolev space H1 per(0,2π) can be defined as the set of functions in L2(R/2πZ,Cp) whose Fourier coefficients form a sequence in h1, where h1={{u}/ uk∈Cp,X k∈Z (1 + k2)kukk2<+∞}. In the same way its dual space H−1 per(0,2π) is isomorphic to h−1, where the isomorphism C:H−1 per(0,2π)→h−1is again given by Cn(T) = 1 2πhT, e−intifor all T∈H−1 per(0,2π). In addition one has the useful property τ1C(T) = C(e−it T) for all T∈H−1 per(0,2π). Now return to equation (17) and consider T=C−1({u}) and S=C−1({f}). One obtains the equivalent problem in H−1 per(0,2π) (z−e−it)T=S. (18) If |z| 6= 1 then (18) has the unique solution T= (z−e−it)−1S, hence z /∈σ(τ1). If z=eiθ is chosen on the unit circle, T= 2πδ−θis a solution for S= 0, corresponding to an eigenvector {u}={einθ}of τ1. 2 17 As it follows from lemma 1, σ(DF(0)) consists of the union of σ(L) with the unit circle. Consequently DF(0) has the property of spectral separation, i.e. the hyperbolic part of its spectrum is isolated from the unit circle. Moreover the centre subspace of DF(0) is simply Xc×H. With these spectral properties at hand, we now apply the centre manifold theorem of [Jam03] which states the following. Theorem 2 There exists a neighbourhood Ωט Λof (Y, S) = 0 in D×Hand a map ψ∈Ck(Xc×H, Dh)(with ψ(0,0) = 0,Dψ(0,0) = 0) such that the manifold M={(Y, S)∈D×H / Y =Yc+ψ(Yc, S), Y c∈Xc} has the following properties. i) Mis locally invariant under F, i.e. if (Y, S)∈ M∩(Ω ט Λ) then F(Y, S)∈ M. ii) If {(Y, S)}is a solution of (16) such that (Yn, Sn)∈Ωט Λfor all n∈Z, then (Yn, Sn)∈ M for all n∈Z(i.e. Yh n=ψ(Yc n, Sn)) and (Yc n, Sn)satisfies the recurrence relation in Xc×H Yc n+1 =˜ f(Yc n, Sn), Sn+1 =τ1Sn,(19) where ˜ f(Yc, S) = πc[L+N(., γ(δ0S))](Yc+ψ(Yc, S)). iii) Conversely, given a solution {(Yc, S)}of (19) such that (Yc n, Sn)∈Ωט Λfor all n∈Z, consider Yn=Yc n+ψ(Yc n, Sn). Then (Yn, Sn)defines a solution of (16) lying on M. iv) If L+N(., λ)commutes with a linear isometry T∈ L(X)∩L(D)then Tψ(Yc, S) = ψ(TY c, S)and T˜ f(Yc, S) = ˜ f(TY c, S). The manifold Mis called a local Ckcentre manifold for (16). It is locally invariant under F(as stated by property i) ) and the linear isometries of (16). Property iv) expresses the invariance of Munder the linear isometry T×Iof (16). Now the proof of theorem 1 follows directly from theorem 2. Since Eis continuously embedded in H,ψdefines a Ckmap from Xc×Einto Dh. In theorem 1 we choose Λ as a ball of centre 0 in Esuch that Λ ⊂˜ Λ and γ=Ion Λ. Then problems (13) and (16) are equivalent for all {λ} ∈ Λ, with Sn=τn{λ}, and properties i)-ii)-iii) of theorem 1 are directly deduced from properties ii)-iii)-iv) of theorem 2. In addition, since (0, τn{λ}) is a solution of (16) for all {λ} ∈ Λ it follows ψ(0, τn{λ}) = 0 (by property ii) of theorem 2), and consequently ψ(0,{λ}) = 0. 18 2.3 Application to the Klein-Gordon lattice 2.3.1 Reduction result In this section we apply the reduction theorem 1 to the inhomogeneous Klein-Gordon lattice considered in section 2.1. We recall that the inhomogeneous system (6) has been reformulated as a nonautonomous map in a loop space given by expression (11). All parameters (sequences of heterogeneities and frequency shift µ) are embedded in the multicomponent parameter {λ}= (µ, {²},{η},{γ},{κ})∈E=R×(`∞(Z))4. The problem has exactly the general form (13) (in the particular case in which the first component of {λ}is constant) and consequently the reduction theorem 1 can be applied to (11). This yields the reduction result for the original system (6) stated below (theorem 3). It is straightforward to check that system (6) has the reduction properties i) and ii) described in this theorem since the equivalent system (11) satisfies properties i) and ii) of theorem 1 (see remark 1 p.15). However there remains to compute the explicit forms (21) and (24) of the recurrence relations given below. These expressions do not simply correspond to the two-dimensional mapping (14) rewritten as a second order recurrence relation. In addition we rewrite (14) in normal form, i.e. we perform a polynomial change of variables which simplifies (14) by keeping only its essential terms. These computations will be the object of the next three sections 2.3.2, 2.3.3 and 2.3.4. Property iii) below is equivalent to property iii) of theorem 1, where the symmetry Tis the half period time shift which satisfies T|Xc=−I. Theorem 3 Fix ω2=ω2 c+µin equation (6), where ωc=√4+Ω2or ωc= Ω (in that case we further assume Ω>2/√3). There exist neighbourhoods U,Vand Wof 0in H2 #,Eand Rrespectively, and a Ckmap φ:R2×E→H2 #(with φ(0,{λ}) = 0, Dφ(0,0) = 0) such that the following holds for all {λ} ∈ V. i) All solutions of (6) such that yn∈ U for all n∈Zhave the form yn(t) = βncos t+ [ φ(βn−1, βn, τn{λ}) ](t). For ωc= Ω,βnsatisfies a recurrence relation βn+1 −2βn+βn−1=Rn(βn−1, βn,{λ}) (20) where Rn:W2×V → Ris Ck. The principal part of Rnreads Rn(α, β, {λ}) = (Ω2ηn(1+γn)−(Ω2+µ)²n+Ω2γn−µ)β+κn(β−α)+B β3+h.o.t., (21) B=Ω2 8(V(4)(0) −5 3(V(3)(0))2).(22) 19 For ωc=√4 + Ω2one has βn+1 + 2βn+βn−1=Rn(βn−1, βn,{λ}),(23) with Rn(α, β, {λ}) = (Ω2ηn(1+γn)−(4+Ω2+µ)²n+Ω2γn−µ)β+κn(β−α)+ ˜ B β3+h.o.t., (24) ˜ B=Ω2 8(V(4)(0) + (V(3)(0))2(Ω2 16 + 3Ω2−2) ).(25) In both cases, higher order terms in Rnare O(k(α, β)k3k{λ}kE+k(α, β)k5)and non-local in {λ}. ii) If βnis a solution of problem (20) or (23) (respectively for ωc= Ω or ωc= √4+Ω2), such that βn∈ W for all n∈Z, then yn(t) = βncos t+φ(βn−1, βn, τn{λ}) satisfies equation (6). iii) The functions φand Rnhave the following symmetries φ(−α, −β, {λ}) = Tφ(α, β, {λ}), Rn(−α, −β, {λ}) = −Rn(α, β, {λ}), where Tdenotes the half period time shift [Tφ(.)](t) = [φ(.)](t+π). A possible way of computing the reduced recurrence relations (20) and (23) would be to consider the equivalent autonomous mapping (15) and use a classical computation scheme for centre manifolds of autonomous systems (see e.g. [Van89] for a description of the method). The first step consists in computing the Taylor expansion of the reduction function ψup to a given order. This can be done using a nonlocal equation for Yn(obtained by expressing Ynin (15) as a function of N(Yn, δ0Sn)) and computing the Taylor coefficients of ψby induction (see [Van89]). The second step is to compute the reduced recurrence relation (19) which is completely determined by ψ. In the next three sections we shall use a different method yielding simpler computations. Firstly we compute the expressions (21) and (24) in the autonomous case {²}={η}={γ}={κ}= 0, using the method of [Jam03]. Then, using a symmetry argument, we deduce how the leading order part of the reduced equation is modified by the nonautonomous terms of (11). To end this section we point out a generalization of theorem 3. As it follows from the analysis of section 2.1, the dimension of the centre space Xcof Aωis twice the number of multiples of ωlying within the band [Ω,(4 + Ω2)1/2]. More precisely, if Ω≤ωp ≤(4 + Ω2)1/2for all integers p∈ {p0, . . . , p1}, with no additional multiples entering the band, then the centre space is spanned by the corresponding Fourier modes (cos(pt),0), (0,cos(pt)). As above the following reduction result follows from theorem 1. 20 Theorem 4 Consider ωc>0such that ωcp∈[Ω,(4 + Ω2)1/2]for all integers p∈ {p0, . . . , p1}, with no additional multiples in this interval. Fix ω2=ω2 c+µin equation (6) and note N=p1−p0+ 1. Consider the subspace Hcof H2 #spanned by the NFourier modes cos(p0t), . . . , cos(p1t)and its complementary subspace H⊥ c consisting of orthogonal Fourier modes. There exist neighbourhoods U,Vof 0in H2 #,Erespectively, and a Ckmap φ:R2N×E→H⊥ c(with φ(0,{λ})=0, Dφ(0,0) = 0), such that for all {λ} ∈ V, all solutions of (6) such that yn∈ U for all n∈Zhave the form yn(t) = p1 X p=p0 [β(p) ncos (pt) ] + φ(β(p0) n−1, β(p0) n, . . . , β(p1) n−1, β(p1) n, τn{λ}).(26) Moreover, all small amplitude solutions of (6) are determined by a finite-dimensional recurrence relation obtained by projecting (6) on Hcand using the ansatz (26). In the following sections 2.3.2, 2.3.3 and 2.3.4, we compute the explicit forms of the reduced recurrence relations given in theorem 3. 2.3.2 Homogeneous case near the lower phonon band edge In this section we restrict our attention to the case when Ω >2/√3 and ω≈ω2= Ω. We consider the autonomous case when {²}={η}={γ}={κ}= 0. Equation (11) now reads Yn+1 =L Yn+N(Yn, µ), n ∈Z(27) where L=AΩis given by (9) and N((z, y), µ) = ³0, µ d2y dt2+W(y)´, with W(y)=Ω2(V0(y)−y). System (27) is a reformulation of the equations of motion for the homogenous Klein-Gordon lattice ω2d2yn dt2+ Ω2V0(yn) = yn+1 −2yn+yn−1, n ∈Z.(28) As in the nonautonomous case (11), system (27) is invariant under the symmetry T Y =Y(·+π). Moreover, the invariance yn→y−nof (28) implies that (27) is reversible with respect to the symmetry R(z, y) = (y, z), i.e. if Ynis a solution then also RY−n. In other words, if Yand [L+N(., µ)](RY ) are in some neighbourhood of 0 in Done has (L+N(., µ)◦R)2Y=Y. Lastly, due to the existence of the additional symmetry T, it is worthwhile to notice that TR defines an other reversibility symmetry. 21 In what follows we use the notations introduced in section 2.2. We recall that the spectrum of L=AΩon the unit circle consists in a double non semi-simple eigenvalue +1, and the associated two-dimensional invariant subspace Xcis spanned by the vectors Vz= (cos t, 0), Vy= (0,cos t), with L|Xc=µ0 1 −1 2 ¶ in the basis (Vz, Vy). For µin some neighbourhood Λ of 0, (27) admits a Cktwodimensional local centre manifold Mµ⊂D(which can be written as a graph over Xc), locally invariant under L+N(., µ) (see [Jam03], theorem 1 p. 32). One can write Mµ={Y∈D / Y =aVz+bVy+ψ(a, b, µ),(a, b)∈R2},(29) where ψ∈Ck(R2×Λ, Dh) and ψ(a, b, µ) = O(k(a, b)k2+k(a, b)k|µ|). Moreover, Mµis invariant under Tand R(see [Jam03], theorem 2 p. 34 and section 5.2). In the sequel we use the notations P∗(y) = 1 πR2π 0y(t) cos t dt,Pcy=P∗(y) cos t and H2 h={y∈H2 #/ P∗(y) = 0 }. The spectral projection πcon Xcreads πc(z, y) = (Pcz, Pcy) and we have Dh=H2 h×H2 h. Since Mµis invariant under Rand Vz,Vyare exchanged by R, we have the symmetry property Rψ(a, b, µ) = ψ(b, a, µ). Consequently, the function ψhas the form ψ(a, b, µ) = (ϕ(b, a, µ), ϕ(a, b, µ)) (30) with ϕ∈Ck(R2×Λ, H2 h). Since Mµis invariant under Tand T|Xc=−Iwe have in addition Tϕ(a, b, µ) = ϕ(−a, −b, µ).(31) For µ≈0, the centre manifold Mµcontains all solutions Ynof (27) staying in a sufficiently small neighbourhood of Y= 0 in Dfor all n∈Z. Their coordinates (an, bn) on Mµare thus given by a two-dimensional mapping which determines all small amplitude solutions when µ≈0. The reduced mapping is given by µan+1 bn+1 ¶=fµµan bn¶(32) where fµµa b¶=µb, −a+ 2b+r(a, b, µ)¶,(33) r(a, b, µ) = −µ b +P∗W(bcos t+ϕ(a, b, µ)).(34) One obtains equation (32) using the fact that zn=ancos t+ϕ(bn, an, µ), yn=bncos t+ϕ(an, bn, µ) 22 for Yn= (zn, yn)∈ Mµand applying P∗to equation (27) (one has P∗ϕ= 0 and P∗◦d2 dt2=−P∗on H2 #). Since the reduced mapping inherits the symmetries of (27) [Jam03], fµcommutes with T|Xc=−Iand thus r(−a, −b, µ) = −r(a, b, µ). Moreover, fµis reversible with respect to the symmetry R(a, b) = (b, a), i.e. (fµ◦R)2=I. This yields the identity r(a, b, µ) = r(−a+ 2b+r(a, b, µ), b, µ). This imposes the following structure for the Taylor expansion of rat (a, b, µ) = 0 r(a, b, µ) = −bµ +c1b3+c2ab2−1 2c2a2b+O(|b|(|a|+|b|)4+|b|(|a|+|b|)2|µ|),(35) where coefficients c1, c2have to be determined. Note that r(a, 0, µ) = 0 (see [Jam03] p.53 for details). For determining the unknown coefficients of (35), we first compute the leading order terms in the Taylor expansion of ψat (a, b, µ) = 0. This can be done using the fact that Mµis locally invariant under L+N(., µ) (see [Jam03], theorem 1 p. 32). For (a, b)≈0, this yields πh[L+N(., µ)] ( (a, b) cos t+ψ(a, b, µ) ) = ψ(fµ(a, b), µ) (36) or equivalently ϕ(−a+ 2b+r(a, b, µ), b, µ) = ϕ(a, b, µ),(37) ϕ(b, −a+ 2b+r(a, b, µ), µ) = (Ω2d2 dt2+2+Ω2)ϕ(a, b, µ)−ϕ(b, a, µ) + (1 −Pc) [µd2 dt2+W] (bcos t+ϕ(a, b, µ)). (38) Thanks to the symmetry property (37), the Taylor expansion of ϕat order 2 takes the form ϕ(a, b, µ) = ϕ011bµ −1 2ϕ110a2+ϕ110ab +ϕ020b2+ h.o.t.(39) By an identification procedure we now compute the coefficients ϕpqr in (39), using (38) and the expansion W(y) = Ω2(1 2V(3)(0) y2+1 6V(4)(0) y3+O(y4) ).(40) Identification at order bµ gives (d2 dt2+ 1)ϕ011 = 0, hence ϕ011 = 0 since ϕ011 ∈H2 h. Identification at order ab leads to ϕ020 =−1 4(Ω2d2 dt2+2+Ω2)ϕ110 (41) 23 and identification at order b2yields −ϕ110 + (Ω2d2 dt2−2 + Ω2)ϕ020 =−1 2Ω2V(3)(0) cos2t. (42) Reporting (41) in (42) gives (d2 dt2+ 1)2ϕ110 =2 Ω2V(3)(0) cos2t and consequently ϕ110 =1 Ω2V(3)(0) (1 + 1 9cos (2t)), ϕ020 =1 2V(3)(0) (−1 2−1 Ω2+ (1 6−1 9Ω2) cos (2t)). As a conclusion, we obtain ϕ(a, b, µ) = 1 Ω2V(3)(0) (1 + 1 9cos (2t)) (ab −1 2a2) +1 2V(3)(0) (−1 2−1 Ω2+ (1 6−1 9Ω2) cos (2t)) b2+ h.o.t.(43) We now compute the two-dimensional mapping giving the coordinates (an, bn) of the solutions on Mµ. Equation (32) can be written an+1 =bn, bn+1 −2bn+bn−1=r(bn−1, bn, µ).(44) Using (34), (43) and (40) yields in equation (35) c1=1 8Ω2V(4)(0) −(19 9+5 6Ω2)1 4(V(3)(0))2, c2=19 18(V(3)(0))2.(45) Lastly, one can write (44) in normal form using the change of variables bn= βn−c2 12β3 n. The normal form of (44) at order 3 reads βn+1 −2βn+βn−1=−µ βn+B β3 n+ h.o.t. (46) with B=c1+c2 2=Ω2 8(V(4)(0) −5 3(V(3)(0))2).(47) This yields the explicit form (20) of the reduced recurrence relation in the autonomous case {²}={η}={γ}={κ}= 0. 24 2.3.3 Homogeneous case near the upper phonon band edge In this section we consider the case ω≈ω1=√4+Ω2, in the autonomous case when {²}={η}={γ}={κ}= 0. Equation (11) takes the form (27), where L=Aω1is given by (9). The spectrum of Lon the unit circle consists in a double non semisimple eigenvalue −1, and the centre space Xcis again spanned by Vz= (cos t, 0), Vy= (0,cos t). For µ=ω2−ω2 1in some neighbourhood Λ of 0, there exists a smooth twodimensional local centre manifold Mµ⊂Dlocally invariant under L+N(., µ), T, Rand having the form (29). The function ψhaving the centre manifold as its graph has the form (30) and shares the property (31). For µ≈0, the centre manifold Mµcontains all solutions Ynof (27) staying in a sufficiently small neighbourhood of Y= 0 in Dfor all n∈Z. Their coordinates (an, bn) on Mµare then given by a two-dimensional mapping, which determines all small amplitude solutions when µ≈0. The operator Lhas the following structure in the basis (Vz, Vy) L|Xc=µ0 1 −1−2¶ and the reduced mapping is given by µan+1 bn+1 ¶=fµµan bn¶(48) where fµµa b¶=µb, −a−2b+r(a, b, µ)¶(49) and ris defined by (34). Since the reduced mapping inherits the symmetries of (27), fµcommutes with T|Xc=−Ihence r(−a, −b, µ) = −r(a, b, µ). Moreover, (48) is reversible with respect to the symmetry R(a, b) = (b, a), which yields the identity r(a, b, µ) = r(−a−2b+r(a, b, µ), b, µ). This implies r(a, 0, µ) = 0 and r(a, b, µ) = −bµ +c1b3+c2ab2+1 2c2a2b+ h.o.t,(50) where the coefficients c1, c2have to be determined. For this purpose, we first compute the leading order terms in the Taylor expansion of ψat (a, b, µ) = 0, using the fact that Mµis locally invariant under L+N(., µ). Equation (36) yields ϕ(−a−2b+r(a, b, µ), b, µ) = ϕ(a, b, µ),(51) 25 In what follows we reformulate the results with respect to the unscaled original system (3). For conciseness we only describe breather bifurcations, but conditions for dark breather bifurcations are easily deduced from theorems 5 and 6. We express the condition ˜ b > 0 of theorem 6 in a different way using the relation ˜ b=b−16 3 (V(3)(0))2 16 + 3Ω2. In addition, as the rescaled potential ˜ Vof (6) is replaced by the original potential Vof (3), coefficients band ˜ bare simply replaced by h=a−2band ˜ h=a−2˜ b. Theorem 7 Consider the Klein-Gordon lattice (3), where the on-site potential V satisfies V0(0) = 0,V00(0) = 1 and m, d, a, k > 0. Assume h=V(4)(0)−5 3(V(3)(0))26= 0and note Ω2=a2d/k,ω2 min =a2d/m,ω2 max = (a2d+ 4k)/m and H=V(4)(0) − 2(V(3)(0))2. i) If h < 0and Ω2>4/3, system (3) admits two families of breather solutions x1 n,x2 n parametrized by their frequency ω(in addition to phase shift), where ω≈ωmin and ω < ωmin. These solutions satisfy x1 −n+1 =x1 nand x2 −n=x2 nand decay exponentially as n→ ±∞. As ω→ωmin, the amplitude of oscillations and the exponential rate of decay are O(|ω−ωmin|1/2). The breather profile is a slow modulation of a linear mode with wavenumber q= 0. ii) If h > 0and Ω2>−16H/(3h), system (3) admits two families of breather solutions x1 n,x2 nparametrized by their frequency ω(in addition to phase shift), where ω≈ωmax and ω > ωmax. These solutions satisfy x1 −n+1(t) = x1 n(t+π/ω),x2 −n=x2 n and decay exponentially as n→ ±∞. As ω→ωmax, the amplitude of oscillations and the exponential rate of decay are O(|ω−ωmax|1/2). The breather profile is a slow modulation of a linear mode with wavenumber q=π. To interpret the conditions on the on-site potential Vin properties i) and ii), it is interesting to note that Vis soft for h < 0 and hard for h > 0 near the origin (i.e. the period of small oscillations in this potential respectively increases or decreases with amplitude). The condition on Ω in property i) corresponds to a nonresonance condition, i.e. it ensures that no multiple of ωlies in the phonon band [ωmin, ωmax] for ω≈ωmin. The condition on Ω in property ii) is of different nature and is equivalent (with the condition h > 0) to fixing ˜ h > 0. Discrete breathers were known to exist in Klein-Gordon lattices for small coupling kafter the work of MacKay and Aubry [MA94]. Theorem 7 considerably enlarges the domain of breather existence, with the limitation that it only provides small amplitude solutions. In particular, it is interesting to note that small amplitude breathers of property ii) exist for all values of kif H > 0. 32 4 Normal form analysis for inhomogeneous lattices According to theorem 3, small amplitude solutions of (6) in H2 #are described (for small inhomogeneities and frequencies close to the phonon band edges) by finitedimensional nonautonomous recurrence relations. In what follows we only consider the case ω≈Ω, the situation when ω≈√4+Ω2yielding similar phenomena. At leading order, the reduced recurrence relation (20) can be approximated by βn+1 −2βn+βn−1= (Ω2ηn(1+γn)−(Ω2+µ)²n+Ω2γn−µ)βn+κn(βn−βn−1)+B β3 n. (69) Different kinds of techniques can be employed to obtain homoclinic solutions of (69). One can use variational methods for asymptotically periodic sequences [PZ01] (see also [Wei99] in the homogeneous case), or proceed by perturbation near an uncoupled limit (also denoted anti-continuous or anti-integrable limit) where Ω and Bare large (see [AA90, Aub95, ABK04] and section 9 of [MA94]). Existence results of localized solutions are also available in [AF88] for disordered defect sequences. Another approach is to start from a known uniformly hyperbolic homoclinic solution in the homogeneous case, which persists for small inhomogeneities by the implicit function theorem, and obtain estimates for defect sizes allowing persistence (see the technique developed by Bishnani and MacKay [BM03]). Interesting related results on the structural stability of discrete dynamical systems under nonautonomous perturbations can be found in [Fra74]. With a different point of view, we develop here a dynamical system technique, valid for a finite number of defects, which allows to analyze bifurcations of homoclinic solutions as defects are varied (see sections 4.1 and 4.2). For an isolated defect we highlight, near critical defect values, bifurcations of new homoclinic solutions (having no counterpart in the homogeneous system) or the disappearance of homoclinic solutions existing in the homogeneous case. Our method is also generalized to a finite number of defects, with the counterpart that (69) is modified by suitable higher order terms depending on the defect sequence (this procedure only provides approximate solutions of (69)). However this does not constitute a strong limitation since the full reduced equation (20) is itself a higher order perturbation of (69). Note that equation (69) is valid (according to theorem 3) for small defect sizes and µ≈0, where the parameter µdetermines for µ < 0 the (weak) degree of hyperbolicity of the fixed point βn= 0 in the homogeneous case. Our analysis does not impose conditions on the relative sizes of these parameters. In order to obtain exact breather solutions of (6) via theorem 3, it would be necessary to proceed in two steps. The first step is the one described above, where exact or approximate homoclinic solutions are obtained for the truncated problem (69). The second step is to show that these solutions persist for the complete equation 33 (20) as higher order terms are added. For this purpose a typical procedure would be to solve (20) using the contraction mapping theorem in the neighbourhood of an exact or approximate solution of (69). However, in section 4.1.4 we analyze tangent bifurcations of homoclinic orbits for which the persistence problem would become extremely complex, since it requires asymptotical techniques beyond all algebraic orders (see section 4.1.4 for more details). We shall not analytically examine the persistence of such bifurcations for the complete equation (20). Instead we shall later compare approximate solutions yn(t)≈βncos t(deduced from (69) and theorem 3) to numerically computed solutions of the original problem (6) (see section 5). This will allow us to study the validity of approximation (69) far from the small amplitude limit and as inhomogeneities become larger. The persistence of homoclinic solutions for the complete equation (20) will be shown for particular homoclinic orbits which appear through a pitchfork bifurcation at the origin when ωreaches a critical value. These orbits correspond to discrete breathers with maximal amplitude at the impurity site (see theorem 8 in section 4.1.5). This part involves standard bifurcation techniques, in contrast with the persistence problem of the above mentioned tangent bifurcations. Note that other interesting bifurcations can exist when impurities act at a purely nonlinear level (see [SKRC01, KKK03] for some examples in spatially discrete or continuous systems). This would correspond to the situation when the on-site potential in (2) has an harmonic part independent on n, whereas higher order terms are inhomogeneous. The subsequent analysis of the reduced recurrence relation would be quite different, and in particular the method developed in sections 4.1 and 4.2 (based on a linear deformation of the unstable manifold) would not apply. 4.1 Case of a single mass defect We start with the simplest case when the coefficients of (69) are constant, except at n= 0 where their value changes. To fix the idea we consider the case of a single mass defect in equation (5), i.e. Dn=d,Kn=k,An=a,Mn= 1 + mn,mn=m0δn0. The case when all lattice parameters are allowed to vary over a finite number of sites will be considered in section 4.2. For equation (6) the above assumption yields ηn=γn=κn= 0 and ²n=m0δn0. Equation (69) reads (recall ω2= Ω2+µ) βn+1 −2βn+βn−1=−(ω2m0δn0+µ)βn+B β3 n.(70) Setting βn−1=αnand Un= (αn, βn)T, equation (70) can be rewritten Un+1 =Gω(Un)−ω2m0δn0µ0 βn¶(71) 34 where Gω(Un) = µβn −αn+ 2βn+ (Ω2−ω2)βn+Bβ3 n¶.(72) One has in particular U1=A(ω, m0)Gω(U0) (73) where the linear transformation A(ω, m0) = µ1 0 −ω2m01¶(74) corresponds to a linear shear. Note that the axis α= 0 consists of fixed points of A(ω, m0). It is worthwhile to notice that the map Gωis reversible under the symmetry R: (α, β)7→ (β, α), i.e. Gω◦R=R G−1 ω. In other words, if Unis a solution of (71) for m0= 0 then R U−nis also solution. This property is due to the fact that equation (70) with m0= 0 has the invariances n→n+ 1 and n→ −n. Obviously the latter invariance still exists for m06= 0. Consequently, for all m0∈R, if Unis a solution of (71) then R U−n+1 is also solution. Now we shall use a geometrical argument to find homoclinic orbits to 0 for equation (70). In the sequel we consider the stable manifold Ws(0) of the fixed point (α, β) = 0 of Gω, and its unstable manifold Wu(0), both existing for ω < Ω. The following result follows immediately. Lemma 4 For 0< ω < Ω, equation (70) possesses an homoclinic orbit to 0if and only if Ws(0) and A(ω, m0)Wu(0) intersect. In addition it is useful to notice that Ws(0) and Wu(0) are exchanged by the reversibility symmetry R. 4.1.1 Linear case As a simple illustration, consider the linear case when Vis harmonic, in which B= 0. Equation (70) reads βn+1 −2βn+βn−1=−(ω2m0δn0+µ)βn.(75) In that case, Ws(0) and Wu(0) correspond respectively to the stable and the unstable eigenspace of a linear mapping in R2. The situation is sketched in figure 4 below. The corresponding stable eigenvalue σ∈(0,1) is given by σ= 1 −µ 2−1 2(µ2−4µ)1/2, µ =ω2−Ω2<0,(76) 35 and Ws(0) is the line β=σα,Wu(0) being the line β=σ−1α. For fixed ω < Ω, A(ω, m0) maps the unstable eigenspace on the stable eigenspace if and only if m0>0 (mass is increased at the defect) and m0=ml(ω),(77) where ml(ω) = 1 ω2(σ−1−σ).(78) 0 0.2 0.4 0.6 0.8 1 1.2 0 0.2 0.4 0.6 0.8 1 A(ω,m0)Wu(0) m0<ml A(ω,ml)Wu(0) A(ω,m0)Wu(0) m0>ml Wu(0) Ws(0) αn βn 0 0.5 1 1.5 2 0 0.5 1 1.5 2 Gω(U0) U1 U0 U2 U−1 Wu(0) Ws(0) A(ω,ml) αn βn Figure 4: Linear case (B= 0). (Left panel) Stable manifold (in the half plane α > 0) and images of the unstable manifold by A(ω, m0) for m0= 0.005, m0=ml and m0=ml+ 0.005. (Right panel) Homoclinic orbit to 0 for m0=ml. In both panels we have fixed Ω = 10, ω= 9.99, which implies ml= 0.0092. Now keeping fixed m0>0, condition (77) can be rewritten ω=ωl(m0), where (for ω < Ω) ω2 l=1 1−m2 0 [ Ω2+ 2 −(4 + m2 0Ω2(Ω2+ 4))1/2], m06= 1,(79) ω2 l=1 2(Ω2+ 2) −2 Ω2+ 2, m0= 1. The solutions of (75) homoclinic to 0 are spanned by βn=σ|n|, and the corresponding solutions of (6) in the linear case read yn(t) = βncos twith ω=ωl(m0). One recovers a classical result, i.e. if mass is increased at the defect then the linear localized mode frequency lies below the phonon band and its frequency is given by ωl. Now let us consider the effects of nonlinear terms. For this purpose we start with the simplest case of a hard potential, i.e. B > 0. The situation when B < 0 is far more complex and will be investigated later. 36 4.1.2 Nonlinear defect modes for hard on-site potentials If B > 0, Ws(0) and Wu(0) do not intersect (except at the origin) for 0 < ω < Ω. Indeed, one can show by induction that |βn|>|βn−1|>0 for any nontrivial orbit on Wu(0), which implies Wu(0) lies inside the sector formed by the lines α=βand α= 0. In the same way, Ws(0) lies inside the sector formed by the lines α=β and β= 0 hence it does not intersect Wu(0). The above property also implies that Wu(0) can be defined (globally) as the graph α=g(β) of an increasing function g, and the same holds true for Ws(0) = R Wu(0) on which β=g(α). For fixed ω∈(0,Ω), the local unstable manifold can be approximated by α=g(β) = σβ +b β3+O(|β|5), with b=σ2(σ2−σ−2)−1B < 0 (this coefficient can be computed by a classical identification procedure, using the fact that Wu(0) is invariant under Gω). Consequently, Ws(0) and Wu(0) have the local shape represented in figure 5. The same situation occurs in the limit ω≈Ω (one can locally approximate the map Gωup to any order in U, µ using the time-one map of an integrable flow [AP90], which allows to determine the shape of Ws(0) and Wu(0) close to Un= 0). In the case when m0≤0, the curves Ws(0) and A(ω, m0)Wu(0) do not intersect (A(ω, m0)Wu(0) remains inside the sector formed by the lines α=βand α= 0). However, Ws(0) and A(ω, m0)Wu(0) intersect if m0>0 provided m0> ml(ω),(80) which means there exists an orbit homoclinic to 0 for equation (70). This property is clear for m0≈mlwhere there exists a unique intersection point (in the half plane α > 0) close to U= 0 due to the local shape of Ws(0) and Wu(0). Moreover, we numerically find a unique intersection point for all values of m0satisfying (80). 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 −0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Ws(0) A(ω,ml)Wu(0) A(ω,m0)Wu(0) m0>ml Wu(0) Gω(U0) U0 U−1 U2U1 αn βn Figure 5: Case B > 0 and ω < Ω. Stable and unstable manifolds (in the half plane α > 0), and image of the unstable manifold by A(ω, m0) for m0=mland m0= 0.05 > ml. We have fixed Ω = 10 and ω= 9.9. 37 Condition (80) is equivalent to ωl< ω < Ω.(81) The amplitude of the homoclinic orbit is O(√ω−ωl) as m0is fixed and ω→ωl (at the limit A(ω, m0)Wu(0) and Ws(0) become tangent at the origin), and its spatial decay rate σis given by equation (76). This homoclinic orbit corresponds for equation (6) to a nonlinear defect mode, i.e. a nonlinear analogue of the above mentioned linear localized mode. This solution can be approximated by yn(t)≈ βncos tfor m0≈0 and the frequency ωvaries with amplitude contrarily to the linear case. In section 4.1.5 we show the persistence of the above homoclinic solution βn for the full normal form (20) and the existence of corresponding small amplitude solutions of (5) (see theorem 8). Alternatively, these solutions could be obtained using an infinite-dimensional version of the Lyapunov centre theorem (see section 4.1.5 for more details). Lastly, let us notice that the above homoclinic orbit possesses the symmetry β−n=βn, or equivalently R U−n+1 =Un. It suffices to check the latter relation for n= 0 to prove it for any n, since both solutions R U−n+1 and Uncoincide if they satisfy the same initial contidion. Since U0lies on the unstable manifold we have α0=g(β0), and in the same way β1=g(α1) since U1lies on the stable manifold. Since by definition α1=β0, this implies α0=β1and thus R U1=U0. Using the properties U1=Gω(U0)−ω2m0(0, β0)Tand R U1=U0we also deduce the relations 2α0= [2 + Ω2−ω2(m0+ 1)]β0+B β3 0(82) 2β1= [2 + Ω2−ω2(m0+ 1)]α1+B α3 1,(83) which are useful in particular for the numerical computation of U0, U1. 4.1.3 Nonlinear defect mode with algebraic decay In the situation of section 4.1.2 (B > 0), the case when m0is fixed and ω→Ω deserves a special attention. Indeed, the homoclinic orbit (αn, βn) converges in this limit towards a solution having an algebraic decay as n→ ±∞. More precisely, if ω= Ω and m0>0 then Ws(0) and A(ω, m0)Wu(0) intersect at a unique point (α1, β1) in the half plane α > 0 (see figure 6). This can be checked analytically for m0≈0 and (α, β)≈0, since the unstable manifold can be locally parametrized by α=β−(B/2)1/2β2+O(|β|3) in the half plane α > 0 (this expansion follows from a classical identification procedure). Using this relation for (α0, β0) in conjunction with (82) we find as m0→0 β0= Ω2(2B)−1/2m0+O(m2 0).(84) 38 Note that in this non-hyperbolic case, the function ghaving the unstable manifold as its graph is not C2at β= 0 (in the half plane α < 0, one has α=β+(B/2)1/2β2+ O(|β|3) on the local unstable manifold). Far from the small amplitude limit, we have also checked numerically the existence and uniqueness of Ws(0) ∩A(ω, m0)Wu(0) in the half plane α > 0. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Ws(0) A(ω,m0)Wu(0) Wu(0) Gω(U0) U0 U−1 U2 U1 αn βn Figure 6: Case B > 0 and ω= Ω. Stable and unstable manifolds (in the half plane α > 0), and image of the unstable manifold by A(ω, m0) for m0= 0.05. We have fixed Ω = 10 in this example. Consequently, there exists a solution of (70) homoclinic to 0 for ω= Ω and m0>0. This solution has an algebraic decay due to the fact that the origin is not any more an hyperbolic fixed point for ω= Ω. One can approximate the solution profile for m0≈0, using the fact that (70) admits at both sides of n= 0 a continuum limit. Indeed, setting βn≈m0β(x), x =m0n, (85) one obtains the following differential equation d2β dx2=Bβ3, x ∈(−∞,0) or (0,+∞), from which we deduce (multiply by β0and integrate) dβ dx =−sign(x) (B/2)1/2β2(86) since β(x)→0 as x→ ±∞. Using (86) and (84) one obtains the following approximation of the homoclinic solution for m0≈0 βn≈m0r2 B(m0|n|+2 Ω2)−1.(87) This yields an approximate solution yn(t)≈βncos tof (6), corresponding to a breather with an algebraic decay and a frequency ω= Ω at the bottom of the phonon band. 39 4.1.4 Case of soft on-site potentials In this section we make some considerations on the case B < 0 (soft on-site potential V) which is far more complex. We fix the parameter ωin equation (70) and let m0 vary. Using a geometrical argument, we show that two (asymmetric) homoclinic solutions of (70) having one hump near the defect site n= 0 disappear through a tangent bifurcation, at a critical value of m0which can be estimated. As we shall later numerically check (see section 5), the same features occur for the Klein-Gordon model which was locally reduced to (20), i.e. to a higher order perturbation of (70). In addition, we show (again using a geometrical argument) that a symmetric solution βnof (70) homoclinic to 0 and centered at n= 0 disappears through a pitchfork bifurcation with −βnfor m0=ml(ω). This bifurcation persists for the full normal form (20) as we shall see in section 4.1.5. In the present section we only study the simplest homoclinic bifurcations that occur in the soft potential case when m0is varied, but an infinity of tangent bifurcations occur in fact due to the complicated structure of the stable and unstable manifolds of the origin. In order to treat the case when m0≈0, we start by fixing m0= 0 and proceed perturbatively. For m0= 0, µ < 0 and B < 0, equation (70) possesses homoclinic solutions to 0. This case has been analyzed in several references with different viewpoints and for different parameter ranges, see e.g. [QX07, HT99, ABK04, PZ01, Wei99, Jam03]. The dynamics of the map Gωis rather complex due to the fact that the stable and unstable manifolds of the origin intersect transversally in general (see figures 7 and 12). This implies the existence of an invariant Cantor set on which some iterate Gp ωis topologically conjugate to a full shift on Nsymbols [AP90], which yields a rich variety of solutions and in particular an infinity of homoclinic orbits to 0. Among these different homoclinic orbits one can point out two particular ones Ui n= (αi n, βi n)T(i= 1,2), corresponding for the Klein-Gordon chain to breather solutions with a single hump near n= 0 (site-centered or bond-centered). These solutions have been described in lemma 2 and theorem 5 in the small amplitude limit. The corresponding homoclinics U1 n, U2 nare reversible, i.e. they satisfy R U2 −n+1 =U2 n (β2 −n=β2 n) and R U1 −n+2 =U1 n(β1 −n+1 =β1 n). In figure 7, the point with label 2 lying on the axis α=βcorresponds to U1 1, and the points with labels 3,1 correspond to U2 0, U2 1respectively. Obviously any translation of Ui ngenerates a breather solution of (6) having its maximal amplitude near a different site. Now let us consider the situation when ωis kept fixed and a small mass defect m0is introduced in (70). As illustrated in figure 7, each of the above solutions is structurally stable. For example, let us consider in figure 7 the intersection points 1,2,3 between Wu(0) and Ws(0). Each of these intersections persists (points 10,20,30in figure 7) as the linear shear A(ω, m0) is applied to Wu(0) for m0≈0 (dashed line in figure 40 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 1’ 1 2 2’ 3’ 3 αn βn Figure 7: First intersection points between the stable and unstable manifolds for parameters ω= 9.9 (µ=−1.99) and B=−75. The dashed line depicts the image of the unstable manifold by the linear shear A(ω, m0) for m0= 0.005. 7). Let us examine the corresponding solutions of (71) and the related breather solutions of the Klein-Gordon model. We denote by ˜ U2 n= (˜α2 n,˜ β2 n)Tthe solution of (71) with initial data ˜ U2 1at the point 10. This solution is homoclinic to 0 according to lemma 4. Repeating an argument of section 4.1.2, one can show that R˜ U2 −n+1 =˜ U2 n, i.e. ˜ β2 −n=˜ β2 n. Consequently, ˜ U2 ncorresponds to an (approximate) breather solution of (6) centered at the defect site n= 0. This solution is a small deformation of the site-centered breather y2 nof theorem 5. Now let us denote by ˜ U1 nthe homoclinic solution of (71) with initial data ˜ U1 1 at the point 20. It corresponds to an (approximate) breather solution of (6), whose profile is a small deformation of the breather y1 ncentered between n= 0 and n= 1 (see theorem 5). Since ˜ U1 1does not belong to the line α=β(it lies at a distance O(|m0|)), the corresponding breather solution is not symmetric any more, which was expected since the atomic masses at n= 0,1 are different. Lastly we note ˜ U3 n= (˜α3 n,˜ β3 n)Tthe homoclinic solution of (71) with initial data ˜ U3 1at the point 30. Since ˜ U3 1is O(|m0|)-close to U2 0(point with label 3), ˜ U3 nis a small deformation of the solution U2 n−1existing for m0= 0. In other words, ˜ U3 1 corresponds to a small deformation of the breather y2 n−1centered at n= 1. The mass defect at n= 0 breaks the mirror symmetry of the solution, since its amplitude has only the imperfect symmetry ˜ β3 −n+1 −˜ β3 n+1 =O(|m0|) for n6= 0. A more delicate question concerns the continuation and the possible bifurcations of the above homoclinic solutions as m0is further varied. The evolution of ˜ U1 n,˜ U2 n,˜ U3 n depends on the structure of the homoclinic windings near U1 1, U2 1, U2 0. Numerically we find that the lobes formed near these points by the stable and unstable manifolds 41 where θn= Ω2(ηn+ηnγn+γn)−ω2²n+κn,ω2= Ω2+µ. In the sequel we shall note ε=k{θ}k`∞(Z)+k{κ}k`∞(Z). Setting βn−1=αnand Un= (αn, βn)T, equation (98) can be rewritten Un+1 =Fn(Un),(99) Fn(α, β) = µβ −(1 + κn)α+ (2 + θn−µ)β+Bβ3¶.(100) Noting F=Gωfor simplicity (see definition (72)), one can observe that Fn= (I+Tn)F+O(|κn||β|3),(101) Tn=µ0 0 θn+ (µ−2) κnκn¶. Note that higher order terms are absent from equation (101) if κn= 0. Since Fn=Ffor |n| ≥ n0+ 1 one has the following property. Lemma 5 Fix µ < 0and denote by Ws(0),Wu(0) the stable and unstable manifolds of the fixed point U= 0 of F. Consider the nonlinear map G=Fn0◦Fn0−1◦ ··· ◦ F−n0◦F−2n0−1. Equation (99) possesses an homoclinic orbit to 0if and only if Ws(0) and G(Wu(0)) intersect. Lemma 5 is hard to use for analyzing homoclinic solutions since it involves a nonlinear transformation Ginstead of a linear one as in lemma 4. However one can recover the linear case when replacing Fnby a suitable approximation ˆ Fn, equal to Fnup to higher order terms. This is possible thanks to property (104) of lemma 6 below. In the sequel we note L=DF(0) = µ0 1 −1 2 −µ¶. Lemma 6 Consider the collection of maps ˆ Fn(−n0≤n≤n0) defined by ˆ Fn=AnF◦A−1 n−1,(102) where A−n0−1=Iand for n≥ −n0 An=LnLn−1··· L−n0L−n−n0−1, Ln= (I+Tn)L. (103) The map ˆ Fnis a leading order approximation of Fn, i.e. ˆ Fn=Fn+O(εk(α, β)k3). Moreover one has the property ˆ Fn0◦ˆ Fn0−1◦ ··· ◦ ˆ F−n0=A F2n0+1,(104) where A=An0reads A=Ln0Ln0−1···L−n0L−2n0−1=I+O(ε).(105) 48 Proof. First we note that the sequence Ansatisfies A−n0=I+T−n0and An+1 = (I+Tn+1)L AnL−1(106) for all n≥ −n0−1. It follows for −n0≤n≤n0 ˆ Fn= (I+Tn)L An−1L−1F◦A−1 n−1.(107) Now let us note that An=I+O(ε). Moreover, the following identity holds true for any parameter-dependent matrix M∈M2(R) with kMk=O(ε) F◦(I+M) = L(I+M)L−1F+O(εk(α, β)k3).(108) Consequently one has also F=L(I+M)L−1F◦(I+M)−1+O(εk(α, β)k3). Using this property in equation (107) leads to ˆ Fn= (I+Tn)F+O(εk(α, β)k3). Using (101) this yields ˆ Fn=Fn+O(εk(α, β)k3), therefore ˆ Fnis a leading order approximation of Fn. Property (104) follows directly from the definition of ˆ Fn. 2 It is worthwhile stressing that A=DG(0), where Gis the nonlinear transformation introduced in lemma 5. Now we fix in addition ˆ Fn=F=Fnfor |n| ≥ n0+ 1. According to lemma 6 we have also ˆ Fn=Fn+O(εk(α, β)k3) for |n| ≤ n0. In the sequel we approximate system (99) by the new one Un+1 =ˆ Fn(Un).(109) Property (104) implies the following result, since Wu(0) is invariant under F2n0+1. Lemma 7 Fix µ < 0and denote by Ws(0),Wu(0) the stable and unstable manifolds of the fixed point U= 0 of F. Equation (109) possesses a solution Unhomoclinic to 0if and only if Ws(0) and A(Wu(0)) intersect, where the matrix A=I+O(ε)is defined in lemma 6. The intersection point corresponds to Un0+1. Consequently, as in section 4.1 one recovers the problem of finding the intersection of Ws(0) with the image of Wu(0) under the (near-identity) linear transformation A. Note that A=I+T0in the single defect case n0=0. Here we shall not attempt to relate the bifurcations of breather solutions of (6) with the properties of the inhomogeneities, via an analysis of homoclinic solutions of 49 (109). This question will be considered in future works using the simplification provided by lemma 7. As for a single defect, for B < 0 one can expect multiple tangent bifurcations between (deformations of) site-centered and bond-centered breathers as inhomogeneities are varied, due to the winding structure of Wu(0) and Ws(0). It is now interesting to compute the leading order contribution of the sequence of inhomogeneities to the matrix A. This is the object of the following lemma. Lemma 8 The matrix Aof lemma 6 takes the form A=I+M+O(ε2+ε|µ|), where M=µM11 M12 M21 M22 ¶, M11 = 2n0 X n=0 n(n+ 1)ρn0−n−n κn0−n, M12 = 2n0 X n=0 −n2ρn0−n+n κn0−n, M21 = 2n0 X n=0 (n+ 1)2ρn0−n−(n+ 1) κn0−n, M22 = 2n0 X n=0 −n(n+ 1)ρn0−n+ (n+ 1) κn0−n, ρn= Ω2(ηn+γn−²n). Proof. Since Tn=O(ε) it follows from definition (105) A=I+ 2n0 X n=0 LnTn0−nL−n+O(ε2).(110) Now we use the expansions Tn=Mn+O(ε2+ε|µ|), Mn=µ0 0 ρn−κnκn¶, L=Lc+O(|µ|), Lc=µ0 1 −1 2 ¶, to obtain A=I+M+O(ε2+ε|µ|), M = 2n0 X n=0 Ln cMn0−nL−n c, 50 where Ln c=µ−n+ 1 n −n n + 1 ¶. Then simple computations lead to the coefficients of Mprovided above. 2 Interestingly, lemma 8 shows that the influence of the inhomogeneities on the set of homoclinic solutions depends at leading order (via the matrix I+M) on algebraically-weighted averages of {κ}and {ρ}. Now let us return to the original parameters mn, dn, an, kndescribing the lattice inhomogeneities (see equations (5) and (6)), with mn=dn=an= 0 for |n| ≥ n0+1, kn= 0 for n≤ −n0−1 and n≥n0. Let us note ˜ε=k{mn/m}k`∞+k{dn/d}k`∞+ k{an/a}k`∞+k{kn/k}k`∞.One obtains κn=kn−1−kn k+O(˜ε2), ρn=rn+O(˜ε2), where rn= Ω2(dn d+ 2an a−mn m) is a linear combination of the on-site potential and mass defect impurities. Some coefficients of Mcan be simplified since 2n0 X n=0 κn0−n=O(˜ε2), 2n0 X n=0 −n κn0−n=1 k n0−1 X n=−n0 kn+O(˜ε2). Noting Ik= n0 X n=−n0 nkrn, J0=1 k n0−1 X n=−n0 kn, one finally obtains A=I+˜ M+O(˜ε2+ ˜ε|µ|) with ˜ M=µ˜ M11 −˜ M11 +n0I0−I1 ˜ M11 + (n0+ 1)I0−I1−˜ M11 ¶ and ˜ M11 =n0(n0+ 1)I0−(2n0+ 1)I1+I2+J0. Consequently, the matrix Adepends (at leading order in ˜εand µ) on the average values I0, J0of rn,kn/k, and on the weighted averages I1, I2of rn(with linear and quadratic weights respectively). Since Tr( ˜ M) = 0, it follows Tr(A) = 2+O(˜ε2+˜ε|µ|) and Det(A) = 1 + O(˜ε2) (Det(A) is independent of µdue to identity (105)). As a consequence, in order to study the spectrum of Afor ˜ε, µ ≈0 (and determine to which type of linear transformation it corresponds) it would be necessary to compute the quadratic terms in (˜ε, µ) in its expansion. 51 5 Numerical results We have performed numerical computations in order to check the range of validity of the analysis of section 4.1, and in particular if discrepancies appear for large amplitude solutions or if parameters (m0, ω) are moved away from (0,Ω). More precisely, we have computed breather solutions of the Klein-Gordon lattice ω2(1 + mn)d2yn dt2+ Ω2V0(yn) = yn+1 −2yn+yn+1 (111) with a single mass defect mn=m0δn,0and periodic boundary conditions y−N(t) = yN(t). In general we have used a lattice with 101 particles, except for the computations of breathers with algebraic decay (case ω= Ω) where 401 particles have been considered. The computations have been compared with homoclinic orbits to 0 of the two-dimensional map (71). For the numerical computations we have always fixed Ω = 10 (recall Ω is the lower phonon band edge for the infinite system). This can be done taking, for instance, k= 0.01, and d= 1 in the original problem (6). For the potential Vwe have chosen a polynomial of degree 4 with V00 (0) = 1. 5.1 Hard potentials To start we have considered the simplest case of a hard potential, i.e. a potential with a strictly positive hardening coefficient B(see definition (22)). We have chosen V(x) = x2 2+x4 4,(112) for which B= 75. In this case, the reduced map (71) possesses a unique orbit homoclinic to 0 in the sector α > 0, β > 0, for m0>0 and ωl< ω < Ω. An example of this homoclinic orbit is shown in figure 5 for a frequency ω= 9.9 (µ=−1.99) and a mass defect m0= 0.05. In figure 8 (left panel) we compare the approximate solution yn=βncos tobtained with this homoclinic orbit (circles) with the exact breather profile computed with the standard numerical method based on the anti-continuous limit [MA96] (continuous line). The agreement is excellent even if the solution profile is very localized. Indeed, as one computes the eigenvalues σ, σ−1(equation (76)) of the linearized map (75) with m0= 0, one obtains σ≈0.27, which implies a strong spatial localization visible in figure 8. The accuracy of the centre manifold reduction (a priori expected for σ≈1) is surprisingly good in this parameter regime. The breather solution can be continued for decreasing frequencies up to ωl≈ 9.8369, which is the frequency of the linear defect mode at which the breather solution bifurcates. Figure 8 (right panel) compares again the numerically computed breather profile and the approximate solution obtained with the homoclinic orbit, 52 −10 −5 0 5 10 0 0.05 0.1 n yn(0) −10 −5 0 5 10 0 1 2 3 4 5 6x 10−3 n yn(0) Figure 8: Comparison between the profile of a breather solution (continuous line) of the Klein-Gordon system (111) with hard potential (112) and the approximate solution yn=βncos t(circles) constructed with the homoclinic orbit of (71). We have considered a mass defect m0= 0.05. In the left panel we have chosen a frequency ω= 9.9 (µ=−1.99). In the right panel we have fixed ω= 9.837 (µ=−3.23) very close to ωl(note the change of scale for the vertical axis). but now very close to this bifurcation point (at ω= 9.837, i.e. µ=−3.23). We still observe an excellent agreement. Note that the oscillations amplitudes are very small, but the solution is still strongly localized. For increasing frequencies the continuation path ends up at the lower edge of the phonon band ω= Ω (µ= 0). For this particular frequency value the breather solution (see continuous line in figure 9, left panel) presents an algebraic decay which is very well described by approximation (87). This approximation fails to describe the maximum amplitude of the oscillation β0for these parameter values. This is not surprising since β0is not small, and βnvaries rapidly near n= 0, hence m0should be further decreased to attain the domain of validity of the ansatz (85) near the solution centre. However, the value of β0obtained from the exact homoclinic orbit of (71) fits very well the maximum amplitude of the breather solution, as it is shown in figure 9, right panel. Note that the agreement is very good even for very large amplitudes or very large mass defect i.e. within a surprisingly large parameter range for a local theory. It is interesting to remark that the accuracy of this fit depends on the symmetry of the potential V(x) we have chosen. Figure 10 shows what happens if we add to the polynomial potential (112) a cubic term x3/6 that breaks its symmetry. The range of validity of our leading order approximation reduces significantly. A similar result was obtained in [SRJCA04] for breather solutions in spatially homogeneous 53 −25 −20 −15 −10 −5 0 5 10 15 20 25 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 n yn(0) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 m0 y0(0),y1(0) Figure 9: Left panel : breather solution at the lower edge of the phonon band ω= Ω (µ= 0) for a mass defect m0= 0.05 and the symmetric potential V(x) = x2/2+x4/4. The continuous line corresponds to the numerically computed breather solution. The circles represent approximation (87) of the homoclinic orbit that fits very well the algebraic decay of the breather tails. Right panel : the continuous lines represent the amplitude of the breather solution at n= 0 and n= 1 versus mass defect. The circles correspond to the homoclinic solutions of the nonlinear map (71) (the upper plot represents β0and the lower plot β1). Fermi-Pasta-Ulam lattices. Obviously the agreement would be improved by taking into account the Taylor expansion of the reduction function φ(see theorem 3) and computing the normal form at a higher order. 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 m0 y0(0),y1(0) Figure 10: Same computation as in figure 9, right panel, but now for the asymmetric potential V(x) = x2/2 + x3/6 + x4/4. Finally we have numerically studied the spectral stability of the breather solu- 54 tions by finding the eigenvalues of the Floquet operator, which gives us the evolution of any small perturbation over one period [Aub97]. We have checked that all breather solutions in the gap ωl< ω < Ω are spectrally stable, at least for the value of the frequency parameter Ω = 10 we have considered. 5.2 Soft potentials In the case of soft potentials (when the coefficient Bdefined by (22) is strictly negative), the situation is far more complex due to the much more intricate structure of the intersections between the stable and unstable manifolds. Therefore one expects a richer bifurcation scenario as parameters (breather frequency, mass defect) are varied. Our computations have been performed with the symmetric potential V(x) = x2 2−x4 4,(113) for which B=−75. Let us recall some basic features of the analysis performed in section 4.1.4, in order to compare the results with numerical computations. For the (truncated) reduced mapping (71) with m0= 0, figure 7 shows some intersections of stable and unstable manifolds emanating from the saddle point at the origin, for a frequency value ω= 9.9<Ω. Iterating the map with an initial condition U1at the homoclinic point with label 1, we obtain an homoclinic orbit which corresponds to a one-site breather centered at n= 0. With an initial condition U1at the homoclinic point with label 2, the corresponding breather is a two-site breather with maximal amplitude at n= 0 and n= 1. An initial condition U1at the homoclinic point with label 3 (symmetric of point 1 respect to the line α=β) corresponds to a one-site breather centered at site n= 1. The dashed line of figure 7 depicts the image of the unstable manifold by the linear shear A(ω, m0) for m0= 0.005. As m0increases A(ω, m0)Wu(0) moves further down so that the intersection points 2’ and 3’, corresponding to homoclinic orbits of the inhomogeneous problem, get closer and closer. So there exists a critical value of m0for which these intersection points collide and then disappear. In fact we have checked numerically that this tangent bifurcation occurs at a critical value m0∈(0.00963,0.00964) for problem (71). This critical value can be approximated using equations (91)-(92), which yields m0≈0.009632 in the present case. These results correspond very precisely to a breather bifurcation numerically observed in the Klein-Gordon chain (at a critical value m0∈(0.00963,0.00964)) and depicted in figure 11. The upper branch of figure 11(a) represents the energy of a breather solution corresponding to point 2’. For m0≈0, the breather has a maximal amplitude at sites n= 0,1. A profile of this breather for m0= 0.0093, close to the bifurcation 55 point, is shown in figure 11(b), where the amplitude is now much larger at n= 1. The lower branch of figure 11(a) represents a one-site breather centered at n= 1 and corresponds to point 3’. Its profile for m0= 0.0093 is shown in Figure 11(c). We have numerically computed the Floquet spectra of these breather solutions for the parameter values of figure 11. The solutions on the lower branch are spectrally stable, whereas the solutions on the upper branch are unstable. As in section 5.1, we have also computed one-site breathers centered at the mass defect, corresponding to point 1’ in figure 7. Again we have found an excellent agreement between the numerically computed breather profiles and the approximate solutions obtained using the map (71). As expected from the analysis of section 4.1.4, these breathers survive up to m0=ml(ω), i.e. up to a much higher value of m0 than the families 2’, 3’ described above. 0 0.002 0.004 0.006 0.008 0.01 0.025 0.03 0.035 0.04 m0 Energy (a) −25 −20 −15 −10 −5 0 5 10 15 20 25 0 0.05 0.1 0.15 0.2 0.25 n yn(0) (b) Upper branch −25 −20 −15 −10 −5 0 5 10 15 20 25 0 0.05 0.1 0.15 0.2 0.25 n yn(0) (c) Lower branch Figure 11: Tangent bifurcation between breather solutions numerically computed in a Klein-Gordon chain with a soft potential. The chain presents a mass defect m0at n= 0, and the bifurcation occurs as m0is increased. In the left panel, the breathers energies E=Pn∈ZΩ2V(yn(0)) + (yn+1(0) −yn(0))2/2 are depicted versus m0(the breathers are even in twith frequency ω= 9.9). For m0≈0, the upper branch represents a two-site breather centered between sites n= 0 and n= 1. The lower branch represents a one-site breather centered at n= 1. The breathers profiles close to the bifurcation point are plotted in the right panels (the value of m0is marked with a dashed line in the left panel). A part of the intersecting stable and unstable manifolds is shown in the left panel of figure 12. Due to their complicated windings, new intersection points appear between A(ω, m0)Wu(0) and Ws(0) as m0is chosen in certain windows of the parameter space, giving rise to new homoclinic solutions of (70). An example is shown in the region marked with a rectangle (see the details in the right panel of figure 12). For some value of m0∈(0.01064,0.01065), a new intersection point between A(ω, m0)Wu(0) and Ws(0) appears. As m0is 56 further increased, this inverse tangent bifurcation gives rise to two new homoclinic points 5’ and 6’. Correspondingly, we have numerically checked that an inverse tangent bifurcation occurs in the Klein-Gordon chain at a critical value of m0∈ (0.01064,0.01065), giving rise to new breather solutions which do not exist in the homogeneous chain. The point 4’ in figure 12 also exists for m0= 0. Returning to figure 7, it is obtained by applying the inverse map G−1 ωto the point with label 2. In the homogeneous limit m0= 0, this homoclinic point corresponds consequently to a two-site breather centered between n= 1 and n= 2. As figure 12 shows, an increase of the mass defect m0moves point 5’ against point 4’ until they collide and disappear through a new tangent bifurcation. This tangent bifurcation is also numerically found in the Klein Gordon chain at critical value of the mass defect very close to the theoretical one (in both cases one obtains m0≈0.01268). Figure 13 shows the bifurcation diagram of the numerically computed breathers corresponding to homoclinic points 4’, 5’, 6’ (left panel), and gives their profiles for a given value of m0in the right panel. A numerical computation of Floquet spectra shows that all these breathers are unstable. −0.5 0 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 αn βn 0.06 0.08 0.1 0.12 0.14 0.16 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 6’ 5’ 4’ αn βn Figure 12: Emergence of new intersection points between A(ω, m0)Wu(dashed curve) and Ws(drawn with a full line) as the mass defect is increased. The figure corresponds to m0= 0.012 and ω= 9.9. The right panel shows a zoom of the left panel over the region marked with a rectangle. The new homoclinic points 5’ and 6’ correspond to new breather solutions of the Klein-Gordon lattice. The point with label 4’ corresponds to a two-site breather, which exists in the homogeneous lattice and persists for m0≤0.012. As a conclusion, we have seen that the truncated normal form (70) allows one to predict with a high precision certain breather bifurcations in the Klein-Gordon chain, which occur as the mass defect m0is varied. These bifurcations depend on 57 [SM97] J.A. Sepulchre and R.S. MacKay. Localized oscillations in conservative and dissipative networks of weakly coupled autonomous oscillators, Nonlinearity 10 (1997), 679-713. [SM98] J.A. Sepulchre and R.S. MacKay. Discrete breathers in disordered media, Physica D 113 (1998), 342-345. [ST88] A.J. Sievers and S. Takeno. Intrinsic localized modes in anharmonic crystals, Phys. Rev. Lett. 61 (1988), 970-973. [Sir05] Y. Sire. Travelling breathers in Klein-Gordon lattices as homoclinic orbits to p-tori, J. Dyn. Diff. Eqs. 17 (2005), 779-823. [SKRC01] A.A. Sukhorukov, Yu.S. Kivshar, J.J. Rasmussen and P.L. Christiansen. Nonlinearity and disorder : classification and stability of nonlinear impurity modes, Phys. Rev. E 63 (2001), 036601. [Swa99] B.I. Swanson et al. Observation of intrinsically localized modes in a discrete low-dimensional material, Phys. Rev. Lett. 82 (1999), 3288-3291. [TP96] J.J.L. Ting and M. Peyrard. Effective breather trapping mechanism for DNA transcription, Phys Rev E 53, n.1 (1996), 1011-1020. [Van89] A. Vanderbauwhede. Centre manifolds, normal forms and elementary bifurcations, Dynamics Reported 2 (U. Kirchgraber and H.O. Walther, eds), John Wiley and Sons Ltd and B.G. Teubner (1989), 89-169. [vECLP06] Titus S. van Erp, Santiago Cuesta-Lopez, and Michel Peyrard. Bubbles and denaturation in DNA, Eur. Phys. J. E 20 (2006), 421-434. [VMZ03] L. V´azquez, R.S. MacKay, M.P. Zorzano (Eds). Localization and Energy Transfer in Nonlinear Systems, Proceedings of the Third Conference (San Lorenzo de El Escorial, Spain 17 - 21 June 2002), World Scientific (2003). [Wei99] M.I. Weinstein. Excitation thresholds for nonlinear localized modes on lattices, Nonlinearity 12 (1999), 673-691. 64