Breathers as Metastable States for the Discrete NLS equation

We study metastable motions in weakly damped Hamiltonian systems. These are believed to inhibit the transport of energy through Hamiltonian, or nearly Hamiltonian, systems with many degrees of freedom. We investigate this question in a very simple model in which the breather solutions that are thought to be responsible for the metastable states can be computed perturbatively to an arbitrary order. Then, using a modulation hypothesis, we derive estimates for the rate at which the system drifts along this manifold of periodic orbits and verify the optimality of our estimates numerically.

in a system of N rotators, coupled to their nearest neighbors, there are regions of phase space in which the transport is very slow and they give theoretical and numerical estimates indicating that the slow transport is due to the metastability of periodic breathers in the system.We note, that breathers and the transport of energy in Hamiltonian lattices, have also been among the topics of R. de la Llave's research ( [3], [4]).
Here, we investigate this question further in a very simple model in which the breathers can be computed perturbatively to an arbitrary order.Then, using a modulation hypothesis see Eq.( 13) we derive estimates for the rate at which the system drifts along this manifold of periodic orbits and verify the optimality of our estimates numerically.
We now recall in more detail the situation considered in [2].The authors start from a Hamiltonian system of N rotators coupled to their nearest neighbors.They choose initial conditions in which almost all of the energy is the system is concentrated on one of the rotators at the end of the system, and they impose a very weak damping at the opposite end.Since it is easy to prove that the energy eventually tends to zero, it is forced to flow from the end of the system where it is initially located to the end where it dissipates.The authors find that the dissipation rate is very slow and their numerics seem to indicate that for a very long period, the solution remains near an approximately periodic, spatially localized, state reminiscent of the discrete breathers often found in Hamiltonian lattices [5], [6], [7].
In the present work we look again at this problem in a finite, discrete nonlinear Schrödinger equation: Here (∆u) j = u j−1 − 2u j + u j+1 , with obvious modifications for j = 1 or N .We will focus primarily on the case N = 3, for simplicity, but in principle, our methods apply to systems with arbitrarily many degrees of freedom, and we plan to return to the consideration of the general case in a future work.We will choose initial conditions for this system in which essentially all of the energy is in mode u 1 , and will add a weak dissipative term to the last mode as in [2] by adding to Eq.(1) a term of the form iγδ N,j u j , i.e., we add dissipation to position N , at the other end from the energetic mode.Eventually, this will lead to the energy of the system tending to zero, but we are interested in what happens on intermediate time scales.
If our initial conditions are chosen so that u 1 (0) = L, and all other u j (0) = 0, then we expect that at least initially, the coupling terms between the various modes will play only a small role in the evolution and the system will be largely dominated by the equation for u 1 : with solution u 1 (τ ) = e iL 2 τ L-i.e., we have a very fast rotation with large amplitude.With this in mind, we introduce a rescaled dependent variable and rewrite the equation in a rotating coordinate frame by setting: Then wj satisfies We now add dissipation by adding a term which acts on the last variable, with γ ≥ 0, Rearranging, and dividing by L 2 gives Finally, we define ε = L −2 , and rescale time so that τ = εt.Setting w(t) = w(τ ), we arrive finally at We will study Eq.( 3) for the remainder of this paper.We will also sometimes rewrite this equation in the equivalent real form by defining w j = p j + iq j , which yields the system of equations, for j = 1, . . ., N : Note that if γ = 0, this is a Hamiltonian system with: (5) 2. The existence of breathers.Based on the rescaling and rewriting the equation in a rotating frame, we might expect that when γ = 0, that is, when we have no damping, Eq.( 3) has a fixed point, close to w 1 = 1, and we will prove that this is the case below.Furthermore, this fixed point is spatially localized-i.e., its components w j decay rapidly with j.Thus, it is an example of the breather solutions often studied in Hamiltonian lattices.Breathers often come in families and that is the case here too-we will prove that when γ = 0 Eq.( 3) has a family of periodic orbits with frequency close to zero which are also spatially localized.
Finding a periodic solution of the form Eq.(2) (i.e., a fixed point in the rotating coordinate system) is reduced to finding roots of the system of equations One can easily compute, for fixed N (not too large), that there is a solution with q * j = 0, and p * j given as a series in ε.For example, for N = 3 the series takes the form: Remark 1.Since Eq.( 3) is invariant under complex rotations w j → e iϑ w j , we actually have a circle of fixed points (when γ = 0).However, these are the only fixed points with |w 1 | ≈ 1.
Remark 2. Since the fixed points in Theorem 2.1 are real valued, they correspond to a fixed point of the system of equations Eq.( 4) of the form (p * , q * = 0).
Numerics, for example for the case of ε = 0.01 and N = 4, shows that for initial conditions close to these states the periodic solutions emerge naturally, and after a transient, the scaling of the different components indeed matches those predicted by the theory, i.e., each oscillator has energy about ε 2 compared to the previous one, see Fig. 1.

