Stability of high-energy solitary waves in Fermi-Pasta-Ulam-Tsingou chains

Stability of high-energy solitary waves
in Fermi-Pasta-Ulam-Tsingou chains

Michael Herrmann111University of Münster, Germany,    Karsten Matthies222University of Bath, United Kingdom,
July 15, 2019

The dynamical stability of solitary lattice waves in non-integrable FPUT chains is a longstanding open problem and has been solved so far only in a certain asymptotic regime, namely by Friesecke and Pego for the KdV limit, in which the waves propagate with near sonic speed, have large wave length, and carry low energy. In this paper we derive a similar result in a complementary asymptotic regime related to fast and strongly localized waves with high energy. In particular, we show that the spectrum of the linearized FPUT operator contains asymptotically no unstable eigenvalues except for the neutral ones that stem from the shift symmetry and the spatial discreteness. This ensures that high-energy waves are linearly stable in some orbital sense, and the corresponding nonlinear stability is granted by the general, non-asymptotic part of the seminal Friesecke-Pego result.

Our analytical work splits into two principal parts. First we refine two-scale techniques that relate high-energy wave to a nonlinear asymptotic shape ODE and provide accurate approximation formulas. In this way we establish the existence, local uniqueness, smooth parameter dependence, and exponential localization of fast lattice waves for a wide class of interaction potentials with algebraic singularity. Afterwards we study the crucial eigenvalue problem in exponentially weighted spaces, so that there is not unstable essential spectrum. Our key argument is that all proper eigenfunctions can asymptotically be linked to the unique bounded and normalized solution of the linearized shape ODE, and this finally enables us to disprove the existence of unstable eigenfunctions in the symplectic complement of the neutral ones.

Keywords: Hamiltonian lattices waves, high-energy limit of FPU or FPUT chains, nonlinear orbital stability, asymptotic analysis,
MSC (2010): 37K60, 37K40, 70H14, 74H10

1 Introduction

This paper deals with the long-time dynamics in a prototypical example of Hamiltonian mass-spring systems, namely the nonlinear atomic chain introduced by Fermi, Pasta, Ulam, and Tsingou in [FPU55]. A key question in this context is whether coherent states persist or to which extent energy is transferred to higher modes leading to the process of thermalization. Various types of coherent motion exist for Fermi-Pasta-Ulam-Tsingou (FPUT) chains (exact breathers, different types of traveling waves) but it is in general not known whether solutions that start near such coherent states remain coherent for all times. Pioneering work for solitary waves was done by Friesecke and Pego in a series of papers [FP99, FP02, FP04a, FP04b]. They set up a general framework for nonlinear orbital stability and established the crucial linear stability for small-amplitude, long-wavelength traveling waves with speed just above the speed of sound, for which the Korteweg-de Vries (KdV) equation is the appropriate scaling limit. There are many extensions of their results as detailed below, but these still deal with situations near the KdV limit or the special situation of the Toda lattice, which is completely integrable. In this paper we study high-energy solitary waves, which are localised on an essentially single moving particle close to a hard-sphere collision, for a wide class of inter-action potentials far away from any obvious continuum approximation. We will provide orbital stability results for open sets of initial data near a universal family of high-energy traveling waves.

The FPUT chain describes an infinite chain of masses which are separated by nonlinear springs and interact according to Newton’s law of motion via


Here, denotes the position of particle at time and describes the potential for the nearest neighbor forces. Using atomic distances and velocities with

we obtain the first-order formulation


In this paper we study coherent behavior in the form of traveling waves, which satisfy the ansatz


where the profiles functions and have to be determined in dependence of the wave speed and as functions in , the space variable in the co-moving frame. Notice that the additional shift in the distance component has been chosen for convenience. It guarantees the existence of even wave profiles and allows us to write the traveling wave equation as


where all spatial differences are expressed in terms of the centered nabla operator


Moreover, eliminating the velocities in (1.4) we find the second-order advance-delay-differential equation


with discrete Laplacian


