A Unifying Perspective: Solitary Traveling Waves As Discrete Breathers And Energy Criteria For Their Stability

TL;DR

This paper offers a unified framework for analyzing the spectral stability of solitary traveling waves in Hamiltonian lattices, linking eigenvalue problems with energy criteria and providing explicit stability change conditions.

Contribution

It introduces an energy-based spectral stability criterion for solitary waves in Hamiltonian lattices, connecting eigenvalues with energy dependence on wave velocity and providing explicit unstable eigenvalue calculations.

Findings

Stability change linked to energy functional monotonicity.

Explicit unstable eigenvalues near critical velocity.

Validated criterion with analytical and numerical examples.

Abstract

In this work, we provide two complementary perspectives for the (spectral) stability of solitary traveling waves in Hamiltonian nonlinear dynamical lattices, of which the Fermi-Pasta-Ulam and the Toda lattice are prototypical examples. One is as an eigenvalue problem for a stationary solution in a co-traveling frame, while the other is as a periodic orbit modulo shifts. We connect the eigenvalues of the former with the Floquet multipliers of the latter and based on this formulation derive an energy-based spectral stability criterion. It states that a sufficient (but not necessary) condition for a change in the wave stability occurs when the functional dependence of the energy (Hamiltonian) $H$ of the model on the wave velocity $c$ changes its monotonicity. Moreover, near the critical velocity where the change of stability occurs, we provide explicit leading-order computation of the unstable eigenvalues, based on the second derivative of the Hamiltonian $H"(c_0)$ evaluated at the critical velocity $c_0$. We corroborate this conclusion with a series of analytically and numerically tractable examples and discuss its parallels with a recent energy-based criterion for the stability of discrete breathers.

Equations

$$\frac{du}{dt}=p,\frac{dp}{dt}=F(u)=-\frac{\partial \mathcal{H}}{\partial u},$$

$$\frac{dU}{dt}=J\nabla \mathcal{H}(U),$$

$$U=\left(\begin{array}[]{c}{u}\\ {p}\end{array}\right),\quad J=\left(\begin{array}[]{cc}0&I\\ -I&0\end{array}\right).$$

$$u_n(t)=\hat{u}(\xi ),p_n(t)=\hat{p}(\xi ),\xi =n-ct,$$

$$\lambda Z=\mathcal{L}Z$$

$$\mathcal{L}=c{\partial }_{\xi }+J{\nabla }^2\mathcal{H}(\hat{U}),$$

$$Z=\left(\begin{array}[]{c}{W}\\ {P}\end{array}\right),\quad\hat{U}=\left(\begin{array}[]{c}{\hat{u}}\\ {\hat{p}}\end{array}\right),\quad J\nabla^{2}\mathcal{H}(\hat{U})=\left(\begin{array}[]{cc}0&I\\ F^{\prime}(\hat{u})&0\end{array}\right).$$

$${\mathcal{L}}^{\ast }=(-{\nabla }^2\mathcal{H}(\hat{U})J-c{\partial }_{\xi })=-J^{-1}\mathcal{L}J$$

$$\begin{array}{rl}0& =⟨J^{-1}{e}_0,e_1⟩=\int (-J^{-1}{\partial }_{\xi }\hat{U})\cdot ({\partial }_c\hat{U})d\xi \\ & =\int \frac{1}{c}\nabla \mathcal{H}(\hat{U})\cdot \frac{\partial \hat{U}}{\partial c}d\xi =\frac{1}{c}\int \frac{\partial \mathcal{H}(\hat{U})}{\partial c}d\xi =\frac{1}{c}{H}^{\prime }(c),\end{array}$$

$$\lambda =\pm \sqrt{\frac{H^{\prime \prime }(c_0)}{{\alpha }_1{c}_0}(c-c_0)}+O(|c-c_0|)$$

$$\begin{split}H&=\sum_{n=-\infty}^{\infty}\biggl{[}\frac{p_{n}^{2}}{2}+V(u_{n+1}-u_{n})\\ &+\sum_{m=-\infty}^{\infty}\frac{\Lambda(m)}{4}(u_{n}-u_{n+m})^{2}\biggr{]}.\end{split}$$

$$\Lambda (m)=\rho ({\mathrm{e}}^{\gamma }-1){\mathrm{e}}^{-\gamma |m|}(1-{\delta }_{m,0}),$$

$$V(u)=\frac{u^2}{2}-\frac{u^3}{3},$$