Figure 1.
Numerical illustration of the dynamics with (weak) dissipation.The parameters are N = 4, γ = 0.2, and ε = 0.01.Shown is the "energy" of the degree of freedom i, i = 1, . . ., 4. Note that the transients vanish after some time and then the energies settle at about ε −2(i−1) .Here we define them as p 2 i + q 2 i .Note also that the dissipation is so slow that no decrease can be observed in the graph of p 2 1 + q 2 1 over the time scale considered.
Proof.If we insert w j (t; ϕ) = e itϕ p * j (ϕ) into Eq.(3), and take real and imaginary parts, we find that the amplitudes p * of these periodic orbits are solutions of Setting p 0 j = δ j,1 , we have F j (p 0 ; 0, 0) = 0 , for all j.Furthermore, the Jacobian matrix at this point is the diagonal matrix which is obviously invertible.Thus by the implicit function theorem, for (ε, ϕ) in some neighborhood of the origin, Eq.( 8) has a unique fixed point p = p * (ε, ϕ), and since F depends analytically on (ε, ϕ), so does p * (ε, ϕ).
As in the case of the fixed points, the invariance under complex rotations generates a two parameter family of periodic orbits Remark 3. Solutions similar to those considered here have been considered by many authors, including for the case of infinity many nodes.See for example [8].
3. A modulation approach to describe the motion along the periodic family.In this section we derive equations that describe the motion of the system near the family of periodic orbits we found above when the dissipation is included in the equations of motion.We focus specifically on the cases N = 2, 3, and 4 sites in this section.We show that these orbits form a "metastable manifold" in the sense of [9] in that once a trajectory comes close this set of states it remains near them for a very long time, moving along the manifold of periodic solutions at a very slow rate, see Fig. 2. We calculate both the rate at which one drifts along this family, and the dissipation rate that this induces in the system and show that both of these predictions are in accord with our numerics.3) slides along the cylinder.It actually does not converge to it but will stay at some small, finite, distance from it.So the cylinder is Lyapunov stable in the sense of [10], at least as long as ε stays small.
Consider the linearization of the system Eq.( 4) around the periodic orbit w * j (t; ϕ) = e itϕ p * j (ϕ).We continue with N = 3. Setting where the p * j are the solutions of Eq.( 6).Similarly, Similar expressions hold for other values of N .
Remark 4. By the covariance of the equations, similar formulas hold for ϕ ∼ 0 and arbitrary ϑ, but we continue with ϕ = 0 and ϑ = 0.
We next argue that the spectrum of M consists of a double eigenvalue 0, and 2N − 2 imaginary eigenvalues, and we have verified this for N = 2, . . ., 6.
The discussion starts with ε = 0.Then, we have the "hatted" quantities We set and study first the spectrum of M .The spectrum of M will then be close to that of M , and changing ϑ leaves the spectrum unchanged.
This comes from the invariance of Eq.( 1) under u → e iϑ u.That is to say, we use the fact that e iϑ p * j is a fixed point of the system of equations for all ϑ, differentiate with respect to ϑ, and set ϑ = 0. Taking real and imaginary parts of the resulting equation shows that v (1) is an eigenvector with eigenvalue zero.The analog of e (2) is then found as follows: From the form of M we see that 1) .
Note that B −1 exists when ε is small enough.Therefore, M has an eigenvalue 0 and is of Jordan normal form, in the subspace spanned by v (1) , v (2) .The vector v (2) is mapped by M onto v (1) , which means that v (1) and v (2) span the zero eigenspace.As with v (1) , v (2) comes from the fact that we have a family of fixed points/period orbits, this time depending on ϕ.As we showed in Theorem 2.2, there is a family of fixed points of Eq.( 8) depending on ϕ.Differentiating this equation with respect to ϕ and setting ϕ = 0 leads to the form of v (2) .The other eigenvalues are (of course) purely imaginary, and perturbation theory (for N = 3) shows that they are all simple for small enough |ε|.If one studies the eigenvalues of M 2 , one finds that they are By analyticity, there is an open set of |ε| where this is negative.And the eigenvalues of M are ± the square roots of the λ k , and so they are purely imaginary.
We next study what happens when one adds dissipation on the coordinate N .In the standard basis, when N = 3, the dissipation is just given by −γε times a matrix C: In the full space, we have the 2N × 2N matrix Assume now that the eigenvalues are all simple.Then, for the nonvanishing ones, one expects that adding dissipation of the form qN = −γεq N and ṗN = −γεp N moves these eigenvalues into the left half plane, by an amount proportional to γ.We illustrate this, numerically in Fig. 3, for N = 3 with the two top lines.While we do not have a proof, numerical studies show clearly that the spectrum, for N = 3, 4, when γ > 0, has the following properties: 1. Eigenvalues which are non-zero when γ = 0 acquire negative real parts proportional to γ.We have verified this for N = 2, 3, 4 numerically.2. The zero eigenvalues which appear when γ = 0 remain very small when γ is non-zero.We address these eigenvalues using modulation theory, below.
We illustrate these observations in Fig. 3.The real parts are shown in Fig. 3, where one sees that one (pair of) eigenvalues has a positive real part, which is very small, while the others have negative real parts.The imaginary part of the first eigenvalue is about 1.3 • 10 −5 when N = 3 while the other eigenvalues have imaginary parts a few percent less than 1 (not shown).
Our numerics indicate that as long as γ is small, solutions remain near the family of periodic orbits for very long times, even though no true periodic orbits survive this perturbation.We note that when γ = 0, p * (ϕ) is no longer a fixed point of Eq.( 4), so the relevance of the spectrum of M − γε C for the dynamics of Eq.( 4) is not immediately apparent, but our approach is similar to that of [11] where as long as the "drift" along the family of periodic orbits is slow, we can "freeze" the coefficients in the linearized vector field to understand the dynamics over some finite interval.
To consider the evolution of solutions of Eq.( 3) (or equivalently, Eq.( 4)) for initial conditions near the family of periodic orbits, we write the solution as To emphasize, if γ = 0, we know that we have solutions of this form with ϕ and ϑ constant and W j (t) ≡ 0 for all t.Motivated by the modulation theory approach to understand perturbations of coherent structures (see, e.g.[12], or [11]), we split the various evolution equations.For this, we consider the adjoint eigenvectors of J∇H x=X * (ϑ,ϕ) , corresponding to the 0 eigenspace.Denote them n (k) (ϑ, ϕ) with k = 1, 2. They are defined by n (1) , ∂ ϕ X * = n (2) , ∂ ϑ X * = 0 .
Due to the special and simple form of M = J∇H| x=X * , we can write down the adjoint eigenvectors that span the zero subspace very simply, at least when the parameter ϑ equals 0. As we found earlier, the zero eigenspace of M is spanned by v (1) = (0, 0, 0, , and if we write out M in its block form, we found x * = B −1 p * .
Because of the block form of M and the fact that A and B are symmetric, we have A 0 But then, since we know from the computation of the eigenvectors of M that Bp * = 0 (where p * is the vector corresponding to the first three components of v (1) ), we can check immediately that Thus, ñ(1) and ñ(2) span the zero eigenspace of the adjoint operator.If we normalize these vectors so that Eq.( 14) holds, we have n (1) = ν 1 ñ(1) and n (2) = ν 2 ñ(2) where in the case of N = 3, the normalization constants can be shown to satisfy ν j = 2 + O(ε), for j = 1, 2, using the facts that p * 1 = 1+O(ε) (with all other components of p * j of higher order in ε) and that B −1 = ( B) −1 + O(ε), where B has the form Eq.( 11).Thus, For ϑ = 0, the linearization of the vectorfield J∇H| x=X * can be conjugated to the ϑ = 0 case by a rotation in the (q j , p j ) plane, so one can derive the corresponding eigenvectors by taking advantage of this rotation.
Inserting this form of the solution into Eq.(3), and expanding, we find, writing ϕ for ϕ(t), ϑ for ϑ(t) and p * j (ϕ) for the composition p * j (ϕ(t)): Recall that p * j is a periodic orbit of the γ = 0 equations, so This leads to the cancellation of many terms in Eq.( 17) and we are left with We have defined W so that it lies (at least initially) in the rapidly damped subspace of the linearized equations, and hence we expect it to decay quickly to a small value and remain small.We plan to return to this point in future research (using the methods of [11]) but for the moment we set W ≡ 0 and focus on how the parameters ϕ(t), ϑ(t) evolve-i.e., how we move along the family of periodic orbits, assuming that we remain very close to it.If we set W ≡ 0, in Eq.( 18) we are left with (t φ(t To isolate the equations for θ and φ we separate the real and imaginary parts of this equation, and rewrite it in vector form as X = (p, q), where as above, q represents the imaginary part of the complex vector, and p its real part.We also recall that the vectors v (1) and v (2) which span the zero eigenspace of the linearization of the differential equation are written in this imaginary/real decomposition as v (1) = (0, 0, 0, p * 1 , p * 2 , p * 3 ) and v (2) = ((∂ ϕ p * ) 1 , (∂ ϕ p * ) 2 , (∂ ϕ p * ) 3 , 0, 0, 0).Thus, if we multiply Eq.( 19) by "i", we see that in the case of N = 3, Eq.( 19) is equivalent to (t φ + θ)v (1) + φv (2) = (0, 0, −εγp * 3 , 0, 0, 0) .Applying the projection operators n (1) | and n (2) | to this equation we obtain: , where the next to last equality used the formula for n (2) derived in Eq.( 16) and the last equality used the form of p * 3 derived in Theorem 2.2.Here, the approximate equalities mean that we have retained only the lowest order terms in ε in these expressions.
Note that this gives us an easy way to check the progress along the family of periodic states.Recall that our modulation hypothesis implies that One particularly easy thing to check numerically is the point at which the first component of the computed solution is purely real (i.e., when q 1 (t) = 0).If we denote the k th time at which this occurs as T k , then from Eq.( 20) we see that From this, we see immediately that Numerical computations of T k exhibited in Fig. 4 show clearly this scaling with √ k indicating that our modulation hypothesis is capturing the drift along the family of periodic orbits.However, Eq.( 21) indicates that this rotation of the phase should begin immediately, while our numerics indicates that there is some initial transient, and the √ k scaling appears only after some time.We are not yet sure what the origin of this transient behavior is, and plan to investigate it further in future work.
As a second check of the validity of this scenario we compute the decay of energy in the system as captured by the 2 norm of the solution.Define (p j (t) 2 + q j (t) 2 ) .
The undamped γ = 0 dynamics preserve this energy.If we differentiate with respect to time, the only contribution comes from the damping applied to the N th mode and we find Ė(t) = −γε(p N (t) 2 + q N (t) 2 ) .
Once again, we check this numerically in Fig. 5.

4.
Conclusions.In a simple model, we have shown that one can compute a family of breather solutions perturbatively.We have then shown the link of these breathers to the extremely slow rate of energy dissipation observed in such lattices of coupled nonlinear oscillators, by using a modulation hypothesis to compute the rate at which trajectories drift along the family of breathers in the limit of small dissipation, and we have then linked the predictions of the modulation theory to the rate of energy dissipation.The predictions of the theory agree very well with numerical experiments performed on this system./ log(ε).If m t decays like exp(−ctγε s ), then the calculation will lead to k = s.Indeed, we see that the decay rate is ε 2N −1 .The calculations shown were done for γ = 0.2.

Figure 2 .
Figure 2. The cylinder illustrates the set of periodic solutions of Eq.(9), with ϕ changing (very little) from left to right and the circle illustrating the angle ϑ.The spiral illustrates the way a time-dependent solution of Eq.(3) slides along the cylinder.It actually does not converge to it but will stay at some small, finite, distance from it.So the cylinder is Lyapunov stable in the sense of[10], at least as long as ε stays small.

Figure 3 .
Figure 3.The figure shows the γ dependence of the real parts of the eigenvalues, for N = 3 and ε = 0.01.The three curves are linear with intercept 0 and slopes 3.2 • 10 −10 , 0.0027, and 0.00727.Note that the first eigenvalue has an extremely small positive real part, while the others are stable.

Figure 4 .
Figure 4.This graph illustrates the behavior of p 1 (t), for N = 3 and several values of ε. and γ = 0.2ε.One measures the downcrossing times T k of the k th downcrossing of p 1 through 0. The theory predicts that X = T k /T k−1 − 1 /( k/(k − 1) − 1) = 1.As noted in the text, the transient behavior is not yet understood.

Figure 5 .
Figure 5.This graph illustrates the decay properties of the 2 norm as a function of time, for various values of N and ε.The horizontal axis is ε and the vertical axis is an estimate of the decay rate, obtained as follows: If m t is the 2 norm at time t and m t that at time t , then we compute k = k(ε) = log log(mt/m t ) γ•(t −t)