The existence of different types of waves (periodic, homoclinic, heteroclinic) has been proven for many potentials by several authors and using different methods. We refer to [Her10, HM15] for an overview but emphasize that rigorous results concerning the uniqueness and the dynamical stability of FPUT-type lattice waves exist only for the few integrable examples and in the universal KdV regime, see [Tes01] for an overview concerning the Toda chain, and [DHM06] for the harmonic case and the hard-sphere limit. The KdV limit has been introduced formally in [ZK65] and studied rigorously in [FP99, FP02, FP04a, FP04b]. The stability of a single solitary wave with respect to the natural energy norms is proven in [Miz09] and the linear stability part has been simplified in [HW13b]. The analysis has been been extended to interacting KdV waves in [HW08, HW09] and -solitary waves are studied in [Miz11, Miz13]. Further approaches to the stability problem are given [HW13a, KP17] and the orbital stability of large amplitude waves is shown in [MP08, BHW12] for the Toda lattice. Modulation equations and the related concept of long- but finite-time stability of KdV-type solutions are discussed in [SW00] for FPUT chains, and in [CBCPS12, GMWZ14] for other atomistic systems. Finally, more general existence proofs for KdV waves can be found in [FML15, HML16, HW17].

Notice also that the mathematical theory of traveling waves in spatially continuous Hamiltonian systems is much more developed. First, traveling waves in such systems are determined by ODEs instead of advance-delay-differential equations and hence often known almost explicitly. Secondly, the stability investigation benefits on the linear level from the rather extensive knowledge about linear differential operators (with non-constant coefficients) and on the nonlinear side from the shift symmetry, which gives rise to a further conserved quantity according to Noether’s theorem and provides easy control concerning accelerations of a given wave. Consequently, the orbital stability of solitary waves is now well-understood for a huge class of dispersive PDEs, see for instance [Ang09] for an overview.

1.1 Setting of the problem and main results

In this paper we suppose that the interaction potential possesses a repulsive algebraic singularity as illustrated in the left panel of Figure 1.3 and rely on the following standing assumption.

Assumption 1 (interaction potential).

The potential is a strictly convex -function on the semiopen interval . Moreover, it admits a global minimum at with

and posses at a normalized algebraic singularity of order . The latter means

and hence

for some real-valued exponents and and some constant . For technical reasons we also assume as well as .

Prototypical examples of Assumption 1 are


with and free parameter , as well as the Lennard-Jones-type potential


with free parameter corresponding to , where the classical choice is . Further examples and a more detailed discussion concerning and is given below in §1.2, and we mention that can be extended to a smooth and strictly convex function on .

We further restrict our considerations to waves that propagate with high speed and must hence carry a large amount of energy


In this asymptotic regime, the wave profiles and from (1.4) can be approximated by the piecewise constant functions


which are – up to elementary transformations – naturally related to the traveling waves in the hard-sphere model of elastically colliding particles, see Figures 1.1 and 1.2 for numerical examples and Figure 1.3 for a schematic representation of the limit dynamics. The high-energy limit of traveling waves in singular potentials has been studied in [FM02] within a variational setting and in [Tre04] by fixed-point arguments. Moreover, the authors of the present paper refined the underlying asymptotic analysis in [HM15] and showed that the fine structure of the wave profiles is determined by an asymptotic shape ODE, which provides more accurate approximation formulas for the wave profiles as well as the leading order terms in the scaling laws for all other quantities of interest. Related results for non-singular but rapidly growing potentials are discussed in [Her10, Her17].

However, [HM15] does not guarantee that the exact wave profiles depend smoothly on the wave speed but, roughly speaking, only that the leading order terms in the expansion around do so. Similarly, the asymptotic estimates control the approximation error but not the properties of . Since both issues are essential for our spectral analysis, we reinvestigate in §2 the high-energy limit on the level of the nonlinear profile functions and present, as explained below, a novel and more robust asymptotic approach to the shape ODE, which also provides an analogue to the local uniqueness results from [Tre04, HM17]. Our main findings are formulated in Theorem 7, Theorem 9, Lemma 8, and Lemma 15, and can informally be summarized as follows.

Main result 1 (family of high-energy waves parameterized by speed ).

Under Assumption 1 there exists a family of solitary solutions to the traveling wave equations (1.4) such that the following properties are satisfied for all sufficiently large :

  1. Qualitative behavior for fixed . The profiles and are even and exponentially decaying. Moreover, takes values in and is non-negative.

  2. Local uniqueness. The solution branch is isolated in natural function spaces.

  3. Smooth parameter dependence. The map yields a smooth curve in some Sobolev space with exponentially growing weights for .

  4. Accurate approximations. The profiles as well as and can be approximated by the unique solution to an asymptotic ODE initial value problem and the solution space of its linearization.