$$V(u)=\{\begin{array}{ll}\frac{u^2}{2},& |u|\le u_c\\ \frac{\chi }{2}(|u|-u_c{)}^2+u_c|u|-\frac{u_c^2}{2},& |u|>u_c,\end{array}$$

$$\epsilon (t)=1-\frac{\bar{H}(t)}{H},$$

$$\eta (t)=\frac{X(t)-X(0)}{tc}-1,$$

$$\dfrac{dY}{dt}=J_{1}\nabla{\cal H}(Y),\quad J_{1}=\left(\begin{array}[]{cc}0&S^{-1}-I\\ I-S&0\end{array}\right).$$

$$U(t;c)=\hat{U}(\xi;c)=\left(\begin{array}[]{c}{\hat{u}(\xi)}\\ {\hat{p}(\xi)}\end{array}\right),\quad\xi=n-ct,$$

$$U(\xi,t)=\left(\begin{array}[]{c}{\hat{u}(\xi)}\\ {\hat{p}(\xi)}\end{array}\right)+\epsilon e^{\lambda t}\left(\begin{array}[]{c}{W(\xi)}\\ {P(\xi)}\end{array}\right)=\hat{U}(\xi)+\epsilon e^{\lambda t}Z(\xi)$$

$$\begin{array}{rl}⟨e_2,J^{-1}{e}_0⟩& =⟨e_2,J^{-1}\mathcal{L}{e}_1⟩=⟨J^{-1}\mathcal{L}{e}_2,e_1⟩\\ & =⟨J^{-1}{e}_1,e_1⟩=0,\end{array}$$

$${\alpha }_1=⟨e_0,J^{-1}{e}_3⟩=-⟨e_1,J^{-1}{e}_2⟩\ne 0.$$

$${\mathcal{L}}_0^3{e}_3={\mathcal{L}}_0^2{e}_2={\mathcal{L}}_0{e}_1=e_0=-{\partial }_{\xi }\hat{U}(\xi ;c_0).$$

$$K_{jk}=⟨J^{-1}{e}_j,{\mathcal{L}}_1{e}_k⟩,L_{jk}=⟨J^{-1}{e}_j,{\mathcal{L}}_2{e}_k⟩.$$

$$A=\left(\begin{array}[]{cccc}0&0&0&0\\ 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\end{array}\right)$$

$$\mathcal{L}({\partial }_{\xi }\hat{U}(\xi ;c))=0$$

$$\mathcal{L}({\partial }_c\hat{U}(\xi ;c))=-({\partial }_{\xi }\hat{U}(\xi ;c)).$$

$$0$$

$${\mathcal{L}}_0{U}_1$$

$$2{\mathcal{L}}_0{U}_2+{\mathcal{L}}_1{U}_1$$

$$\begin{array}{rl}& H^{\prime \prime }(c_0)=\nabla H\cdot 2U_2+U_1\cdot {\nabla }^{2}H{U}_1\\ & =\int [c_0{J}^{-1}(e_0)\cdot 2U_2+U_1\cdot J^{-1}({\mathcal{L}}_0-c_0{\partial }_{\xi })U_1]d\xi \\ & =\int [c_0{J}^{-1}{\mathcal{L}}_0{U}_1\cdot 2U_2+U_1\cdot (J^{-1}{\mathcal{L}}_0{U}_1-c_0{J}^{-1}{\partial }_{\xi }{U}_1)]d\xi \\ & =\int U_1\cdot (J^{-1}{e}_0)d\xi +c_0\int U_1\cdot J^{-1}({\mathcal{L}}_{0}2U_2-{\partial }_{\xi }{U}_1)d\xi \\ & =c_0\int U_1\cdot J^{-1}(-2{\partial }_{\xi }{U}_1-{\mathcal{L}}_1{U}_1)d\xi \\ & =-c_0\int {\mathcal{L}}_0{e}_2\cdot J^{-1}2{\partial }_{\xi }{U}_{1}d\xi -c_0\int e_1\cdot J^{-1}{\mathcal{L}}_1{e}_{1}d\xi \\ & =-c_0\int e_2\cdot (-1)J^{-1}{\mathcal{L}}_{0}2{\partial }_{\xi }{U}_{1}d\xi -c_0\int e_1\cdot J^{-1}{\mathcal{L}}_1{e}_{1}d\xi \\ & =c_0\int e_2\cdot J^{-1}{\mathcal{L}}_{1}2e_{0}d\xi -c_0\int e_1\cdot J^{-1}{\mathcal{L}}_1{e}_{1}d\xi \\ & =-c_0(2K_{20}-K_{11}).\end{array}$$

Full text

I Introduction

Solitary traveling waves (STWs) are ubiquitous in Hamiltonian lattice dynamical systems with intersite interactions. They arise in the model at the very foundation of nonlinear science, namely the Fermi-Pasta-Ulam (FPU) lattice FPU , as well as in the Toda lattice toda , one of the key systems of interacting particles, and, arguably, the most significant integrable one. In addition to their theoretical relevance in the above models, they constitute the most generic, robust and often experimentally tractable excitation in nonlinear systems, in particular, in granular crystals nesterenko ; sen ; review and other materials.

Given the relevance of STWs in theoretical, numerical mertens ; eilbeck ; english and experimental nesterenko ; sen studies, it is natural to be concerned about their stability. This may be accessible in some special cases, such as the Toda lattice wayne , or the FPU problem in the low-energy (near-sonic) regime pegof3 ; pegof4 , where specialized techniques become available due to the system’s integrability (or proximity to it). Nevertheless, from a physical perspective, it would be desirable to have a more general criterion that would be intuitive as well as straightforward to test. This is especially important given that in a number of studies mertens1 ; mertens2 ; at , the possibility of unstable STWs has been demonstrated.

In the present work, we offer such a criterion (a sufficient yet not necessary condition) by establishing that a change in the monotonicity of the STW’s energy (Hamiltonian $H$) dependence on the velocity $c$ will result in a change in its (spectral) stability. In other words, we establish that when, for a critical velocity $c_{0}$, it happens that $H^{\prime}(c_{0})=0$, a pair of eigenvalues associated with the traveling wave vanish, entailing the potential for instability. While this criterion first appeared in pegof3 , where it was motivated by the study of the FPU problem in the near-sonic limit, here we provide both a concise proof, and also a definitive leading-order calculation for these two near-zero eigenvalues to explicitly show why (and when) instability appears. We also systematically test the criterion numerically in a broad array of physically relevant cases.

Equally important in our approach is the fact that we provide a generalized perspective of the problem of the stability of STWs in a Hamiltonian lattice. In the frame traveling with the solution, the stability leads to a standard eigenvalue problem. Yet, here, motivated by earlier works such as floria , we also propose a complementary approach, where the solution is viewed as a periodic orbit of the map involving (a) running the solution for a period of $h/c$, where $h$ is the lattice spacing, rescaled to unity below and (b) shifting back by one lattice site. In light of this periodicity, Floquet analysis can be brought to bear and will turn out to yield coincident stability conclusions about instabilities produced by the criterion put forth. Furthermore, this perspective enables a unification of the lattice STWs in such Hamiltonian systems through their consideration as discrete breathers. Here the effective frequency $\omega$ is proportional to their velocity $c$ according to $\omega=2\pi c/h$. This, in turn, directly connects the criterion we analyze with a recently established criterion for the spectral stability of discrete breathers dmp . We emphasize here that the unifying connection of STWs with breathers does not impose any a-priori restrictions on the nature of their decay of at infinity.

The paper is organized as follows. In Sec. II we formulate the problem, analyze the properties of the linear operator associated with a STW and prove the energy-based stability criterion. We also describe the behavior of the relevant eigenvalues near the critical velocity, based on the derivation presented in Appendix A. In Sec. III we discuss an alternative perspective for the spectral stability, which is associated with the Floquet analysis. Our results are corroborated by numerical examples in Sec. IV, with further details provided in Appendix B. We summarize our findings and discuss some open questions in Sec. V.

II Stability Analysis in the Co-traveling Frame and the

Energy Criterion

We consider a rescaled Hamiltonian system of the form

[TABLE]

where ${\cal H}$ denotes the Hamiltonian energy density of the system, $u(t)$ and $p(t)$ are infinite-dimensional vectors denoting the displacement and particle velocity values on the lattice, with components $u_{n}$ and $p_{n}$, respectively. In a more compact notation, Eq. (1) can be written as

[TABLE]

where

[TABLE]

We assume the existence of STWs for a continuous interval of velocities. These are localized solutions of the form

[TABLE]

where $c$ denotes the velocity of the wave, and $\xi$ is the co-traveling frame variable (note that $\hat{p}(\xi)=-c\hat{u}^{\prime}(\xi)$), with finite energy (see Appendix A for more details). Linearization about the STW in the co-traveling frame, with $u(\xi,t)=\hat{u}(\xi)+\epsilon e^{\lambda t}W(\xi)$ and $p(\xi,t)=\hat{p}(\xi)+\epsilon e^{\lambda t}P(\xi)$ for small $\epsilon$, then yields the eigenvalue problem

[TABLE]

for the linear operator

[TABLE]

where

[TABLE]

Solving the problem in Eq. (3) provides information about the stability of the STW, through the spectrum of the linearization operator $\mathcal{L}$, with adjoint

[TABLE]

(note that $J\mathcal{L}$ is self-adjoint). Given the time translation symmetry, an important feature of $\mathcal{L}$ is the existence of an eigenvector $e_{0}=-\partial_{\xi}\hat{U}$ associated with eigenvalue $\lambda=0$. The corresponding generalized eigenfunction is $e_{1}=\partial_{c}\hat{U}$, i.e., $\mathcal{L}e_{1}=e_{0}$. In other words, the spectrum of $\mathcal{L}$ always contains a double eigenvalue at zero. Moreover, by symmetry, the algebraic multiplicity of the zero eigenvalue can only be even.

The presence of an additional instability presupposes the increase of the algebraic multiplicity of the [math] eigenvalue. Since the kernel of $\mathcal{L}$ is one-dimensional, an algebraic multiplicity higher than two (i.e., at least four) implies that there exists $e_{2}$ such that $\mathcal{L}e_{2}=e_{1}=\partial_{c}\hat{U}$. Since $J^{-1}e_{0}=-J^{-1}\partial_{\xi}\hat{U}$ is in the kernel of $\mathcal{L}^{*}$, this yields the solvability condition

[TABLE]

where $H=\int{\cal H}d\xi$ is the conserved Hamiltonian of the system, and $\langle\cdot,\cdot\rangle$ denotes the relevant inner product.

As soon as $c$ deviates from the critical velocity $c_{0}$ satisfying $H^{\prime}(c_{0})=0$, the above solvability condition fails (e.g. assuming $H^{\prime\prime}(c_{0})\neq 0$), and hence two eigenvalues start to move away from zero and can possibly emerge on the real axis. Thus the condition $H^{\prime}(c_{0})=0$ constitutes a threshold for instability of STWs, as per the concise proof above and detailed numerical considerations below extending the formulation of pegof3 . In fact, by computing the leading-order approximation of these two near-zero eigenvalues near $c_{0}$ one can reveal the trend of their motion. Suppose, as will be typically the case when the stability changes (including examples in Sec. IV below), that the generalized kernel of $\mathcal{L}$ is exactly four-dimensional at $c_{0}$, with $\mathcal{L}^{3}e_{3}=\mathcal{L}^{2}e_{2}=\mathcal{L}e_{1}=e_{0}=-\partial_{\xi}\hat{U}$. Then, as shown in Appendix A, the pair of eigenvalues of $\mathcal{L}$ responsible for the change of stability will be given by

[TABLE]

for $c$ is near $c_{0}$, where nonzero $\alpha_{1}$ is defined in (12) in terms of generalized eigenvectors.

In Sec. IV we numerically verify the theoretical predictions (and test the validity of Eq. (6)), showing that a change of the monotonicity of $H(c)$ will constitute a sufficient (but not necessary) condition for the transition from stability to instability, or vice versa, depending on the sign of $H^{\prime\prime}(c_{0})\alpha_{1}c_{0}$.

III A Complementary Perspective: Floquet Analysis of the

time $T=h/c$ map

Let us now envision anew the case of a STW on a lattice. Over the period $T=h/c$ (below we again set $h=1$), the STW $\hat{U}$ moves over by one lattice site. However, due to the integer shift invariance of the lattice, the configuration has to be identical to the one with which we started. This means that upon running for a period and shifting back using the shift operator $S$ such that $Su_{n}(t)=u_{n-1}(t)$, we generate a periodic orbit on the lattice floria . Thus, a fixed point of this operation consisting of (a) run for $T=1/c$ and (b) shift, is a discrete breather (DB) i.e., a localized time-periodic solution aubry ; FlachPR2008 by construction with frequency $\omega=2\pi c$. Yet, at the same time the resulting profile constitutes a lattice STW.

Two important consequences of this complementary perspective are as follows. (1) The fixed point operation discussed above has a corresponding monodromy matrix aubry ; FlachPR2008 ; arnold whose eigenvalues are the Floquet multipliers (FMs) of the relevant periodic orbit. These FMs determine the stability of the periodic orbit (i.e., in this case of the STW), as do the eigenvalues of co-traveling problem computation. Hence, one should expect that an instability manifested through an eigenvalue crossing zero should be accompanied by a FM $\mu$ crossing unity, due to the well known relation $\mu=e^{\lambda T}$ between the multipliers and eigenvalues ourjesus . (2) Given the intimate connection of lattice STWs and DBs, an immediate correlation emerges between the criteria for stability change of discrete breathers, such as $H^{\prime}(\omega)=0$ that was recently established in dmp and the stability of lattice STWs discussed here (and also in pegof3 ). Observing that for lattice STWs, $\omega=2\pi c$, an alternative derivation of the latter from the former is, in fact, immediate.

IV Numerical Corroboration

We now test the above prediction in a set of numerical examples with the generalized Hamiltonian of the form

[TABLE]

Here $V(u)$ is a generic potential governing the nonlinear interactions between nearest neighbors, and $\Lambda(m)$ are the coefficients of all-to-all linear long-range interactions, which decay as $|m|\rightarrow\infty$; in the absence of such interactions, $\Lambda(m)=0$. For instance,

[TABLE]

with $\rho>0$ and $\gamma>0$, corresponds to the Kac-Baker interactions, and $\Lambda(m)=\rho|m|^{-s}(1-\delta_{m,0})$ with $s=5$ ($s=3$) corresponds to the dipole-dipole (Coulomb) interactions between charged particles on a lattice. In principle, the methodology can capture nonlinear long-range interactions, but here we consider linear ones for simplicity.

As our first example, we consider the analytically tractable and well known case of the Toda lattice toda where $V(u)=e^{-u}+u-1$ while $\Lambda(m)=0$, which has a one-soliton solution of the form $u_{n}(t)=\log\left[\cosh(\kappa(n-ct-1))\mathrm{sech}(\kappa(n-ct))\right]$, where $\kappa$ is the unique positive solution of $c\kappa=\sinh(\kappa)$. The resulting Hamiltonian can be computed explicitly for the single soliton family: $H=\sinh(2\kappa)-2\kappa$, leading to $H^{\prime}(c)=2(\cosh(2\kappa)-1)\partial_{c}\kappa>0$, resulting in generically (spectrally) stable solitary waves in the Toda lattice. This is also in tune with the nonlinear stability of the solitary waves in this case, which has been explored, e.g., in mizu .

A second famous example consists of the $\alpha$-FPU case FPU , where

[TABLE]

while $\Lambda(m)=0$. In this case too, as identified via the methods of mertens ; eilbeck ; english ; jphysa (see Appendix B, for details on numerical simulations) and shown in Fig. 1, the family of STWs numerically features $H^{\prime}(c)>0$, in full agreement with their identification as stable. Similar conclusions hold for the highly experimentally relevant solitary waves of granular crystals nesterenko ; sen ; review .

Arguably, these cases, while interesting from the prototypically nonlinear and experimental perspectives, are perhaps somewhat less exciting from the point of view of our criterion as they do not feature a stability change. Hence, we turn to some examples which, while more exotic from the point of view of practical applications, have been argued to be of interest and, additionally, feature a change of stability, which is especially relevant in the context of this work. The first such case that we will consider concerns the Kac-Baker interactions that have been argued to be of relevance for modeling Coulomb interactions in DNA molecules in mertens2 . In this case, we maintain the potential in Eq. (9) of the FPU case, but add long-range interactions with the kernel in Eq. (8). Fig. 2 showcases the power of the stability criterion and illustrates the complementary nature of the co-traveling steady state and the periodic orbit FM calculation approaches. It can be seen that $H^{\prime}(c)$ becomes negative (the top panel of Fig. 2) for $1.6709<c<1.6937$, for our chosen values of $\gamma=0.17$, $\rho=0.0172$ selected in tune with mertens2 . For this very interval of velocities, an eigenvalue of the operator $\mathcal{L}$ crosses through $\lambda=0$ and acquires a positive real part (dots in the bottom panel of Fig. 2). In fact, it can be shown pegof3 that the stability problem in the co-traveling frame also possesses eigenvalues $\lambda+i(2\pi jc)$, where $j\in\mathbb{Z}$. Finally, the solid curve in the bottom panel of the Fig. 2 showcases the FM calculation associated with the time $T=1/c$ map of the corresponding periodic orbit, transformed (in order to compare with the steady state eigenvalue approach) according to the relation $\lambda=\log(\mu)/T$. Confirming the complementary picture put forth, we find that in this case a FM pair crosses through $(1,0)$ and into the real axis for the exact same parametric interval.

To connect with the theoretical analysis of Eq. (6), the inset of Fig. 2 shows the dependence of $\lambda^{2}$ with respect to $c-c_{0}$, which, according to Eq. (6), must be linear in the vicinity of $c_{0}\approx 1.6937$ with the slope $\beta=H^{\prime\prime}(c_{0})/(\alpha_{1}c_{0})=-2.9794\times 10^{-4}$. Our numerical calculations yield $\beta=-3.0383\times 10^{-4}$; the mismatch of $\sim 2\%$ is likely due to the fact that $\alpha_{1}$ in Eq. (6) cannot be computed at the precise value of $c_{0}$ in the numerical setup. A similar agreement was also found in the vicinity of the other critical point at $c_{0}\approx 1.6709$.

As our final example, it is interesting to explore a case where the relevant theory does not directly apply due to limited regularity. As such an example, we consider an FPU model with the potential of the form

[TABLE]

which allows construction of explicit solitary waves at , and $\Lambda(m)=0$; here $\chi>1$ and $u_{c}>0$. In this case the potential possesses only one continuous derivative, and hence the calculation of eigenvalues $\lambda$ and FMs $\mu$ is less straightforward to justify, given the relevant jump discontinuities. Nevertheless our detailed computations, in line with the numerical results and stability conjecture in at , are in a clear agreement with the criterion put forth analytically in this work. Namely, $H^{\prime}(c)>0$ in this case too corresponds to dynamical stability, while $H^{\prime}(c)<0$ leads to the manifestation of instability.

In order to qualitatively measure the instability, we have defined two diagnostic quantities. The first of them is the energy dispersion, given by

[TABLE]

where $\bar{H}(t)$ is the energy at the nine central sites of the STW. In the case of a stable propagating wave, $\varepsilon(t)\sim 10^{-4}$. The other quantity is the relative velocity change defined as

[TABLE]

with $X(t)$ being the energy center of the STW.

The top panel of Fig. 3 shows the curve $H(c)$ for the FPU model with the potential of Eq. (10) and parameters $\chi=4$ and $u_{c}=1$; the bottom panels of this figure display the dependence of $\varepsilon_{\infty}\equiv\varepsilon(2000T)$ and $\eta_{\infty}\equiv\eta(2000T)$ with respect to $c$. In accordance with our stability criterion, $\varepsilon_{\infty}\sim 10^{-4}$ in the region for which $H^{\prime}(c)>0$, confirming a stable propagation. In the region with $H^{\prime}(c)<0$ there are three intervals of high dispersion, as measured by corresponding values of $\varepsilon$, and two intervals where the dispersion drops to low values. The region of low dispersion corresponds to STWs whose velocity is higher than the initial one (indeed, higher than the critical one and hence reverting to the stable propagation regime). Fig. 4 shows the evolution of unstable STWs in two cases, corresponding to high ($c=1.025$) and low ($c=1.034$) dispersion. In the former, linear waves are continuously being created and the STW degrades with time; in the latter, a linear wave is expelled from the STW, which transforms into a wave with a different (now in the stable regime of $c>c_{0}$) velocity. Note that in addition to demonstrating instability of waves with $c<c_{0}$, these results suggest the potential bistability between dispersive waveforms and STWs with $c>c_{0}$.

V Conclusions and Future Challenges

In summary, in this work we have presented a unified perspective connecting the stability of lattice solitary traveling waves with that of discrete breathers of an appropriate map involving running for the time associated with moving by one lattice site and shifting back. We have also concisely established a (sufficient but not necessary) criterion for the change in spectral stability of the Hamiltonian lattice STWs that seems to be in very good agreement with numerical observations and to constitute a natural extension of a criterion recently put forth for the spectral stability of discrete breathers. The specific eigenvalue responsible for the instability was theoretically identified and favorably compared to detailed numerical computations.

Nevertheless, there are numerous problems that remain open for future consideration. One relevant issue concerns the fact that the FM computation leads to as many multipliers as lattice points, while the computation of eigenvalues for a STW involves a partial differential equation (PDE). While the latter will capture the lattice instabilities, it may also feature instabilities absent on the lattice, which are a by-product of this PDE’s ability to resolve scales smaller than $h$. Hence, a more systematic connection between the spectra of the two problems (and of the instabilities that each may feature) is of paramount importance. Observe also that while this work dealt with families of STWs parameterized by velocity, in some cases such entities occur for isolated velocity values jcm ; annav , potentially being members of a wider family encompassing waveforms with non-vanishing tails. It would be interesting to explore whether our considerations can be extended to such cases. Another question is that of going to the continuum limit: our proof did not directly use the underlying lattice nature of the system (only its time reversal invariance). On the other hand, in the continuum limit, symmetries (like Galilean or Lorentz invariance) may arise. Future work will involve reconciling these two features in a consistent continuum limit picture, as well as connecting our criterion with well-established existing stability criteria, such as vk ; gss ; bara , in continuum systems. Finally, analysis of the stability of lattice STWs in systems with limited regularity, such as our last example, also merits future consideration.

Acknowledgements.

J.C.-M. thanks financial support from MAT2016-79866-R project (AEI/FEDER, UE). A.V. acknowledges support by the U.S. National Science Foundation through the grant DMS-1506904. P.G.K. gratefully acknowledges support from the Alexander von Humboldt Foundation, the Greek Diaspora Fellowship Program, the US-NSF under grant PHY-1602994. as well as the ERC under FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096).

