Metastable systems are common throughout multiple disciplines including biology, physics, machine learning and chemistry. In this paper we reexamine the existence and properties of metastable states in (finite) lattices of coupled nonlinear oscillators. We are specifically interested in a finite discrete nonlinear Schrödinger equation (DNLSE) chain in which energy is forced to flow through the system from one end to the other. Stochastically perturbed systems of this type have been used as models to understand the microscopic origins of heat transport and steady state behavior in non-equilibrium statistical systems [1]. A key question in such studies is how long it takes for the system to converge to its invariant measure. This controls how long it takes for the system to exhibit ergodic properties -in some cases making that time scale so long that these properties are e↵ectively unobservable [2] [3] [4] [5] [6]. While statistical mechanics methods have been applied to such systems, it has been been observed that the presence of breathers in the system prevents the definition of a well-defined thermodynamic temperature for the state. [7] [8]. In both numerical studies and theoretical investigations, it has been observed that the solution can be "trapped" for long times in small regions of the phase space. This long-time trapping of the solution in restricted regions of phase space is the metastability we will be studying, and in this paper we look for geometric structures in the phase space, and an analytical mechanism that can help explain this trapping.
In recent work, it has been proven that for chains of oscillators subject to very weak, localized damping, there are open sets of initial data in which this metastable behavior can be approximated for very long times by "breathers", i.e. temporally periodic, spatially localized solutions of the undamped equation [9] [10] [11]. The proof used a modulation approach in which the solution of the damped system was decomposed into a breather whose frequency and phase were allowed to vary (slowly) in time, and a correction term which was proven to remain small for very long times. A similar e↵ect is present in the experimental technique of self-localization, where damping is turned on at both ends of a dNLSE chain for a short period of time to prepare the system into a localized state after which the dissipation is turned o↵ [12] [13]. Thus, we expect that while employing this technique, the system is in fact in a metastable state similar to the one studied above and it is this metastability which causes the phenomenon of self-localization.
In the present paper we prove the existence of a new family of breather type solutions which give even better approximations to the metastable states.
We are specifically interested in modeling the types of states that have been found to inhibit the convergence of models of heat conduction toward their invariant measure [4] [5] [14]. These states are ones in which the motion of the system is approximately periodic, and there is a strong localization e↵ect of the system. In additional the dissipation in the system typically only occurs near the ends, either through an explicit dissipative term, or through the coupling of the end of the system to a heat bath. This localization of the dissipation is important, and as pointed out in [5](p.1001) is quite di↵erent from models in which the dissipation acts equally on all sites of the system. This new family of breather type solution is constructed, not by regarding the damped system as a perturbation of the undamped case, but instead by adding an additional very small perturbation to the damped system. We prove that this perturbed system has breather solutions, and we then analyze the stability of these solutions. We show that these new localized states (which we dub damped and driven breathers) reproduce previous predictions for the metastable behavior, but also provide more accurate approximations to the metastable states themselves. An additional advantage of these solutions is that they allow us to use standard techniques for analyzing the stability of fixed points in studying the metastable states.
There has been an immense amount of work on breathers in general nonlinear lattices, including the addition of damping and driving terms. Such systems arise naturally in optical cavities, for example, [15], [16]. The existence and properties of the breathers in such systems have been studied both numerically and analytically in these papers and others such as [17], [18], [19]. However, in these models the dissipation is spread throughout the system and this creates very di↵erent issues and e↵ects, both physically and mathematically. Indeed, one of the first rigorous investigations of breathers included the possibility of dissipation [20], but again, the types of systems considered in that work do not include those with localized dissipation of the type studied here.
Much closer to our work is that of Khomeriki [21], and Maniadis, et al [22]. In the first of these papers, the author studies a discrete NLS equation with driving and one end and damping at the other, very similar to the model we study. The author finds that for su ciently large driving, one can create moving solitons in band gaps for this model. For small driving amplitudes one finds stationary breathers similar to those analyzed here, but the focus of [21] is not on the sort of metastability issues considered in this paper nor is the stability of those solutions analyzed. Maniadis, et al, also consider a chain of nonlinear Hamiltonian oscillators driven at one end, though in this case they focus on a discrete Klein-Gordon chain. They also find numerically that one can generate a large family of moving breather type solutions in this model by varying the driving force, and they track the types of solutions generated by a numerically implemented bifurcation analysis. Again, the focus of this paper is not of the type of metastable behavior that we are interested in. Finally, one very interesting possible connection of our work to current research that was pointed out by a referee is to PT-symmetry breaking in nonlinear lattices [23]. If one chooses the strength of the damping and driving to be equal, then the linearization of our problem satisfies the PT symmetry property. Quantum mechanical systems with this type of perturbation present at two sites as a sort of "impurity" have been studied in [24] and [25]. There, the authors find that if the strength of the perturbation is su ciently large, one has PT-symmetry breaking and the appearance of exponentially growing or decaying states. One could then ask whether there is a bifurcation that occurs when the driving and damping di↵er from one another and whether or not this could be used to give an alternate proof of the existence of the sort of states we study. Since in our particular situation, the values of the driving and damping di↵er by several orders of magnitude, we haven't pursued this approach here, but it seems like an interesting question for future investigation. Impurity states have also been studied in systems without PT-symmetry (see for example, [26]) where they have again been shown to give rise to localized breathing modes. However, the mechanism in this case is rather di↵erent than the one we observe since it arises from exponentially localized eigenstates of the linearized problem created by the potential well associated with the impurity, whereas the states we study are perturbations of inherently nonlinear states.
We conclude the introduction with an overview of the remainder of the paper. Our goal is to understand the trapping of energy in localized regions of the lattice which is observed in studies of heat transport in coupled lattice of nonlinear Hamiltonian oscillators. These "trapped" regions correspond to states in which the energy is highly localized in particular spatial regions, and exhibits a very slow dissipation rate. With this goal, we first prove the existence of a class of spatially localized, temporally periodic solutions in a system subject to weak damping and driving. We then derive precise estimates on these states which we show reproduce aspects of the metastable or trapped states observed in the numerical studies of such systems. The states we construct occur in continuous families and we then analyze the stability of this family. In particular, we prove that the linearization of the system at each of these fixed points has one zero eigenvalue, corresponding to a symmetry of the system, one very small positive eigenvalue, and all the remaining eigenvalues have strictly negative real part. We then argue that these stability properties explain why solutions are trapped near these states for long periods -in particular, the eigenvalues with negative real part imply that nearby solutions are attracted toward the family of fixed points, while the small positive eigenvalue leads to a very slow drift along that family. Finally, we present numerics which show that solutions which start close to some member of one of these families of states do remain close to that family for very long times.
Following [9] [11], we work with a very specific system of coupled oscillators, namely the discrete nonlinear Schrödinger (dNLS) equation:
Note that if ✏ = 0, this system becomes N , uncoupled, nonlinear oscillators, and we then have trivial, localized, periodic solutions in which one site rotates with non-zero angular frequency, and all other sites have zero amplitude. This is sometimes referred to as the anti-integrable limit.
We add a driving term to the first site of our system, and a damping term to the last site by modifying the equation as
The addition of damping and driving to the equation allows one to study energy flow from one end of the system to the other. The existence of breathers in the dNLS equation with localized forcing but distributed damping was studied in [27]. Throughout the paper we will be referring to the = 0, = 0 case as undamped, 6 = 0, = 0 as damped, and 6 = 0, 6 = 0 as damped and driven. We emphasize again that the choice of our particular type of damping and driving is motivated by the desire to understand the very slow evolution of models of heat conduction toward their equilibrium state and on the basis of numerical experiments the most slowly decaying states appear to be those in which the energy is localized as far as possible from the sites at which the damping occurs. However, we believe that if one chose to impose the damping and driving at other sites, one would also observe localized periodic solutions of the type we construct -however they would correspond to more rapidly decaying states in the heat conduction models.
If we now translate to a rotating frame of reference z j (t) = e i!t ⇣ j , then
Finally, we will often want to work with real variables so we decompose ⇣ j = p j + iq j . Then (3) can be written as the system of di↵erential equations:
Past research on systems of this type has focused mainly on the Hamiltonian case in which the damping and driving are both absent -i.e. = = 0. In this case, there is an extensive theory investigating the existence and properties of "breathers". These spatially localized but temporally periodic solutions are known to exist for very general types of lattices of nonlinear Hamiltonian oscillators including the NLS equation, [28] [29] [30].
For the particular case of (4)-( 5), one has the following result:
Theorem 1 (see [11]; Theorem 1.6) If = = 0, there exists ✏ 0 > 0 and ⌦ > 0 such that for all |✏| < ✏ 0 , and |! 1| < ⌦, the equations (4)-( 5) have a family of fixed points (p ⇤ j (✏, !), 0) which vary smoothly with ✏ and !, and such that p
), j = 2, . . . , n, and q ⇤ j = 0.
Remark 1 Note that when we say that (4)-( 5) has a family of fixed points of the form (p ⇤ j (✏, !), 0), we mean that the q j components of the fixed point are zero.
Remark 2 Due to the invariance of (3) under complex rotations ⇣ j ! e i✓ ⇣ j , each of the fixed points constructed in Theorem 1 corresponds to a circle of fixed points.
Remark 3 Because the fixed points in Theorem 1 are constructed with the aid of the Implicit Function Theorem, for each (✏, !), they are the unique fixed points in a neighborhood of p 1 = 1, p 2 = p 3 = • • • = p N = 0 (up to the complex rotation noted in the previous remark).
Remark 4 Note that each of the fixed points constructed in Theorem 1 corresponds to a periodic solution with angular frequency ! for (2). If we fix ✏ and consider the family of solutions as a function of !, combined with the circle of solutions given by the complex rotations we see that in the phase space of (2) we have a cylinder filled with periodic orbits.
2 Existence and properties of damped and driven breather for NLS
2.1 The implicit function theorem and the existence of breathers for the damped and driven system
Coupled systems of nonlinear oscillators like (2) have often been studied as models for heat transport in solid state matter. In these studies, the first and last oscillator are coupled to stochastic heat baths with di↵erent temperatures, to force an energy flow through the system, and a damping term is added to insure that the total energy in the system does not grow without bound. It has been observed that the convergence to the invariant measure in such systems is often extremely slow due to the presence of metastable states in the phase space. In prior work, the energy flow through the system is enforced by placing a localized damping term at one end of the system and placing all the initial energy of the system at the opposite end [31] [9]. One can prove that the energy in the system tends to zero as t ! 1 and hence it must flow through the system in order to dissipate in this way. In [11] it is proven that for very long times during this dissipative process the system remains in a neighborhood of the cylinder of breather solutions for the undamped system.
In this paper we argue that a new family of solutions for a system subject to both damping and (very weak) driving provides even better approximations to the metastable states than do the undamped breathers. The existence of this family of breathers follows from the following theorem.
Theorem 2 For ✏ and su ciently small and all ! in a neighborhood of ! = 1, there exists a unique value of such that the equations (4)-( 5) have a stationary solution in a neighborhood of the undamped breather. This family of fixed points depends smoothly on , ! and ✏.
Remark 5 Note that stationary solutions of (4)-( 5) correspond to periodic solutions with angular frequency ! of (2).
The proof of this theorem follows from an application of the Implicit Function Thoerem. Details are provided in Appendix A.
We wish to show that the damped and driven breathers constructed in the previous subsection provide a better approximation to the metastable states that occur during energy dissipation in these systems than do the undamped breathers. To do so, we need accurate approximations to these solutions which we derive in this subsection. Because we proved the existence of these breathers with the aid of implicit function theorem, we know that they are analytic functions of the parameters ✏, and !, and hence we can derive power series approximations for them. We describe those approximations in this section and provide the details of their derivation in Appendix A.
In Appendix A, we define functions f j (p 2 , p 3 , . . . p N , q 2 , q 3 , . . . , q N ; ✏, !, , p 1 , q 1 ) whose solutions give the breather solutions. For example, f 2 (p 2 , p 3 , . . . , p N , q 2 , q 3 , . . . , q N ; ✏, !, , p 1 , q 1 ) = ✏(p 3 + p 1 p 2 ) !p 2 + (p 2 2 + q 2 2 )p 2 -i.e. the component of the vector field corresponding to q2 .
In order to derive approximations to the breathers, we consider f j (p 2 , p 3 , . . . , p N , q 2 , q 3 , . . . , q N ; ✏, !, , p 1 , q 1 ) = 0 (6)
for j = 2, . . . , N. We want to approximate the solutions of these equations when ✏ and are small and p 1 and ! are near 1.
Note the equations for f 2 , . . . f N 1 , f N +2 , . . . , f 2N 1 do not involve or , and hence are the same as those for the undamped breathers. Thus we can construct power-series approximations for these equations exactly as in the undamped case -i.e., p
. , while all of the q ⇤ j are zero to this order in ✏. (See Section 2 of [9] for more details.) Note that because the implicit function theorem guarantees that the solution is analytic and unique, if we can find a consistent power series approximation to the solution to a given order in ✏, , . . . , then it must correspond to the actual solution. In order to find the leading order terms in p ⇤ N and q ⇤ N , we must look simultaneously at the equations for f N and f 2N . We find that (to lowest order in ✏ and ) they satisfy:
Inserting our leading order expansion for p ⇤ N 1 , and inverting the matrix
we find:
Using this result for q ⇤ N , we construct the leading order terms for q ⇤ j for j = N 1, N 2, . . . , 2. For instance, consider the equation
Inserting the leading order expression for q ⇤ N from (9), we find
or
Continuing in this way, we find the leading order expression
for j = 2, . . . , N 1.
If we now substitute the leading order expression for q ⇤ 2 into the equation f N +1 = 0, (see (74) in Appendix A for the form of f N +1 ) we calculate that the unique value of up to leading order is:
Remark 6 The last equality used the fact that 2N 2 is always even and it shows that the value of is positive for any choice of N , and also extremely small when N becomes large. The positivity of is important and natural, since in order to have a steady state in the damped system we expect to need to inject energy at some point, which corresponds to a positive value of .
For the undamped breathers, the fixed points all have q ⇤ j = 0 which means that the z j all lie on the real line within the complex plane -although by rotational invariance, they can be rotated to any line through the origin. In contrast, the damped and driven breathers calculated in the previous section have q ⇤ j 6 = 0 for j > 1, when q 1 = 0. Therefore, the components of the damped and driven breathers do not lie on the same line through the origin. This is the twist of the breathers which we will explore in this section.
An appropriate form of analyzing the twist is by studying the norm and the complex phase of the equations z j (t) = r j (t)e i' j (t) . The utility of such a representation in the study of periodic solution of the undamped dNLS equation was noted at least as long ago as [32]. Since we are interested in studying the energy of the system, we can look at the energy of each site individually, defined by E j (t) = 1 2 |z j (t)| 2 . The system can be represented using these energies and the complex phase. Furthermore, due to rotational invariance we can use the di↵erences in the phases rather than the absolute phases to represent the system. The breather solutions correspond to stationary points of this representation. Calling j (t) = ' j+1 (t) ' j (t) for j = 1, 2, . . . , N 1, the system can be written as follows:
Remark 7 Note that an additional advantage of the energy-phase representation is that it uses the rotational invariance of the system to reduce the dimension of the system of equations from 2N to 2N 1.
Remark 8 In (15), we set the first term equal to zero if j = 1, and the second term equal to zero if j = N .
This energy-phase representation of the breathers is especially useful when we wish to compare the metastable states which emerge when the damped system is allowed to decay with our breather solutions. In particular, in the undamped, n-site case all the breathers have sin j = 0 for all j. This means, for example, that if we use the rotational invariance to set q ⇤ 1 = 0, then all the q ⇤ j will be equal to zero. However, for nonzero , the equations will lead to sin j 6 = 0.
For the particular breather calculated in Section 2.2, the phase is calculated from ' j = arctan(
. . . The phases for this breather will be:
while ' 1 = 0 as a consequence of choosing q 1 = 0. Using this, we find that j = ' j+1 ' j to leading order is simply:
In particular, we obtain the interesting following result for the last phase di↵erence:
Theorem 3 For a fixed > 0, ! > 0, the phase di↵erence between the last two oscillators, N 1 , is up to leading order independent of the number of sites. The expression for this phase is:
Remark 9 As we will see in the numerical experiments in Section 5, the metastable states we observe during the energy decay in these systems, also exhibit a non-zero "twist" of this type, whose value is very close to the prediction of (21).
In this section, we examine the stability of the damped and driven breathers constructed above, using the fact that they are fixed points of the equations ( 4)- (5). Define F (p 1 , . . . , p N , q 1 , . . . , q N , ; ✏, !, )
where the component functions f j are defined in (65)-( 67)-(72) in Appendix A. We recall that our goal in this stability analysis is to show that the stability properties of the family of breathers constructed in the previous section is consistent with the trapping of trajectories observed in studies of heat conduction.
If we linearize F at the the fixed points (p ⇤ , q ⇤ ) = (p ⇤ 1 , . . . , p ⇤ N , q ⇤ 1 , . . . , q ⇤ N ), the Jacobian matrix takes the form
Here, A (!,✏) and B (!,✏) are N ⇥ N tri-diagonal matrices and D (1) and D (2) are N ⇥ N diagonal matrices. More precisely, the matrix A (!,✏) has ✏ on the sub-and super-diagonal, while the diagonal elements are ! a j ✏ (p ⇤ j ) 2 3(q ⇤ j ) 3 , j = 1, 2, . . . , N, where a 1 = 1, a N = 1 and all other a j = 2. The matrix B (!,✏) has sub-and super-diagonal elements equal to ✏, while the diagonal elements are ! + a j ✏ + 3(p ⇤ j ) 2 + (q ⇤ 1 ) 2 . The diagonal matrix D (1) has diagonal elements D
(1) jj = 1,j N,j 2p ⇤ j q ⇤ j , while D (2) has matrix elements D
(2) jj = 1,j N,j + 2p ⇤ j q ⇤ j . Note that from the approximations to p ⇤ j and q ⇤ j derived in the preprint, we see that
). Thus, the diagonal part of the Jacobian is dominated by the dissipation coming from .
We begin by recalling a few key facts that we already know about the spectrum of the linearization which we will need in order to approximate the "small" eigenvalue. The proofs of these facts can be found in [11] 1. Due to the invariance of the original system of equations under complex rotations, there is a zero eigenvalue with eigenvector
2. The linearization also has 2N 2 simple eigenvalues near ±i which all have negative real parts that are O( ). Denote these eigenvalues by 3 , 4 , . . . , 2N and the associated eigenvectors by v (j) , j = 3, 4, . . . , 2N . Note that in [11] this is proven for the linearization about the breather with = 0. However, since all of these eigenvalues are simple, this remains true for su ciently small by simple, non-degenerate perturbation theory estimates.
We also identify another vector
which corresponds to the tangent vector in the direction of changes in frequency of the breathers. While not an eigenvector, there are two important properties of ṽ(2) that we use in what follows:
1. The set of vectors {v (1) , ṽ(2) , v (3) , v (4) , . . . , v (2N ) } form a basis for R 2N . This follows from the fact that they are small perturbations of the eigenvectors and generalized eigenvector of the linearization of the undamped breathers studied in [11].
2. We know that when the linearization at the breather acts on ṽ(2) we have
This follows by di↵erentiating the equation for the damped and driven breather with respect to !.
From our approximations to the breather we have
From the form of ṽ(2) , and the expression for @ ! p ⇤ 1 we have 0
Remark 10 To simplify the following expressions, we choose ! = 1 for the remainder of this section, but the case of general ! is very similar.
Then, if we expand (p ⇤ 1 , 0, . . . , 0) T in terms of the basis {v (1) , ṽ(2) , v (3) , v (4) , . . . , v (2N ) }, we have
with ↵ 2 = 2 + . . . , and all the other ↵ j at least O(✏) or smaller. Note that the coe cients ↵ j are, in principle, computable, at least approximately, since we have explicit formulas for v (1) and ṽ(2) , and since the remaining v (j) 's correspond to simple eigenvalues, they can be computed approximately by ordinary non-degenerate perturbation theory, at least to leading order.
We now denote the small, non-zero eigenvalue as 2 , and write its eigenvector as
Note that we compute below that the µ j 's are all small and go to zero as ! 0. This is expected since for = 0, we know that the linearization at the breather has a twodimensional zero eigenspace, but only one zero eigenvector, v (1) [11].
Now consider the eigenvalue equation for 2 :
Insert the expansion from (30) into both sides of this equation and use the fact that on the left hand side we know that D (p,q) F (p ⇤ ,q ⇤ ) v (1) = 0, D (p,q) F (p ⇤ ,q ⇤ ) v (j) = j v (j) , for j = 3, . . . , 2N , and D (p,q) F (p ⇤ ,q ⇤ ) ṽ( 2) is given by (26). Thus, we have
We now equate the coe cients of the various basis vectors on the two sides of this equation.
From the coe cients of the vector ṽ(2) , we have
Inserting our approximate expressions for ↵ 2 and (@ ! ), we find the approximate expression for the small, positive eigenvalue:
(Recall that we have set ! = 1.) While they are not really important, we can also solve for the remaining coe cients µ j . Thus, equating the coe cients of v (1) in (33), we find
while
Remark 11 Note that for N = 2, the formula in (34) agrees with the computation for the two-site breather in energy-angle coordinates in (61)
From this expression, we obtain µ j << µ 2 for j 3. Thus, from (30), the small positive eigenvector is (up to very high order) in the span of the vectors v (1) , ṽ(2) which corresponds to the tangent space of the cylinder of breathers, with v (1) corresponding to rotations around the cylinder and ṽ(2) to the direction along the family of breathers. Hence, the positive eigenvalue pushes solutions to slide along the family of breathers, causing ! to decrease over time. Movement perpendicular to the cylinder is negligible until ✏ ⇠ !, at which point v (2) has a significant component in the normal direction (µ j ⇠ µ 2 ) forcing the solution to finally exit the neighborhood of the cylinder.
We next consider how the spectral picture derived in the previous subsection can be used to model the metastable evolution in the damped (but not driven) system of oscillators. We will see that we reproduce the picture obtained from the modulation approach of [11], by considering the evolution near one of the fixed points represented by the damped and driven breathers. However, as we discuss, these damped and driven breathers provide a more accurate approximation to the metastable state than the undamped breathers used as an approximation in the prior work.
To simplify the exposition, we introduce the following notation. We set u = (p, q) T , where p and q are N -vectors. Likewise, the damped and driven breather will be denoted by u ⇤ = (p ⇤ , q ⇤ ) -depending on the context, we may also include the dependence of the breathers on the parameters -i.e. u ⇤ = u ⇤ ( , !). We also introduce two diagonal matrices which represent respectively the driving and damping -B( ) is a 2N ⇥ 2N matrix with (1, 1) and ((N + 1), (N + 1)) element equal to and all other elements zero, and ( ) is 2N ⇥ 2N matrix with (N, N ) and (2N, 2N ) element equal to and all other elements zero. Finally, we let H(u) denote the Hamiltonian of the undamped and undriven lattice of oscillators. We can then write the damped, but not driven, evolution as
The damped and driven breathers u ⇤ are solution of the system of equations
Here, J is the standard, symplectic matrix for Hamiltonian systems, which takes the form
, with 0 representing an N ⇥ N matrix of zeros, and 1 an N ⇥ N identity matrix.
Recall that u ⇤ is the ✓ = 0 member of a family of fixed points of (38) of the form
Because of the equivariance of the linearized equation under these complex rotations, the spectrum of the linearization at the breather solution is independent of ✓, and hence, without loss of generality, we can assume that ✓ = 0.
We will consider a solution of (37) with initial conditions near the breather. We will expand the solution with respect to the basis of vectors {v (1) , ṽ(2) , v (3) , . . . , v (2N ) }, where as in Section 3, v (1) = ( q ⇤ , p ⇤ ) is the eigenvector of the linearization about the damped and driven breather with eigenvalue zero, v (j) , with j = 3, 4, . . . , 2N are the eigenvectors with eigenvalues near ±i (but all with real parts ⇠ O( )) and ṽ(2) = @ ! u ⇤ . We know from our discussion of the spectrum of the linearization about the damped and driven breather that the eigenvector with small positive eigenvalue is approximately v (1)
We take an initial condition for (37) of the form
We assume that w 0 is:
1. small, and 2. has no component in the v (1) or ṽ(2) direction.
The second of these conditions can be enforced by adjusting the angle ✓ in (39) and the frequency !. Because of the equivariance of the linearization about the damped and driven breather with respect to the complex rotations by ✓, the spectrum of the linearization is independent of ✓ (and depends only weakly on !) so we will assume that we are linearizing about the breather with ✓ = 0 and ! = 1.
We look for a decomposition of the solution similar to that in (40) -that is, we write the solution of (37) as
Inserting this into (37) we find
The goal in rewriting the equation for ẇ in this form is to note that from the definition of u ⇤ , the expression in (42) vanishes, while if we apply Taylor's Theorem to expand (43) in w, the linear terms in w just give L u ⇤ w, the linearization about the damped and driven breather. Thus, we have
where Ñ (w) collects the nonlinear terms resulting from applying Taylor's Theorem to (43), and is O(w 2 ).
We can write the solution as
We will ignore the contributions of ↵ j (t), with j 3, since we expect the components of the solution in those directions to be contracting since the eigenvalues corresponding to these eigendirections all have negative real part.
To isolate the evolution of ↵ 1 (t) and ↵2 (t), we also introduce the vectors
for constants c 1 and c 2 determined below. Note that if we consider the transpose of L u ⇤ , we have
and
The first vector on the right hand side of (48) is the 2N vector whose first N components are all zero, and whose second N components have (@ ! )p ⇤ 1 in the (N + 1) st component and then zero in the remaining components. We choose the constants c 1 and c 2 so that ñ(1)
• v (1) = n (2) • ṽ(2) = 1. Note that from the definitions of these vectors, and the asymptotic formulas we derived for p ⇤ and q ⇤ , we find that for ! = 1, c 1 = 2 + O(✏),
• ṽ(2) and n (2) • v (1) are both O( ✏ 2N 1 ).
Thus, if we take the dot product of the equation ( 45) with ñ(1) , we have (ignoring the nonlinear terms, for the moment) ↵1 = ↵ 1 hñ (1) , L u ⇤ v (1) i + ↵2 hñ (1) , L u ⇤ ṽ( 2) i (50)
hñ (1) , Bu ⇤ i hñ (1) , B(↵ 1 v (1) + ↵2 ṽ( 2) )i ↵2 hñ (1) , ṽ(2) i
, ṽ(2) i hñ (1) ,
Here we used the form of ñ(1) and u ⇤ to bound |hñ (1) , Bu ⇤ i| by O(✏ 4N 2 ).
If we now repeat this procedure, taking the dot product of (45) with n (2) , we find ↵2 = ↵ 1 hn (2) , L u ⇤ v (1) i + ↵2 hn (2) , L u ⇤ ṽ(2) i ↵1 hn (2) , v (1) i (52) hn (2) , Bu ⇤ i hn (2) , B(↵
If we ignore the high order correction terms, we can solve these equations and we find
and
By adjusting the breather u ⇤ about which we are perturbing, we can insure that ↵ 1 (0) = ↵2 (0) = 0, so we have
, and ↵2 (t) = 2 t .
(56)
Inserting the value of , we have
Recall that our parameter corresponds to ✏ and our ↵2 corresponds to in the paper [9]. (This last equivalence follows from the fact that in [9] is the change in the frequency of the breather, while our ↵2 corresponds to the displacement in the direction @ ! u ⇤ -i.e. also the change in frequency of the breather.) Then for N = 3, the case considered in that reference, the formula for (t) derived in that work agrees exactly with our formula for ↵2 .
To compare the expression for ↵ 1 (t) with that of ✓(t) in [9], note that ↵ 1 (t) is the coe cient of v (1) , the eigenvector corresponding to rotations in the complex plane u ! e i✓ u. (This is not the same ✓ as in [9]). In [9], the quantity corresponding to a complex rotation of this sort is t (t) + ✓(t) (see equation (13) of that reference). If we take the formulas ✓(t) and (t) from [9] and compute, we obtain
which again, agrees exactly with our formula for ↵ 1 (t) when N = 3. (Again, recalling that our should be replaced by ✏ to compare with the results of [9].) This is only an approximate calculation -in particular, we have ignored the contributions of the nonlinear terms, but in [11] it was shown rigorously that in the case of damping alone, the nonlinear e↵ects could be rigorously incorporated without changing the predictions of the linear approximation and we expect that the same is true here.
Table 1: Numerically computed spectrum of the linearization about the breather, when ✏ = 0.03, = 0.0035, and ! = 1. 0.00257707 + 0.987135i 0.00257707 0.987135i 0.000922942 + 0.914551i 0.000922942 0.914551i 3.43099 ⇤ 10 8 + 0.i 2.28602 ⇤ 10 9 + 0.i
In this section, we present numerical experiments which illustrate the results of the preceding sections. We will compute the spectrum of the linearization about the breather, showing that it has the form described in the previous section, and we will compute the dynamics of the metastable state of the damped system. In addition to showing that solutions of the damped system, which start close to the damped and driven breathers remain close to this family for a long time we will also follow such solutions until their disappearance and develop a conjecture about how such families of metastable states terminate. We focus on systems with a small number of sites (specifically, N = 2, 3) because the very slow drift along the family of breathers means that the computational time required increases very rapidly with N . This makes it unfeasible to compute the solution up until the disappearance of the metastable state for larger N . In addition, by concentrating on these small systems we are able to perform many computations explicitly, which helps to elucidate the numerical results. All numerical computations are done with Mathematica.
We begin by computing the spectrum of the linearization of equations ( 4)-( 5) about the breather solution. Table 5 contains the spectrum of the linearization around the breather solution with ✏ = 0.03, = 0.0035, and ! = 1. Note that as expected, there are two eigenvalues near i, and two more near i, all with negative real part, one with small, positive real part and one eigenvalue which is zero to within the errors of our computations (i.e. ⇠ O(10 9 )),
As we saw in the previous section, the crucial eigenvalue for understanding the metastable behavior of the system is the small positive eigenvalue, so we focus on that. In Figure 1, we plot the numerically computed value of this eigenvalue for N = 3, = 0.0035, and for five di↵erent values of ✏, ✏ k = 0.03 + 0.002 ⇤ k, k = 1, 2, . . . , 5, and compare these values to the values predicted by equation (34). As we see, the computed values are quite close to the theoretically predicted values, and the approximation improves as ✏ becomes smaller. We next compute the full time evolution of the three site system to examine how long the metastable behavior persists, and how it terminates. For these computations we found it easier to work in the energy-phase representation of the system, (15), since this reduces the number of equations we need to solve by one. We choose as initial conditions for our numerics the approximation to the damped and driven breather computed in Section 2.2. We also need some way of determining what member of the breather family we are close to at any given time in the course of the the computation.
We proceed as follows. Choose a value of ! 0 to determine the initial values for the equations. Then set p 1 = p ! 0 . and approximate the initial values for the other variables using the equations of motion and our knowledge of the breather solution.
As the solution evolves, it moves along the family of breathers, and we need to determine how the frequency (which determines which breather we are close to) changes. If we are very close to a damped and driven breather we can consider ˙ j ⇡ 0 and so all the 'j are the same. For this reason, we can define !(t) to be any of the 'j and we should obtain approximately the same result. For simplicity, we define it as !(t) = '1 (t) and in turn we can calculate it from:
Notice that ! 0 6 = !(0) because of the fact that we calculate an approximation to p 1 from ! 0 and then use it to calculate approximations for E 1 (0), E 2 (0), 1 (0). In fact
The first quantity we examine are the graphs of 1 (t) and 2 (t) for the three site system. We recall that for the undamped breathers, all of the j (t) are zero, while the damped and driven breathers have non-zero "twist". Figure 2 graphs 1 (t) for the three-site system with two approximations. The blue curve represents the computed value of j (t). From Theorem 3, we know that the lowest order approximation to the twist of the breather is j (t) ⇡ !(t) , and this is plotted in orange, using the approximation !(t) ⇡ '1 (t). We obtain a good approximation to the computed twist, showing that the damped and driven breathers provide a better approximation to the metastable states than the undamped breathers. Finally, recalling that !(t) ⇡ 2E 1 (t) we have also plotted (in green) the quantity E 1 (t) . Interestingly, this quantity seems to provide an even better approximation to the numerically computed twist though at the moment we have no theoretical explanation of why this should be the case.
Speaking colloquially, our numerics show that the metastable state slides along the family of damped and driven breathers. This is consistent with our computations of the spectrum of the linearization which showed that the span of the eigenspace of the zero eigenvector and the eigenvector with small positive eigenvalue was very close to the tangent plane of the family of damped and driven breathers. Moreover, this allows for the modeling of the dynamical evolution of the metastable state of the damped system from the damped and driven breathers (especially the time evolution of the phases) for long times. We believe that the period for which this approximation is accurate is connected to the bifurcations The energy of each oscillator in an eight-site system. The parameters chosen were = 0.05, ✏ = 0.1, and the initial condition is p 1 = 1 and all other p j , q j = 0. of the undamped system.
In Figure 3, we present numerics for a system with N = 8 sites, which show that the e↵ects we have described persist for systems with larger numbers of oscillators. Again we note that once the system has fully settled in the metastable state, we observe an initial, almost flat state where the system moves along the family of damped and driven breathers for a very long time. The system eventually reaches a point where the energy of the first oscillator decays suddenly while that of the other oscillators rise to meet it until they have the same energy, after which they reach the thermalized state and decay exponentially. The sudden decay in the energy is significantly more dramatic for the eight sites than it was for two or three sites.
Thus, we conjecture that the general evolution of these metastable states is as follows.
When the trajectory is "captured" by one of the damped and driven breathers, the system begins to evolve by drifting along this family until the family of breathers disappears through bifurcation, at which point the system may be captured by another attractive family of approximately periodic orbits. There may, of course be other, possibly large, parts of the phase space where these metastable solutions do not exist and the evolution of the system is chaotic. In particular, it is not clear what properties of the periodic orbit in the undamped Hamiltonian system result in the solution becoming attractive once the system is subjected to damping.
We now turn to a discussion of the two-site damped and driven system where many of the phenomena discussed above can be illustrated by explicit computations. A detailed derivation is found in Appendix B, while in this section we will mostly show the results.
The family of breathers of the Damped and Driven system related to the metastable state we have studied is analytically calculated to be:
We can compute the spectrum of the linearization around these breathers, first writing out the Jacobian and then using Mathematica to derive series approximations. The small positive eigenvalue is:
The zero eigenvalue, 1 , is absent in the energy-phase representation, because it corresponds to uniform phase rotation at all sites, and the energy-phase representation considers only phase di↵erences between adjacent sites and hence is insensitive to phase rotations of the whole system. The leading order term in the small positive eigenvalue 2 has the same form as we derived in equation (34), if we set N = 2 in that formula, as expected.
In Appendix B, we derive an approximation to the time evolution of the metastable state using the Ansatz that the solution maintains the relation E 1 (t) E 2 (t) = p E(t) 2 ✏ 2 resulting in an expression for the phase shift:
Note that this is only an approximation to the solution -this value of will not satisfy the equation for ˙ . However Figure 4 shows that it is a very accurate approximation when
This result can then be used to calculate an approximation to the time evolution of the energy with initial energy E(0) = E 0 :
Here W 1 is the 1 branch of the product logarithm, and ⌧ is the time for which lim t!⌧ E(t) =
|✏| and is given by:
If (63) was exact, ⌧ would be the time when E = |✏|, corresponding to the bifurcation for the undamped system, therefore the solution would be valid up until ⌧ . Since this is not the case, ⌧ instead serves as an approximation for the time at which (63) stops working. In fact, the approximation stops working a small amount of time before ⌧ . This is expected since (62) does not make sense for E > q ✏ 2 2 4 due to the domain of the arcsin.
Figure 4 shows the e↵ectiveness of the approximation. E 0 was chosen so that it coincided with the initial energy of the graphs in the previous sections, leading to ⌧ = 4700.8. The energy was graphed for longer than ⌧ to identify the behaviour after it leaves the breather, while the phase di↵erence was graphed for less than ⌧ since the expression for the approximation explodes due to the singularity in the arcsin. After this, the behavior of the solution is governed either by (84) or (85). Note that by approximating the energy of the solution by the expression for the energy of the breather in (87) we are assuming that this approximation moves along the family of breathers and from Figure 4 we again see that this scenario accurately reproduces the observed numerical behavior of the solution of the system.
The simplicity of the two-site system allows for a more global investigation of the phase space. Once again, the analysis is better suited to using the energy-phase equations rather than expressing the system in the p, q coordinates, and the energy-phase equations for the two-site system give the simple system of three ODEs:
For the undriven-undamped, two-site system, (i.e. = = 0), all breather solutions were found in [32] and in our notation they take the form:
when ✏ > 0 with similar expressions for ✏ < 0. For each of these breathers, we can calculate the angular frequency ! from E. The values are ! (1) = E, ! (2) = E + 2✏, ! (3) = ! (4) = 2E + ✏. The third breather is the one close to the metastable states discussed so far.
These breathers continue to the damped and driven setting , and again, one can derive explicit formulas for them by calculating the fixed points of (76) -for example, if 0 < << , and < ✏ 2 , the family of breathers E (3) , continues to a family of the form described by (60) (Note that to compare these formulas with those derived for the general N -site damped and driven breather, one should first fix and ✏ and then choose so that E
(3) 1 = 1/2, since those breathers were chosen to be close to the breather with p 1 = 1 and all other components zero.)
If we linearize energy-phase system system (76) about this breather, the Jacobian matrix takes the form:
Using Mathematica to derive series approximations for the solutions of the characteristic polynomial, we find the eigenvalues of the linearization are the eigenvalue 2 , given by (61) and:
These two eigenvalues are close to ±i!, but with strictly negative real part, proportional to , as expected. An important fact to notice is that if we instead find the linearization about this fixed point using the p, q coordinates, we obtain the same eigenvalues 2 , 3 , 4 , with an added 0 eigenvalue corresponding to the invariance with respect to rotations about the cylinder of breathers.
For |✏| >> > 0 but = 0, the breathers disappear but they can be used to help find approximate solutions of the damped problem, and to understand its metastable behavior.
For example, if we make the Ansatz that E 1 (t) = E 2 (t) for all time, and substitute this into the equations of motion we find exact solutions: where E 0 = E 1 (0)+E 2 (0) is the initial total energy, and ✏ > 0. One has similar expressions when ✏ < 0. These are the states corresponding to (77) and (78) respectively, with only the first of these being stable.
We can use a similar argument to approximate the e↵ects of initial conditions close to the solution E (3) of (77). First, di↵erentiate d dt E
(3) 1 and d dt E
(3) 2 :
d dt
Using Ė = Ė1 + Ė2 = 2 E 2 , and considering E 1 = E
(3)
1 , E 2 = E (3)
Solving for sin we get (62). However, if we calculate ˙ from this expression, it will not match the expression in the equations of motion, meaning that the solution is not exact. This is expected due to the metastability of the state, the solution cannot be described by these equations for all time. Numerically, the above expression provides a very good approximation for long time. Recall that for this breather ! (3) = 2E + ✏, so calculating the leading order in we find = ! + . . . , which is identical to the value of the twist in the damped and driven breather from Theorem 3. This is further evidence of the fact that the damped and driven breathers provide a better approximation to the metastable states than do the undamped breathers.
We also obtain an approximation to the time evolution of the energy through:
If the system starts with initial total energy E 0 , an approximate expression we simply solve the integral:
If we set the upper limit to be ✏ and the right hand side to be ⌧ , we obtain (64), using the indefinite integral:
To obtain (63), it is worth defining:
Then, evaluating the indefinite integral (89) at E(t), E 0 , will lead to an algebraic equation in f (t):
Which can be shown to solve to:
From where (63) is derived, and ⌧ defined as in (64). To get that the branch of the product logarithm, we take the limit: lim