Figure 1.1: Distance profiles in the high energy limit for the Lennard-Jones potential (1.9) with corresponding to and in Assumption 1: Numerical solutions (black, dashed) to the nonlocal equations (1.12) along with the ODE approximations from Lemma 5 (gray, solid) and the tent-shaped limit from (1.11) (dark gray, thin line). All profile functions are plotted against , the unscaled space variable in the comoving frame. The precise approximation results are formulated in Theorem 7 and Lemma 8, where corresponds to , and the approximation error is plotted in Figure 2.3.
Figure 1.2: Normalized velocity profiles corresponding to Figure 1.1, which approach the indicator function of the interval , see (1.11). The numerical scheme to compute the nonlinear waves relies on the maximization of the potential energy under the constraint , is easy to implement, and provides traveling waves depending on the free parameter , see [Her10, HM15] for the details. Our existence result §2 is not based on constrained optimization but on explicit approximation formulas and nonlocal fixed point arguments, and yields waves parameterized by . In the high energy limit, both approaches are equivalent since we have and hence .
Figure 1.3: Left panel. Potential as in Assumption 1. Right panel. Cartoon of a solitary wave in the hard-sphere model, which propagates by ellastic collisions of a single active particle (gray) with a resting background (black). The corresponding distance and velocity profiles are affine transformations of the functions and from (1.11) and hence naturally related to the high-energy waves from Figures 1.1 and 1.2.

The asymptotic analysis of the stability of high-energy waves is presented in §3 and follows the general framework proposed by Friesecke and Pego. In particular, we study the stability of waves in exponentially weighted spaces and in an orbital sense. The latter means that certain initial perturbations might shift or accelerate the wave, and this requires long-time adjustments of the phase variable and the wave speed. All other perturbations, however, remain small, fall behind the wave, and decay exponentially in weighted norms centered around the moving bulk of the wave. We formulate our main result as follows and refer to §1.2 for more details concerning the proof strategy and the key arguments in the asymptotic analysis.

Main result 2 (nonlinear orbital stability in the sense of Friesecke and Pego ).

For any sufficiently large speed , the corresponding high-energy wave from Main Result 1 is stable in the following sense, where

denotes the wave profiles and where the positive constants , , depend on both and : Given FPUT initial data with

for some and , there exist a unique wave speed and a unique shift with

such that the solution to the FPUT chain (1.2) satisfies the -stability estimate

as well as the weighted decay estimate

for all . Here, stands for the integer variable .

1.2 Discussion and asymptotic proof strategy