Appendix A Proof of the leading-order approximation of the near-zero eigenvalues

In this Appendix we prove Eq. (6) in Sec. II, which provides the leading-order approximation of the eigenvalues splitting away from zero at velocities near the critical value $c_{0}$.

First, we observe that while we consider a lattice Hamiltonian system in the displacement form (2), the problem can be alternatively formulated in terms of strain variables $y_{n}(t)=(S^{-1}-I)u_{n}(t)=u_{n+1}(t)-u_{n}(t)$, where we recall from Sec. III that $S$ denotes the shift operator such that $Su_{n}(t)=u_{n-1}(t)$. If the Hamiltonian energy density can be written as $\mathcal{H}(y,p,t)$, we have, for $Y=(y,p)^{T}$,

[TABLE]

In what follows, we focus on the formulation (2), but our arguments also work for Eq. (11).

Suppose Eq. (2) has a family of solitary traveling-wave solutions ${U}(t;c)$ parametrized by the velocity $c$ taking values in some continuous interval. Then

[TABLE]

where $\xi$ is the co-traveling frame variable and $\hat{p}(\xi;c)=-c\partial_{\xi}\hat{u}(\xi;c)$. Considering the ansatz

[TABLE]

with small $\epsilon$ and linearizing around the traveling wave $\hat{U}$, we obtain Eq. (3), where the operator $\mathcal{L}$ and its adjoint $\mathcal{L}^{*}$ are given by Eq. (4) and Eq. (5), respectively.