The key observation for our asymptotic analysis of solitary waves with high energy is that the advance-delay-differential equations (1.4) can be replaced by an asymptotic ODE problem which dictates the fine structure of the distance profile near (’tip of the tent’) and provides in turn all relevant properties of lattice waves in the limit . This idea will be explained in §2.1 on a heuristic level and has already been exploited [HM15]. However, in §2 we present a modified approximation strategy which is likewise based on the asymptotic shape ODE but has the following advantages.

  1. The high-energy waves in [HM15] were constructed as solutions to a constrained optimization problem depending on a small parameter , and the wave speed could not be described a priori but came out implicitly as Lagrangian multiplier. In particular, in the variational setting we were not able to prove that the relevant wave quantities depend smoothly on the wave speed. In the modified approach, however, is a free parameter and we do not rely on optimization techniques. More precisely, we first identify in §2.2, §2.3, and §2.4 for any sufficiently large an approximate solution and show afterwards in §2.5 and §2.6 the existence of a unique exact high-energy wave in a small neighborhood within a certain function space. In this way we do not only establish the smooth parameter dependence, see §2.7 and §2.8, but can also provide accurate asymptotic formulas for the derivatives of the wave profiles with respect to . The latter are naturally related to the neutral Jordan modes in the linearized FPUT equation and feature prominently in the spectral analysis within §3.

    We also mention that the numerical illustrations in Figure 1.1 and Figure 1.2 are computed by a straight forward implementation of the constrained optimization problem. We refer to [Her10] and the caption of Figure 1.2 for the details and recall that the uniqueness results in §2.6 and [HM17] guarantee that both the variational and the non-variational approach provide for large the same family of high energy waves.

  2. In [HM15], the derivation of the shape ODE for the fine structure of the distance profile near (‘tip scaling’) was followed by asymptotic arguments to describe the jump like behavior of the velocity profiles near (‘transition scaling’) and local properties of the distance profile near (‘base scaling’). This strategy hung on explicit matching conditions and required several lengthy computations. The revised approach is more robust, and replaces the study of the advance-delay-differential equations (1.4) by the analysis of equivalent integral equations. More precisely, integrating the traveling wave equations (1.4) we get


    where the convolution is with the indicator function from (1.11). This system can also be written as


    and turns out to be very useful in the asymptotic analysis. In particular, thanks to (1.13) we do not need to introduce the transition and the base scaling anymore and satisfy the hidden matching conditions implicitly, see the discussion in §2.1.

  3. In the present paper we cover a broader class of potentials and distinguishing between the parameters and we gain a better understanding of the different contributions to the error terms. In Assumption 1 we require that both and because otherwise some of the error estimates derived below would no longer be small but could attain values of order or larger. It is not clear, at least to the authors, whether high energy waves are unstable for or , or whether our approximation and stability results are still valid but necessitate a more detailed analysis in the proofs. To ease the notation we also require that . This seems to be surprising at a first glance but gets more plausible if we study the class

    which includes both (1.8) and (1.9) as special cases. In fact, for those potentials the next-to-leading term in the expansion for is singular and regular for and , respectively, but not defined for . The criticality of shows up several times in our asymptotic analysis and manifests in algebraic error estimates involving the exponent . We emphasize that all asymptotic formulas can also be established for but entail logarithmic corrections in the error bounds, see the discussion after Lemma 5 for more details.

The non-asymptotic part of the the Friesecke Pego-Theory is reviewed in §3.1 and guarantees that the nonlinear orbital stability in the sense of Main Result 2 depends only on the spectral properties of the linearized FPU chain. More precisely, inserting the formal ansatz

into (1.2), we obtain after linearization with respect to the perturbations and the dynamical equation

where we scaled the right hand side by for convenience. The operator acts on pairs of spatial functions via


and depends explicitly of the distance profile of the underlying traveling wave. It has been observed in [FP02, FP04a] that the spectrum of the operator is -periodic due to the spatial discreteness of (1.1) and that the existence of a one-parameter family of solitary waves dictates that the kernel is a two-dimensional Jordan block, where proper and cyclic eigenfunctions represent the neutral modes related to shifts and accelerations of the wave, respectively. Moreover, the essential part of this spectrum can be computed explicitly and has strictly negative real part in exponentially weighted spaces that penalizes perturbations in front of the wave. These properties hold under fairly mild assumption on the interaction potential and are not restricted to any asymptotic regime, see the discussion in §3.2.

The main analytical task is to show that the operator does not admit any unstable point spectrum in the symplectically orthogonal complement of the neutral modes, where the underlying symplectic product is provided by the Hamiltonian structure of (1.1), see again §3.1 for the details. This crucial spectral property has been shown in [FP04b] for near sonic waves by rigorously relating the operator to the well-studied properties of solitary waves in the KdV equation and the corresponding linear differential operator. For general FPUT waves, however, no spectral theory is available for the linear advance-delay-differential operator and the stability problem is hence completely open.

Our main result concerning orbital stability is formulated in Theorem 22 and disproves for all sufficiently large the existence of unstable eigenfunctions in the aforementioned symplectic complement. Our proof strategy differs essentially from the asymptotic approach in [FP04b] because in our case is not related to the long-wave length limit and can hence no longer be regarded as a perturbation of the linear differential operator from the KdV theory. Instead we show that any appropriately rescaled eigenfunction solves locally the linearized shape ODE from §2.2 up to small error terms, and that the global behavior of any eigenfunction is completely determined by its local properties, the desired spatial decay, and certain asymptotic matching conditions. The underlying ideas are explained in §3.4 on a heuristic level while the technical details are presented in §3.5 and §3.6.

2 Family of nonlinear high-energy waves

In this section we prove the existence of a smooth family of nonlinear high-energy waves and derive explicit approximation formulas for the wave profiles as well as their derivatives with respect to the wave speed. For this analysis it is convenient to replace the large wave speed by a corresponding small quantity. In this paper we choose


because this quantity determines, as explained in §2.1, the spatial scale of the asymptotic shape ODE and behaves asymptotically like , the small parameter from the variational approach in [HM15].

2.1 Heuristic arguments and overview

As preparation for the rigorous asymptotic analysis derived below, we start with an informal overview on the high energy limit and explain on a heuristic level, why the solution to the nonlinear advance-delay-differential equations (1.4) is asymptotically determined by an ODE initial value problem. Our discussion is, as already mention in §1, based on the arguments from [HM15] but simplifies some computations since the crucial scaling laws and matching conditions are inferred from the nonlinear integral formulation in (1.12).

Distance profile near the origin. The first key observation — at least for unimodal and even wave profiles as in Figure 1.1 that attain their minimum at — is the following. If the is close to , the algebraic singularity of implies


so both the advance and the delay term in the second-order equation (1.6) can be neglected. This gives the approximate ODE

for the ‘tip of the tent’, but this equation is still singular as it involves very large quantities.

Asymptotic shape ODE. In the next step we rescale both the spatial variable and the amplitude of the distance profile by the ansatz




where the scaling parameter will be determined below. Combining this with the properties of from Assumption 1 and using


we get


and arrive at the nonsingular ODE


which can be accompanied by initial conditions. Restricting to even profiles we prescribe


and infer

as heuristic rule from (2.3).

Further scaling law. Our next goal is to compute a more explicit formula for . The integral formulation (1.13) implies

thanks to and the evenness of . In view of (2.3) and (2.6) we obtain




Inserting the asymptotic ODE (2.7) into (2.9) provides

and hence

due to , where the differential operator is defined as


and will be used frequently. By elementary ODE arguments — see Lemma 2 for precise statements — we compute both and and obtain


where we used that and are asymptotically small and large, respectively, as . A precise definition of — and hence also of and — will be given below in Lemma 2.17.

Discussion. Our formal asymptotic analysis from above can be summarized as follows. Rescaling the distance profile by (2.3), where , , and are given as in (2.1)+(2.4)+(2.12), the rescaled distance profile satisfies on the interval the ODE (2.7) along with the initial values (2.8) up to small error terms. Moreover, all these quantities are on the approximate level completely determined by and , and give rise to explicit asymptotic formulas for on the interval .

The final argument is that the local behavior of the distance profile determines also the global one because (2.2) can also be used to justify the approximate differential equations


for the ‘base of the tent’ and

for the ’tail of the tent’, which by evenness hold also for . In principle, these identities allow us — as in [HM15] — to recover the entire distance profile provided that the constants of integration will be determined such that both and are continuous at and , and vanish at . For our asymptotic proofs, however, it is much more convenient to rely on the approximation


where the indicator function of the interval serves as both convolution kernel and cut-off function. In fact, this formula bridges neatly between local and global approximations since the right hand side represents a -function on which depends only on the behavior of on . In particular, (2.14) encodes the above mentioned matching conditions in a concise and implicit way.

We finally mention that the even velocity profile can be approximated either piecewise via the approximate differential equation

or better globally by .

Outlook. In what follows we employ the outlined ideas and prove that the asymptotic ODE problem corresponding to (2.7) and (2.8) determines in fact all asymptotic properties of lattice waves in the limit . Of course, the analytical details will be more involved then the above non-rigorous computations and require explicit bounds for the error terms. In principle we could exploit ODE perturbation techniques as in [HM15], but here we proceed differently and in a more concise way. We first identify in §2.4 almost explicit formulas for the approximate solutions and conclude afterwards in §2.6 the existence of a nearby exact solution by a fixed point argument in a suitable function space. From a technical point of view, the key ingredients in our approach are the properties of the solution to the shape ODE studied in §2.2, the precise definition of the scaling parameters in §2.3, and a careful investigation of a nonlocal but linear auxiliary problem, which will be introduced in §2.5 and may be regarded as a simplified and scaled version of the linearized traveling wave equation (1.13). In this way we can also establish a local uniqueness result for even high-energy waves; this question remained open in [HM15] and could only be answered in [HM17]. Finally, in §2.7 and §2.8 we quantify the smooth parameter dependence and estimate the tail decay, respectively.