Suppose $\hat{U}\in H^{1}(\mathbb{R}^{2})$, so that its partial derivatives in $\xi$ and $c$ are in $L^{2}(\mathbb{R}^{2})$. While this assumption implies that the displacements are localized, for problems with kink-type traveling waves in terms of displacement that tend to nonzero constant limits at infinity, we can use the strain formulation (11), in which case we assume that the traveling wave solution $\hat{Y}\in H^{1}(\mathbb{R}^{2})$, i.e., the strains are localized. One can show that the operator $\mathcal{L}$ is densely defined on $L^{2}(\mathbb{R}^{2})$. By differentiating Eq. (2) in $\xi$ and $c$, respectively, we find that $\mathcal{L}e_{0}=0$ and $\mathcal{L}e_{1}=e_{0}$, where $e_{0}=-\partial_{\xi}\hat{U}$ and $e_{1}=\partial_{c}\hat{U}$ (or, more generally, $e_{1}=\partial_{c}\hat{U}+d_{10}e_{0}$, where $d_{10}$ is any constant), implying that the algebraic multiplicity of the eigenvalue $\lambda=0$ for $\mathcal{L}$ is at least two. Let $c_{0}$ denote the critical velocity such that $H^{\prime}(c_{0})=0$. Then $\langle e_{1},J^{-1}e_{0}\rangle=0$ at this critical value, and there exists $e_{2}$ such that $\mathcal{L}e_{2}=e_{1}$. Since

[TABLE]

we have $e_{2}\in(\ker(\mathcal{L^{*}}))^{\perp}={\rm im}(\mathcal{L})$, so $e_{2}$ belongs to the range of $\mathcal{L}$, and hence there exists $e_{3}$ such that $\mathcal{L}e_{3}=e_{2}$. Assuming that the zero eigenvalue of $\mathcal{L}$ at $c_{0}$ is exactly quadruple, which is the generic case for traveling waves in Hamiltonian lattices due to symmetry, we have

[TABLE]

We now consider a neighborhood of the critical speed $c=c_{0}$ where the derivative $H^{\prime}(c)$ changes its sign. Assuming that $\hat{U}(\xi;c)$ is sufficiently smooth in $c$ near $c=c_{0}$, we have the expansion $\hat{U}(\xi;c_{0}+\epsilon)=U_{0}+\epsilon U_{1}+\epsilon^{2}U_{2}+\dots$ for small enough $\epsilon$, where $U_{0}=\hat{U}(\xi;c_{0})$, $U_{1}=(\partial_{c}\hat{U}(\xi;c))|_{c=c_{0}}$ and $U_{2}=\frac{1}{2}(\partial_{cc}\hat{U}(\xi;c))|_{c=c_{0}}$. Accordingly, the operator $\mathcal{L}$ at $c=c_{0}+\epsilon$ can be written as $\mathcal{L}=\mathcal{L}_{0}+\epsilon\mathcal{L}_{1}+\epsilon^{2}\mathcal{L}_{2}+\dots$. Let $\{e_{0},e_{1},e_{2},e_{3}\}$ be the eigenfunction and generalized eigenfunctions of $\mathcal{L}_{0}$ for $\lambda=0$ such that