2.2 Asymptotic shape ODE

In this section provide detail information on the solution to the ODE initial value problem


which determines — via and as explained in §2.1 — the asymptotic shape of the lattice waves in the high-energy limit . We also study the linearized ODE


because the two linearly independent solutions are important building blocks for the subsequent asymptotic and spectral analysis. Notice that the function is uniquely determined by the parameter from Assumption 1 and see Figure 2.1 for numerical examples.

Figure 2.1: Top row: Numerical solution of the nonlinear shape ODE (2.15): graphs of (left), (center), and (right) for (black), (dark gray), and (light gray). Centre row, bottom row: The respective plots for the even solution and the odd solution of the linearized shape ODE (2.16).
Lemma 2 (asymptotic shape ODE).

The unique solution to the ODE initial value problem (2.15) is even, convex, strictly positive and grows asymptotically linear. Moreover, we have

where both and are asymptotically constant with


Here, the constants , depend only on and the differential operator has been defined in (2.11).


A simple phase plane analysis shows that is even and that both and are increasing for . The equations for and follow by differentiating (2.15), and by direct computations we verify the conservation law

This energy law finally implies the asymptotic conditions for and as well as the claimed error estimates. ∎

Lemma 3 (linearized shape ODE).

The solution space to the linear ODE (2.16) is spanned by an even function and an odd function , which are defined by

and comply with the normalization condition

Moreover, and are asymptotically affine and constant, respectively, with

and satisfy the decay estimates

as well as

for some constant depending on .


We directly compute that and solve (2.16). All other assertions follow from Lemma 2. ∎

We finally mention explicit expressions for , namely

which imply the non-generic property .

2.3 Scaling parameters

Motivated by the heuristic discussion in §2.1 we aim to study the wave profiles and as functions of the rescaled space variable from (2.5) in order to resolve the fine structure of the distance profile near . We therefore specify the -dependence of , , and by the following result, which does not alter the algebraic relations in (2.4) and (2.10) but provides a precise definition for in terms of . This value is asymptotically consistent with (2.12) and ensures that the approximate distance profile constructed in Lemma 5 below has sufficiently nice properties.

Figure 2.2: Top row. The scaled distance profile and the corresponding velocity profile from (2.19). The shaded region indicate the intervals and related to the ‘tip of the tent’ and the ‘base of the tent’, respectively. The explicit approximations and are constructed in Lemma 5 and Lemma 8, respectively, and vanish identically outside and , respectively. Bottom row. The resulting profiles after applying and , see Assumption 1.
Lemma 4 (parameters , , and ).

For all sufficiently large , there exists a unique such that


Moreover, we have


as well as

and depends smoothly on .


The left hand side in (2.17) is a strictly increasing and asymptotically linear function in its argument , attains a positive minimum at the origin, and involves the small prefactor . On the other hand, the right hand side is concave, strictly increasing, and independent of . It also vanishes at the origin and grows sublinearly. We therefore conclude for all sufficiently large that there exists precisely two solutions and to (2.17), which depend smoothly on and satisfy

We set and in view of the asymptotic properties of from Lemma 2 we reformulate (2.17) as

This implies

and yields the desired results by direct computations and thanks to the asymptotic formulas from Lemma 2. ∎

Our subsequent analysis relies on the ansatz


which scales the space variable but not the amplitudes, and the integral equations in (1.13) transform into




Here, is the indicator function of the interval


which corresponds to the interval with respect to the unscaled space variable , and we readily verify


i.e., the convolution kernel in (2.20) is a tent map with height and base length . For later use we also introduce a modified discrete Laplacian by


where denotes an arbitrary weight parameter.

As illustrated in Figure 2.2, the profiles and possess very large amplitudes in the limit . We therefore define


and expect — according to the informal discussion in §2.1 — that


where the limit profiles


can be computed from the unique solution of the asymptotic shape ODE (2.15).

We finally mention that the amplitude of remains bounded as while attains values of order . We decided not to normalize , because this would complicate the notations in §3, but always display normalized velocities in the plots of the numerical data, see for instance Figure 1.2.

2.4 Explicit approximation formulas

We now introduce an approximate distance profile. The existence of a nearby exact solution is proven in Theorem 7 below, which also provides explicit bounds the approximation error.

Lemma 5 (approximate distance profile).

The function