[TABLE]

We then define the following constants:

[TABLE]

Remark 1

If the generalized kernel of $\mathcal{L}_{0}$ is exactly four-dimensional, then only the cases $\lambda\sim\epsilon^{1/2}$ and $\lambda\sim\epsilon$ are possible.

Indeed, this follows from the fact that two of the four eigenvalues of $\mathcal{L}$ are always zero. It suffices to calculate the leading-order terms of the eigenvalues for the perturbed operator $\mathcal{L}$ at $c=c_{0}+\epsilon$. By restricting the operator in the invariant subspace $G_{4}=\text{cl}({\rm span}\{e_{0},e_{1},e_{2},e_{3}\})$, the question reduces to the perturbation of the matrix

[TABLE]

with two constraints that hold for any $c$,

[TABLE]

and

[TABLE]

Note that the characteristic polynomial of the unperturbed matrix $A$ is $\lambda^{4}=0$. For the matrix $A$ with $O(\epsilon)$ perturbation, the characteristic polynomial is $\lambda^{4}+a_{3}\lambda^{3}+a_{2}\lambda^{2}+a_{1}\lambda+a_{0}=0$ where the coefficients $a_{j}$ are at most $O(\epsilon)$. Moreover, due to two existing constraints in Eq. (14) and Eq. (15), two of the eigenvalues are always zero, so we have $\lambda^{2}(\lambda^{2}+a_{3}\lambda+a_{2})=0$. Thus, either $\lambda\sim\epsilon^{1/2}$ (if $a_{2}\neq 0$) or $\lambda\sim\epsilon$ (if $a_{2}=0$).

Here we focus on the case $\lambda\sim\epsilon^{1/2}$ and show below that it requires $H^{\prime\prime}(c_{0})\neq 0$. Since Eq. (14) holds for any $c$, direct calculation shows that

[TABLE]

Moreover, utilizing the fact that (15) is true for any $c$, one can expand both sides in $\epsilon$ and obtain

[TABLE]

Since $H^{\prime}(c_{0})=0$, we can write $H(c)=H(c_{0})+\frac{\epsilon^{2}}{2}H^{\prime\prime}(c_{0})+o(\epsilon^{2})$, where

[TABLE]

Assuming $\lambda=\epsilon^{1/2}\lambda_{1}+\epsilon\lambda_{2}+\epsilon^{3/2}\lambda_{3}+\dots$ and $Z=Z_{0}+\epsilon^{1/2}Z_{1}+\epsilon Z_{2}+\epsilon^{3/2}Z_{3}+\dots$ and substituting these into Eq. (3), we obtain

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

From Eq. (22), we find that $Z_{0}=-\partial_{\xi}\hat{U}(\xi;c_{0})=e_{0}$. Then Eq. (23) suggests that $Z_{1}=\lambda_{1}e_{1}+d_{10}e_{0}$, where $d_{10}$ is a constant. Note that $Z_{2}$ and $Z_{3}$ can be written as

[TABLE]

where $Z_{2}^{\perp}$ and $Z_{3}^{\perp}$ are in $G_{4}^{\perp}$, and $d_{2j}$, $d_{3j}$, $j=0,\dots,3$ are constants. Projecting Eq. (24) onto $J^{-1}e_{0}$ yields $\lambda_{1}\langle J^{-1}e_{0},Z_{1}\rangle+\lambda_{2}\langle J^{-1}e_{0},Z_{0}\rangle=K_{00}$. The left hand side is zero since $H^{\prime}(c_{0})=0$, and one can show that the right hand side vanishes ($K_{00}=0$) upon considering Eq. (17). Projecting Eq. (24) onto $J^{-1}e_{1}$ and recalling Eq. (12), we obtain $d_{23}\langle J^{-1}e_{1},e_{2}\rangle+K_{10}=d_{23}\alpha_{1}+K_{10}=0$, so

[TABLE]

Projecting Eq. (24) onto $J^{-1}e_{2}$ and using (12), we find that $d_{22}\langle J^{-1}e_{2},e_{1}\rangle+K_{20}=-d_{22}\alpha_{1}+K_{20}=-\lambda_{1}^{2}\alpha_{1}$, and thus

[TABLE]

Projecting Eq. (24) onto $J^{-1}e_{3}$, we have

[TABLE]

where we used Eq. (12) and set $\alpha_{2}=\langle J^{-1}e_{3},e_{2}\rangle$. Projecting Eq. (25) onto $J^{-1}e_{0}$ yields $\lambda_{1}d_{23}\langle J^{-1}e_{0},e_{3}\rangle=-\lambda_{1}d_{23}\alpha_{1}=\lambda_{1}K_{01}$, which again yields Eq. (27) since $K_{01}=K_{10}$. Projecting Eq. (25) onto $J^{-1}e_{1}$, we obtain

[TABLE]

Finally, projection of Eq. (26) onto $J^{-1}e_{0}$ yields

[TABLE]

Using the equations (27), (28), (29), (30), (31) along with the fact that $K$ is symmetric, we obtain

[TABLE]

Since two eigenvalues are always zero, this equation should have two zero roots. This implies

[TABLE]

which can also be shown directly using projections of Eq. (17) onto $J^{-1}e_{j}$, $j=0,\dots,3$, and projection of Eq. (18) onto $J^{-1}e_{0}$. We then obtain

[TABLE]

Thus, for $\lambda\sim\epsilon^{1/2}$ it is necessary to have $H^{\prime\prime}(c_{0})\neq 0$, and the behavior of the two eigenvalues splitting away from zero at $c\neq c_{0}$ is described by Eq. (6) in Sec. II.

Appendix B Numerical methods for computing solitary traveling waves

In this Appendix, we describe the numerical procedures we used to compute solitary waves in a lattice with the Hamiltonian in Eq. (7) and analyze their stability. The governing equations corresponding to Eq. (7) are

[TABLE]

where the overdots here and in what follows denote the time derivatives. Since the solitary solutions we consider are kink-like in terms of displacement, it is more convenient to rewrite Eq. (32) in terms of the strain variables $y_{n}=u_{n+1}-u_{n}$, obtaining

[TABLE]

To find solitary traveling wave solutions, we use the procedure followed in Yasuda . To this end, we seek solutions of Eq. (33) in the co-traveling frame corresponding to velocity $c$:

[TABLE]

obtaining the advance-delay partial differential equation

[TABLE]

Traveling waves $\phi(\xi)$ are stationary solutions of Eq. (34). They satisfy the advance-delay differential equation

[TABLE]

Solitary traveling waves are solutions that in addition satisfy

[TABLE]

Following the approach in eilbeck , we assume that $\phi(\xi)=o(1/\xi)$ and $\phi^{\prime}(\xi)=o(1/\xi^{2})$ as $|\xi|\rightarrow\infty$, multiply Eq. (35) by $\xi^{2}$ and integrate by parts to derive the identity

[TABLE]

which imposes the constraint (36) on the traveling wave solutions. Here we assume that $\Lambda(m)$ decays faster than $1/m^{3}$ at infinity, so that the series on the left hand side converges.

To solve Eq. (35) numerically, we introduce a discrete mesh with step $\Delta\xi$, where $1/\Delta\xi$ is an integer, so that the advance and delay terms $\phi(\xi\pm m)$ are well defined on the mesh. We then use a Fourier spectral collocation method for the resulting system with periodic boundary conditions Trefethen with large period $L$. Implementation of this method requires an even number $\mathcal{N}$ of collocation points $\xi_{j}\equiv j\Delta\xi$, with $j=-\mathcal{N}/2+1,\ldots,\mathcal{N}/2$, yielding a system for $\xi$ in the domain $(L/2,L/2]$, with $L=\mathcal{N}\Delta\xi$ being an even number, and the long-range interactions are appropriately truncated. To ensure that the solutions satisfy Eq. (36), we additionally impose a trapezoidal approximation of Eq. (37) on the truncated interval at the collocation points. This procedure is independent of the potential and the interaction range. However, the choices of $\Delta\xi$ and $L$ depend on the nature of the problem. In the particular cases considered in the paper, we used $\Delta\xi=0.1$, $L=800$ for the $\alpha$-FPU lattice with nearest-neighbor potential in Eq. (9) and Kac-Baker long-range interactions with coefficients in Eq. (8) and $\Delta\xi=0.025$, $L=200$ for the FPU lattice with piecewise quadratic short-range interaction potential in Eq. (10) and no long-range interactions.

To investigate spectral stability of an obtained traveling wave $\phi(\xi)$, we substitute

[TABLE]

into Eq. (34) and consider $O(\epsilon)$ terms resulting from this perturbation. This yields the following quadratic eigenvalue problem:

[TABLE]

By defining $b(\xi)=\lambda a(\xi)$, we transform this equation into the regular eigenvalue problem

[TABLE]

for the corresponding linear advance-delay differential operator $\mathcal{M}$. Note that this problem is equivalent to the eigenvalue problem (3) via the transformation $(a(\xi),b(\xi))=(W(\xi),P(\xi)+cW^{\prime}(\xi))$. Spectral stability can be determined by analyzing the spectrum of the operator $\mathcal{M}$ after discretizing the eigenvalue problem the same way as the nonlinear Eq. (35) and again using periodic boundary conditions. A solution is stable when the spectrum contains no real eigenvalues.

An alternative method for determining the stability of the traveling waves is to use Floquet analysis. To this end, we cast traveling waves $\phi(\xi)$ as fixed points of the map

[TABLE]

which is periodic modulo shift by one lattice point, with period $T=1/c$. Indeed, one easily checks that $\hat{y}_{n}(t)=\phi(n-ct)=\phi(n-t/T)$ satisfies $\hat{y}_{n+1}(T)=\hat{y}_{n}(0)=\phi(n)$ and $\dot{\hat{y}}_{n+1}(T)=\dot{\hat{y}}_{n}(0)=-c\phi^{\prime}(n)$. To apply the Floquet analysis, we trace time evolution of a small perturbation $\epsilon w_{n}(t)$ of the periodic-modulo-shift (traveling wave) solution $\hat{r}_{n}(t)$. This perturbation is introduced in Eq. (33) via $y_{n}(t)=\hat{y}_{n}(t)+\epsilon w_{n}(t)$. The resulting $O(\epsilon)$ equation reads

[TABLE]

Then, in the framework of Floquet analysis, the stability properties of periodic orbits are resolved by diagonalizing the monodromy matrix $\mathcal{F}$ (representation of the Floquet operator for finite systems), which is defined as:

[TABLE]

For the symplectic Hamiltonian systems considered in this work, the linear stability of the solutions requires that the monodromy eigenvalues $\mu$ (also called Floquet multipliers) lie on the unit circle. The Floquet multipliers can thus be written as $\mu=\exp(i\theta)$, with Floquet exponent $\theta$.

Note that the two procedures for analyzing spectral stability described above require the potential $V(u)$ to be twice differentiable, as in the case of the $\alpha$-FPU problem considered in Sec. IV. Due to the absence of such regularity in the case of the piecewise quadratic potential in Eq. (10), the examination of stability was performed solely on the basis of direct numerical simulations. Specifically, it was analyzed by means of tracking the dynamics of a slightly perturbed solution $[\{\hat{y}_{n}(0)\},\{\dot{\hat{y}}_{n}(0)\}]$. To this aim, the fourth order explicit and symplectic Runge-Kutta-Nyström method developed in Calvo , with time step equal to $10^{-3}$, was used.