## Abstract

We consider a class of fully nonlinear Fermi–Pasta–Ulam (FPU) lattices, consisting of a chain of particles coupled by fractional power nonlinearities of order *α*>1. This class of systems incorporates a classical Hertzian model describing acoustic wave propagation in chains of touching beads in the absence of precompression. We analyse the propagation of localized waves when *α* is close to unity. Solutions varying slowly in space and time are searched with an appropriate scaling, and two asymptotic models of the chain of particles are derived consistently. The first one is a logarithmic Korteweg–de Vries (KdV) equation and possesses linearly orbitally stable Gaussian solitary wave solutions. The second model consists of a generalized KdV equation with Hölder-continuous fractional power nonlinearity and admits compacton solutions, i.e. solitary waves with compact support. When , we numerically establish the asymptotically Gaussian shape of exact FPU solitary waves with near-sonic speed and analytically check the pointwise convergence of compactons towards the limiting Gaussian profile.

## 1. Introduction

The problem of analysing the response of a nonlinear lattice to a localized disturbance arises in many applications, such as the study of stress waves in granular media after an impact [1,2], the excitation of nonlinear oscillations in crystals by atom bombardment [3,4] or the response of nonlinear transmission lines to a voltage pulse [5]. Several important dynamical phenomena can be captured by the Fermi–Pasta–Ulam (FPU) model [6] consisting of a chain of particles coupled by a pairwise interaction potential *V* . The dynamical equations for a spatially homogeneous FPU chain read as follows:
1.1where is the displacement of the *n*th particle from a reference position. System (1.1) can be rewritten in terms of the relative displacements *u*_{n}=*x*_{n}−*x*_{n−1} and particle velocities as follows:
1.2The dynamical evolution of localized solutions of (1.2) is strongly influenced by the properties of the interaction potential *V* . In its most general form, the interaction potential satisfies
1.3where *α*>1 and *κ*≥0.

In the work of Mielke & Patz [7], the dispersive stability of the zero equilibrium state is proved for *κ*>0 and *α*>4, i.e. for sufficiently weak nonlinearities near the origin. More precisely, the amplitude (i.e. supremum norm) of the solution of the FPU lattice (1.2) goes to 0 when for all initial conditions sufficiently small in ℓ^{1}, where ℓ^{1} denotes the classical Banach space of bi-infinite summable sequences.

By contrast, in many situations nonlinear effects are strong enough to compensate dispersion, yielding the existence of coherent localized solutions of the FPU lattice (1.2), such as solitary waves propagating at constant speed, or time-periodic breathers (see e.g. [6] for a review). The first existence result for solitary waves in a general class of FPU lattices was obtained by Friesecke & Wattis [8], when *V* has the local minimum (not necessarily strict) at the origin and is superquadratic at one side (see also [9] and references therein). In addition, the existence of solitary waves near the so-called long wave limit was established in [10–14] for smooth (*C*^{3}) potentials *V* . More precisely, for *κ*>0 and *V* ′′′(0)≠0 (i.e. *α*=2 in (1.3)), there exist a family of small amplitude solitary waves parametrized by their velocity , where *c*_{s} defines the ‘sound velocity’ of linear waves. These solutions take the form
where and *z*(*η*)=sech^{2}(*η*/2). In particular, these solitary waves decay exponentially in space and broaden in the limit of vanishing amplitude. Equivalently, one has
1.4where *ξ*:=*ϵ*(*n*−*c*_{s}*t*), *τ*:=*ϵ*^{3}*c*_{s}*t*/24 and *y*(*ξ*,*τ*):=*z*(*ξ*−*τ*) is a solitary wave solution of the Korteweg–de Vries (KdV) equation
1.5More generally, the solutions of KdV equation (1.5) yield solutions of the FPU system of form (1.4), valid on a time scale of order *ϵ*^{−3} [15–17]. In addition, the nonlinear stability of small amplitude FPU solitary waves was proved in [10–13,18,19], as well as the existence and stability of asymptotic *N*-soliton solutions [20,21]. These results allow to describe in particular the propagation of compression solitary waves in homogeneous granular chains under precompression [1].

Another interesting case corresponds to fully nonlinear interaction potentials, where *κ*=0 (which corresponds to a vanishing sound velocity, that is, *c*_{s}=0) and *V* has the local minimum at the origin. A classic example is given by the Hertzian potential
1.6with *α*>1, where we denote by *H* the Heaviside step function. This potential describes the contact force between two initially tangent elastic bodies (in the absence of precompression) after a small relative displacement *x* [22]. The most classical case is obtained for *α*=3/2 and corresponds to contact between spheres, or more generally two smooth non-conforming surfaces. More recently, granular chains involving different orders of nonlinearity have attracted much attention, see [23,24] and references therein. In particular, experimental and numerical studies on solitary wave propagation have been performed with chains of hollow spherical particles of different width [25] and chains of cylinders [26], leading to different values *α* in the range of 1.15≤*α*≤1.5 (see also [27] for other systems with *α* close to unity).

The propagation of stationary compression pulses in the FPU lattice (1.1) with potential (1.6) for *α*=3/2 was first analysed by Nesterenko [1]. These results rely on a formal continuum limit and provide approximate solitary wave solutions with compact support. An alternate continuum limit problem has been introduced in [28] for arbitrary values of *α*>1, leading to different (compactly supported) approximations of solitary waves. The existence of exact solitary wave solutions of FPU lattice (1.1) with potential (1.6) follows from the general result of Friesecke & Wattis [8] mentioned previously (see also [29,30]). The width of these solitary waves is independent of their amplitude owing to the homogeneous nonlinearity of the Hertzian potential. In addition, the fully nonlinear character of the Hertzian potential induces a doubly exponential spatial decay of solitary waves [30,31].

While the above analytical results provide useful informations on strongly localized solitary waves, they are not entirely satisfactory for several reasons. First of all, the existence result of Friesecke & Wattis [8] does not provide an approximation of the solitary wave profile, and the approximations available in the literature [1,28] rely on a ‘long wave’ assumption that is not justified (for example, the solitary waves considered in [1] are approximately localized on five particles). In addition, the dynamical properties of solitary waves in fully nonlinear FPU lattices are not yet understood. Indeed, no mathematical results are available concerning their stability, the way they are affected by lattice inhomogeneities or the existence of *N*-soliton solutions. Another interesting problem is to characterize the excitation of one or several solitary waves from a localized initial perturbation [32,33]. For *c*_{s}≠0 and small amplitude long waves, this problem can be partially analysed in the framework of KdV approximation by using the inverse scattering transform methods [34], but such reduction is presently unavailable for fully nonlinear FPU lattices. These questions are important for the analysis of impact propagation in granular media, and more generally for the design of multiple impact laws in multi-body mechanical systems [2,32].

In this paper, we attack the problem by considering a suitable long wave limit of fully nonlinear FPU lattices. We consider FPU lattice (1.1) with the homogeneous fully nonlinear interaction potential
1.7with *α*>1. Obviously, all solutions *u*_{n}≤0 of FPU lattice (1.1) with the potential (1.7) are also solutions of Hertzian FPU lattice (1.1) and (1.6). The problem can be rewritten in terms of the relative displacements in the following way
1.8where we denote and (Δ*u*)_{n}=*u*_{n+1}−2*u*_{n}+*u*_{n−1} is the discrete Laplacian. For approximating the temporal dynamics of (1.8) in a continuum limit, fully nonlinear versions of the Boussinesq equation considered in [1,28] possess serious drawbacks, as they may lead to blow-up phenomena in analogy with the classical ‘bad’ Boussinesq equation [35]. In §2, we numerically show that these models introduce artificial dynamical instabilities with arbitrarily large growth rates, which suggests ill-posedness of these equations [36]. In §3, instead of using a Boussinesq-type model, we then formally derive a logarithmic KdV (log-KdV) equation as a modulation equation for long waves in fully nonlinear FPU lattices, obtained in the limit (§3*a*). The log-KdV equation takes the form
1.9and provides approximate solutions *u*_{n}(*t*)≈*v*(*ξ*,*τ*) of the original FPU lattice (1.8) for , and .

Log-KdV equation (1.9) admits Gaussian solitary wave solutions (§3*b*), which have been previously identified as solutions of the stationary logarithmic nonlinear Schrödinger equation (log-NLS) in the context of nonlinear wave mechanics [37]. Closer to our case, Gaussian homoclinic solutions have been also found to approximate the envelope of stationary breather solutions in Newton's cradle (i.e. system (1.1) and (1.6) with an additional on-site potential) in the limit [38]. In §3*b*, we numerically check that solitary wave solutions of the Hertzian FPU lattice with velocity *v*_{s}=1+*c*(*α*−1) converge towards Gaussian approximations when is fixed and . These solitary waves have velocities close to unity, which corresponds to the value of sound velocity in the linear chain with *α*=1. In addition, we check that the FPU solitary waves are well approximated by the compacton solutions derived in [28] when *α*∈(1,3/2]. To go beyond the stationary regime, we check numerically that the Gaussian approximation captures the asymptotic shape of a stable pulse forming after a localized velocity perturbation in the Hertzian FPU lattice (1.1)–(1.6) with *α*≈1 (§3*c*). Consistently with the above dynamical simulations, we prove in §3*d* the linear orbital stability of Gaussian solitary waves for the log-KdV equation. Our analysis makes use of a suitable convex conserved Lyapunov function, but negative index techniques developed in recent works [39,40] for KdV-type equations would also apply.

The link between Gaussian solitary waves and compactons is made explicit in §4, where we check the pointwise convergence of the compacton solutions of [28] towards Gaussian profiles when . In addition, following the methodology developed in §3*a*, we derive from the fully nonlinear FPU lattice a generalized KdV equation with Hölder-continuous nonlinearity (H-KdV):
1.10When , the H-KdV equation (1.10) is consistent with the FPU lattice in the sense that each solution to this equation ‘almost’ satisfies (1.8) up to a small residual error. Equation (1.10) admits explicit compacton solutions whose form is close to the compactons obtained in [28] with the use of a Boussinesq-type model. When , these solutions converge towards the Gaussian solitary waves studied in §3*b*, and thus they provide an (asymptotically exact) approximation of FPU solitary waves with near-sonic speed. This result sheds a new light on the compacton approximations for FPU solitary waves heuristically derived in the literature [1,28]. Another interest of the H-KdV equation lies in the (non-differentiable) Hölder-continuous nonlinearity *v*|*v*|^{1/α−1}, which allows for the existence of compactons. This type of degeneracy is quite different from the classical feature of compacton equations which incorporate degenerate nonlinear dispersion [36,41].

We finish this paper with a summary of our results and a discussion of several open questions concerning the qualitative dynamics of the log-KdV and H-KdV equations and their connections with fully nonlinear FPU chains (§5).

## 2. Fully nonlinear Boussinesq equation and compactons

Fully nonlinear Boussinesq equations have been introduced in [1,28] as formal continuum limits of FPU chains with Hertzian-type potentials. In the study of Nesterenko [1], the continuum limit is performed on system (1.1) with potential (1.7) describing particle displacements, whereas [28] considers system (1.8) for relative displacements. In what follows, we discuss the continuum limit introduced in [28], which takes a slightly simpler form than the system derived in [1].

The fully nonlinear Boussinesq equation introduced in [28] takes the form
2.1where *u*_{|x=n} denotes an approximation of a solution *u*_{n} of (1.8). The right-hand side of (2.1) is obtained by keeping the first two terms of the formal Taylor expansion of the discrete Laplacian in (1.8)
2.2This truncation is purely formal, but a numerical justification is presented in [28] in the particular case of solitary wave solutions. More precisely, the solitary waves *u*_{n}(*t*)=*w*(*z*), *z*=*n*−*t* of (1.8) are numerically compared with solitary wave solutions *u*(*x*,*t*)=*w*(*z*), *z*=*x*−*t* of (2.1). For this class of solutions, equation (2.1) reduces to a fourth-order ordinary differential equation, which can be integrated twice and leads to
2.3whereas equation (1.8) reduces to the differential advance-delay equation
2.4with (Δ*w*)(*z*)=*w*(*z*+1)−2*w*(*z*)+*w*(*z*−1). The wave velocity can be normalized to unity owing to a scaling invariance of the FPU system (1.1) with homogeneous potential (1.7) (or Hertzian potential (1.6)), namely each solution *u*_{n} generates a one-parameter family of solutions |*v*_{s}|^{2/(α−1)}*u*_{n}(*v*_{s}*t*) with (the same scaling invariance exists in system (2.1)).

According to the numerical computations presented in [28], the solitary wave of the differential advance-delay equation (2.4) is well approximated by the compactly supported solitary wave of differential equation (2.3) for *α*=3/2, and the discrepancy increases with *α*. The compacton solution of (2.3) found in [28] takes the form
2.5where

In what follows, we re-examine the consistency of (1.8) and (2.1), from a dynamical point of view, by analysing the spectral stability of compactons. Linearizing (2.1) at the compacton *w*_{c} in the reference frame travelling with unit velocity, we use the ansatz *u*(*x*,*t*)=*w*_{c}(*x*−*t*)+*U*(*x*−*t*) e^{λt}, where *λ* is the spectral parameter and *U* is the perturbation term. We arrive at the spectral problem
2.6where and **1** denotes the characteristic function. One can note that *k*_{α} and vanish at the endpoints *z*=±*π*/2*B* of the compact support of *w*_{c}. We look for eigenvectors in the Hilbert space
where *H*^{m} denotes a classical Sobolev space [42]. As and *k*_{α}′′′ are discontinuous at *z*=±*π*/2*B*, this yields the condition
2.7This allows us to reduce eigenvalue problem (2.6) to the compact interval [−*π*/2*B*,*π*/2*B*] with boundary conditions (2.7) and approximate the spectrum with the standard finite difference method (we have used second-order difference approximations for derivatives and 2000 grid points). If there exist eigenvalues with Re(*λ*)>0, then the solitary wave is spectrally unstable. If all the eigenvalues are located at the imaginary axis, then the solitary wave is called spectrally stable.

Figure 1 shows the complex eigenvalues *λ* of the spectral problem (2.6) and (2.7) for *α*=1.05,1.2,1.5. Spectrum is invariant under and (note the presence of a small number non-symmetric eigenvalues, which originate from numerical errors). We find the existence of unstable eigenvalues for all values of *α*>1 considered, and the eigenvalues approach the real line far from the origin (this part of the spectrum is not visible in figure 1*a*,*b*). Consequently, these results imply the spectral instability of compacton (2.5) in system (2.1). Note that the usual notion of instability may not be well defined, because the evolution problem (2.1) may not be well posed. Indeed, our numerical results indicate that the spectrum of (2.6) and (2.7) is unbounded in the positive half-plane (in fact at both sides of the imaginary axis), and thus the linearized evolution problem may be ill-posed. We conjecture that ill-posedness occurs also in system (2.1), in analogy with ill-posedness results recently obtained in [36] for certain nonlinear degenerate dispersive equations.

Along these lines, it is interesting to consider the limit case of the spectral problem (2.6) and (2.7) when . As for all and , the limiting spectral problem possesses constant coefficients and is defined on the entire real line. Using the Fourier transform, one can compute the (purely continuous) spectrum explicitly, which yields 2.8This limit case is represented in figure 1 (grey curves). The spectrum being unbounded in the positive half-plane, the corresponding linear evolution problem is then ill-posed.

The above instability phenomena are not physically meaningful because the solitary waves are known to be stable from simulations of impacts in Hertzian chains [1]. In equation (2.8) obtained in the limit , these instabilities occur for short wavelengths (with *k*^{2}>12), whose dynamical evolution cannot be correctly captured by continuum limit (2.1). In the next section, we derive a different asymptotic model free of such artificial instabilities.

### Remark 2.1

A classical model introduced to correct artificial short wavelength instabilities of (2.1) corresponds to the regularized Boussinesq equation
2.9(e.g. [43]). This model has the inconvenience of altering the spatial decay of solitary waves in the case of fully nonlinear interaction potentials. Indeed, looking for travelling wave solutions *u*(*x*,*t*)=*w*(*z*), *z*=*x*−*t* and integrating equation (2.9) twice, one obtains
2.10after setting two integration constants to 0. Equation (2.10) admits non-trivial symmetric homoclinic solutions ±*w*_{α}(*z*) satisfying , corresponding to solitary wave solutions of (2.9). These solutions decay exponentially in space, which is too slow compared with the superexponential decay of the solitary wave solutions of (2.4).

## 3. The log-KdV equation and Gaussian solitary waves

### (a) Formal derivation of the log-KdV equation

In order to pass to the limit for long waves, it is convenient to rewrite (1.8) in the form
3.1where
3.2(uniformly in *u* on bounded intervals) when . For *α*=1, we have *f*_{1}(*u*)=0 for all and system (3.1) reduces to a semi-discrete linear wave equation. In that case, the scaling in (1.4) (with *c*_{s}=1) yields a linearized KdV equation for the envelope function *y*. To analyse the limit , we assume the same type of scaling for the solution *u*, i.e. we search for solutions depending on slow variables and , where *ϵ*>0 is a small parameter. We look for solutions of the form
3.3In contrast with (1.4), the leading term *v*(*ξ*,*τ*) is assumed of order unity and the remainder term in (3.3) is . Of course, owing to the scaling invariance of the FPU lattice (1.8) for *α*>1, solutions with arbitrarily small or large amplitudes can be deduced from any solution of form (3.3).

From scaling (3.3) and using a Taylor expansion and the chain rule, we obtain
3.4and
3.5To evaluate the right-hand side of (3.1), we use the expansion
where the logarithmic remainder term accounts for the possible vanishing of *v*. Setting now and using (3.4), we obtain
3.6With this choice of *ϵ*, the left- and right-hand sides of equation (3.1) have the same order *ϵ*^{4} according to expansions (3.5) and (3.6). Substituting these expansions in (3.1) yields
3.7Then neglecting the higher order terms and integrating with respect to *ξ* leads to
3.8where the integration constant has been fixed to 0 in order to cover the case when . We shall call equation (3.8) the logarithmic KdV (log-KdV) equation. It can be rewritten
3.9where the potential *W* reads
Equation (3.7) shows that the log-KdV equation is consistent with the nonlinear lattice (1.8), i.e. each solution of (3.8) is almost a solution of (1.8) up to a small residual error.

Note that if *v* is a solution of (3.8), so is −*v*. In addition, equation (3.8) admits a non-standard Galilean symmetry involving a rescaling of amplitude, i.e. each solution *v* generates a one-parameter family of solutions
3.10In particular, all travelling wave solutions of the log-KdV equation (3.8) can be deduced from its stationary solutions. This property is inherited from the scaling invariance of FPU system (1.8) (this point will be detailed in §3*b*).

Equation (3.8) falls within the class of generalized KdV equations. Systems in this class possess three (formally) conserved quantities [44], namely the mass
3.11the momentum
3.12and the energy
3.13Well-posedness results for the Cauchy problem associated with (3.9) are known when *W*^{′′} is a *C*^{2} function [44], but the existing theory does not apply to our case, where diverges logarithmically at the origin.

### (b) Stationary solutions

Looking for solutions of log-KdV equation (3.8) depending only on *ξ*, one obtains the stationary log-KdV equation
3.14Integrating once under the assumption , one obtains
3.15This equation can be seen as a (one-dimensional) stationary log-NLS equation [37].

The potential *W* in (3.15) has a double-well structure with the local maximum at *v*=0 (figure 2), hence there exists a pair of (symmetric) homoclinic orbits to 0 and a continuum of periodic orbits. The homoclinic solutions have the explicit form *v*(*ξ*)=±*v*^{0}(*ξ*) with
3.16Note that these (Gaussian) homoclinic solutions decay super-exponentially, but do not decay doubly exponentially unlike the solitary wave solutions of the differential advance-delay equation (2.4) [30,31].

Homoclinic solution (3.16) yields an approximate Gaussian solitary wave solution of FPU lattice (1.8) with velocity equal to unity
3.17with
3.18Figure 3 compares the solitary wave solution of the differential advance-delay equation (2.4) computed numerically with the analytical approximations corresponding to compactly supported solitary wave (2.5) and Gaussian solitary wave (3.18). The numerical approximations of solitary wave solutions of (2.4) were obtained using the algorithm described in Ahnert & Pikovsky [28], based on a reformulation of (2.4) as a nonlinear integral equation and the method of successive approximations (see also [9,31] for variants of this method). Figure 4 shows the relative error (in norms) between solitary wave solutions of the differential advance-delay equation (2.4) and the two approximations (2.5) and (3.18) as a function of *α*. The Gaussian solitary wave provides a worse approximation compared with the compactly supported solitary wave, but both approximation errors converge to zero, when . We have in addition ; hence the absolute errors between the exact solitary wave *w* and the two approximations converge to zero similar to the relative errors plotted in figure 4.

So far, we have computed a solitary wave solution of the FPU lattice (1.8) with unit velocity and have checked its convergence towards the Gaussian approximation (3.17) when . We shall now examine the convergence of solitary waves with velocities different from unity. Using the Galilean invariance of (3.8), homoclinic solution (3.16) yields two (symmetric) families of solitary wave solutions of log-KdV equation (3.8)
3.19parametrized by the wave velocity . These profiles yield the approximate solitary wave solutions of the original FPU lattice (1.8)
3.20where we have set ,
3.21introduced an additional phase shift and used the fact that . One can observe that the width of the approximate solitary wave (3.20) diverges as (*α*−1)^{−1/2} when . Moreover, similar to solitary wave solutions of the FPU lattice (1.8), the wave width remains constant if *α* is fixed and the wave amplitude (or equivalently the wave velocity *v*_{s}) is varied. In addition, approximation (3.20) can be rewritten as
3.22where the renormalized Gaussian profile *w*_{G} is defined in (3.18). One can note that if is fixed . Consequently, approximation (3.22) is close to a rescaling of (3.17) through the invariance *u*_{n}↦|*v*_{s}|^{2/(α−1)}*u*_{n}(*v*_{s}*t*) of (1.8). This observation illustrates why the (non-standard) Galilean invariance (3.10) of the log-KdV equation is inherited from the scaling invariance of the FPU lattice (1.8).

Similarly, the solitary wave solution *w* of the differential advance-delay equation (2.4) yields the family of solitary wave solutions of (1.8)
3.23One can distinguish two contributions to the error between the Gaussian approximation (3.22) and exact solution (3.23), one originating from the profile function *w*_{G} and the other from the wave amplitude. From the numerical results of figure 4, we know that when . In addition, fixing and considering wave velocities (3.21) close to unity when , we have seen that . Consequently, exact solitary wave (3.23) with velocity (3.21) converges uniformly towards Gaussian approximation (3.22) when .

### Remark 3.1

Note that the convergence result above concerns solitary waves with velocities converging towards unity, i.e. the value of sound velocity in the linear chain with *α*=1. This restriction is owing to specific scaling (3.3) assumed for solutions described by the log-KdV equation. On the contrary, exact FPU solitary wave (3.23) with fixed velocity *v*_{s}≠±1 becomes degenerate when , as the wave amplitude goes to 0 if |*v*_{s}|<1 and diverges if |*v*_{s}|>1.

### Remark 3.2

Compactly supported solitary wave (2.5) yields approximate solutions *u*^{c} of FPU system (1.8) taking the form . This approximation can be compared to the exact FPU solitary wave *u* defined by (3.23). The results of figure 4 show that as *α*→1^{+}. Consequently, one can infer that the relative error converges to zero when *v*_{s} is fixed and .

### (c) Formation of Gaussian solitary waves

In §3*b*, we have computed solitary wave solutions of the FPU lattice (1.8) and have checked their convergence towards Gaussian approximation (3.20) when . These results are valid in a stationary regime and for prescribed wave velocities converging towards unity. To complete this analysis, we shall study the formation of solitary waves from a localized perturbation of given magnitude and compare their profiles to Gaussian approximations when *α* is close to one.

In what follows, we numerically integrate FPU system (1.1) with Hertzian potential (1.6) for different values of *α*>1. We consider Hertzian potential (1.6) rather than symmetrized potential (1.7) because of its relevance to impact mechanics. In addition, differential equations (1.1) are easier to integrate numerically owing to the absence of dispersive wavetrains for the above initial condition.

We consider a lattice of *N*=2000 particle with free-end boundary conditions. Computations are performed with the standard ODE solver of the software Scilab. We consider a velocity perturbation of the first particle (at *n*=0), corresponding to the initial condition
3.24Owing to the scale invariance of system (1.1) and (1.6), all positive initial velocities for yield a rescaled solution of the form .

The front edge of the solution evolves into a solitary wave whose profile becomes stationary for large enough times, at least on the timescales of the simulations. When *α* is sufficiently close to unity, one observes that the asymptotic velocity of the solitary wave is close to unity (this property is also true for any initial velocity ). The solitary wave is of compression type (i.e. with *u*_{n}<0), hence it also defines a solution of the FPU lattice (1.8) with symmetrized potential (1.7) and can be compared with approximation (3.20).

The results are shown in figure 5 for *α*=1.01. The parameter *c*≈−2.07 in the Gaussian approximation (3.20) is determined from the relation , where *a* is the exact solitary wave amplitude obtained by integrating (1.1) with initial data (3.24). The approximation of the solitary wave profile is very accurate in the stationary regime, as shown by figure 5*c*. In addition, the measured velocity of the numerical solitary wave solution *v*_{num}≈0.9789 can be compared with the velocity *v*_{app}=1+*c*(*α*−1)≈0.9793 of the approximate solitary wave (3.20), which yields a relative error *E*=|*v*_{app}−*v*_{num}|/*v*_{num} around 0.04%.

Discrepancies appear between the profiles of the numerical solution and the Gaussian approximation for larger values of *α*, as already noted in figures 3 and 4. In addition, we find a relative error *E* between numerical and approximate wave velocities around 7% for *α*=1.222 and 36% for *α*=1.5. As a conclusion, while some quantitative agreement is still obtained for *α*≈1.2 between the numerical solution and the Gaussian approximation, the latter becomes unsatisfactory for *α*=1.5.

### (d) Linear stability of Gaussian solitary waves

The numerical results of §3*c* indicate the long time stability of the solitary wave solutions of system (1.1) with Hertzian potential (1.6) that form after a localized perturbation. Therefore, the stability of Gaussian solitary waves (3.19) appears as a necessary (of course not sufficient) condition to establish the validity of log-KdV equation (3.8) as a modulation equation for the FPU system with Hertzian potential. In this section, we prove the linear orbital stability of solitary waves of the log-KdV equation. We perform the analysis for the stationary Gaussian solution *v*^{0} defined by (3.16). By scaling transformation (3.10), the stability result extends to the entire family of solitary waves e^{c}*v*^{0}(*ξ*−*cτ*) with .

Log-KdV equation (3.9) can be written in the Hamiltonian form
3.25associated with the energy (3.13). The stationary Gaussian solution *v*^{0} is a critical point of the energy *E*(*v*), i.e. *E*′(*v*^{0})=0. The Hessian operator evaluated at this solution reads
3.26and corresponds to a Schrödinger operator with a harmonic potential. Equation (3.25) linearized at *v*^{0} reads
3.27In what follows, we formulate (3.27) as a differential equation in suitable function spaces and derive a linear stability result based on the energy method for KdV-type evolution equations [45]. Hereafter, we assume that the reader is familiar with this method and basic spectral theory.

The spectral properties of *L* are well known [46]. The operator *L* is self-adjoint in with dense domain
Its spectrum consists of simple eigenvalues at integers *n*−1, where (the set of natural numbers including zero). The corresponding eigenspaces are spanned by rescaled Hermite functions . In particular, *L* has a simple eigenvalue −1 with eigenspace spanned by *ϕ*_{0}=*v*^{0}/∥*v*^{0}∥_{2}, a simple zero eigenvalue with eigenspace spanned by , and the rest of its spectrum is bounded away from zero by a positive number. The discreteness of the spectrum comes from the fact that the harmonic potential of the Schrödinger operator *L* is unbounded at infinity, which implies that the embedding is compact (e.g. [46, p. 43–44]) and *L* has a compact resolvent in .

Hereafter, we denote by (⋅,⋅) the usual *L*^{2}-scalar product. The operator inherits a double non-semi-simple zero eigenvalue, with generalized kernel *E*_{0} spanned by the eigenvector and the generalized eigenvector *ϕ*_{0}. The algebraic multiplicity of this eigenvalue is 2 because the equation ∂_{ξ}*Lu*=*ϕ*_{0} has no solution in *D*(*L*) (because ). The double zero eigenvalue is linked with the existence of a two-parameter family of solitary waves of log-KdV equation (3.8) parametrized by the location and velocity of the waves. It induces in general a secular growth of the solutions of (3.27) along the eigenvector , linked with a velocity change of perturbed solitary waves. In order to prove a linear stability result, we thus have to project (3.27) onto the invariant subspace under ∂_{ξ}*L* associated with the non-zero part of the spectrum. Following a classical computation scheme [47], eqns (2.10), (2.11), (2.19), the spectral projection onto *E*_{0} takes the form with
(one can readily check that *P*_{0} commutes with ∂_{ξ}*L* using integration by parts). Note in passing that *a* is well defined because . Now, splitting the solutions *v*(⋅,*τ*)∈*D*(*L*) of (3.27) into
with *y*(⋅,*τ*)=(*I*−*P*_{0})*v*(⋅,*τ*), one obtains the following equivalent system:
3.28and
3.29In order to prove the linear (orbital) stability of the Gaussian solitary wave, we have to show the Lyapunov stability of the equilibrium *y*=0 of linear evolution equation (3.29) for a suitable topology. Let us recall that *P*_{0}*y*(⋅,*τ*)=0, hence *y*(⋅,*τ*) belongs to the codimension-2 subspace of *D*(*L*) defined as
As *L* is positive on *ϕ*_{0}^{⊥} and *D*_{1}⊂*ϕ*_{0}^{⊥}, we can define ∥*y*∥_{L}=(*Ly*,*y*)^{1/2}. Owing to the fact that *L* is positive-definite on and , ∥⋅∥_{L} defines a norm on *D*_{1} (roughly a weighted *H*^{1}-norm). Denote by *H*_{1} the completion of *D*_{1} with respect to the norm ∥⋅∥_{L}. This norm defines a convex conserved Lyapunov function for system (3.29) because
3.30For simplicity, let us choose an initial data in the Schwartz space of rapidly decreasing functions
3.31Owing to property (3.30), we get from the energy method [42] (using the Hilbert basis of Hermite functions for the Galerkin approximation, see e.g. [42], ch. 11.1.2) and standard bootstrapping arguments (e.g. [42, ch. 11.1.4]) a unique global solution of (3.29)–(3.31) which is infinitely smooth in time and space. We have and ∥*y*(⋅,*τ*)∥_{L}=∥*y*_{0}∥_{L} for all , which shows the Lyapunov stability of the equilibrium *y*=0 of (3.29) in *H*_{1}. Therefore, we have proved the linear orbital stability of the Gaussian solitary wave.

The Lyapunov stability of the equilibrium *y*=0 implies the absence of eigenvalue of ∂_{ξ}*L* with positive real part, since such an eigenvalue would lead to exponential growth of the solution along a corresponding eigenvector. In addition, the spectrum of ∂_{ξ}*L* is invariant by as ∂_{ξ}*L* possesses a reversibility symmetry, i.e. anticommutes with the symmetry *v*(*ξ*)↦*v*(−*ξ*). This implies that all the eigenvalues of ∂_{ξ}*L* lie on the imaginary axis. This result contrasts with the instability of the compacton solutions of the fully nonlinear Boussinesq equation numerically analysed in §2. Moreover, it is consistent with the absence of solitary wave instabilities observed in §3*c*.

### Remark 3.3

The absence of eigenvalues of ∂_{ξ}*L* with positive real part could be also obtained from the recent works [39,40]. This result follows from the main theorems in [39,40] if the number of negative eigenvalues of *L* is equal to one and
3.32where . Assumption (3.32) is satisfied as (recall *Lϕ*_{0}=−*ϕ*_{0}) and

## 4. Compacton approximation revisited

The results of figures 3 and 4 indicate that compacton approximations converge towards solitary wave solutions of the differential advance-delay equation (2.4) when . It seems delicate to establish this result directly from the methodology described in §2, where equation (2.3) is heuristically obtained by the truncation of the series expansion (2.2) of the discrete Lapacian in equation (2.4). However, it is instructive to compare analytically the compacton approximation *w*_{c}(*z*) defined in (2.5) and Gaussian approximation (3.18) when *α* is close to unity, because our numerical results indicate that both profiles become very close (figure 3*b*). To check the consistency of (2.5) and (3.18), we note that the compact support [−*π*/2*B*,*π*/2*B*] of approximation (2.5) extends to the entire real line as . Furthermore, let , , and perform the expansions
4.1and
4.2for all fixed . From expansions (4.1) and (4.2), it follows that the renormalized compacton converges towards Gaussian solution (3.16) for any fixed when .

One possible approach for the justification of compactons consists in deriving an asymptotic model consistent with FPU lattice (1.8) and supporting compacton solutions, in analogy with the derivation of the log-KdV equation. Such a model may be free of the artificial instabilities introduced by the Boussinesq equation (2.1), and may lead to a well-posed evolution problem in a suitable function space and for appropriate initial data.

In what follows, we show that a generalized KdV equation with Hölder-continuous fractional power nonlinearity can be derived from (1.8) using the method of §3. As shown subsequently, the limited smoothness of the nonlinear term in this equation allows for the existence of compactons. When , one can rewrite the expansion of the nonlinearity *f*_{α}(*u*) defined by (3.2) in the following form:
4.3uniformly in *u* on bounded intervals. The leading order term at the right-hand side of (4.3) involves a Hölder-continuous nonlinearity *u*|*u*|^{1/α−1}, whereas the function *f*_{α} is more regular (*C*^{1}) in *u*. As it follows from expansion (4.3), the generalized KdV equation
4.4is consistent with system (1.8) at the same order as the log-KdV equation and converges towards the log-KdV equation when . We call equation (4.4) the Hölderian KdV (H-KdV) equation. The above approximation of the logarithmic nonlinearity is reminiscent of results of [38] obtained for a stationary Hölderian NLS equation close to a logarithmic limit.

Note that H-KdV equation (4.4) admits the same three conserved quantities (3.11)–(3.13) with potential *W* replaced by
Moreover, equation (4.4) admits a non-standard Galilean invariance similarly to the log-KdV equation. More precisely, any solution *v* of (4.4) generates a one-parameter family of solutions *ϕ*(*c*)⋅*v* defined by
and parametrized by . One can check that this symmetry reduces to the Galilean invariance (3.10) when .

Let us check that the H-KdV equation admits compacton solutions. Their existence is owing to the non-differentiable Hölder-continuous nonlinearity, in contrast to classical compacton equations where degenerate nonlinear dispersion plays a central role [41]. The stationary H-KdV equation integrated once reads
4.5(the integration constant has been set to 0), or equivalently . This equation is integrable and the potential has a double-well structure with the local maximum at *v*=0. This property implies the existence of a pair of symmetric homoclinic orbits to 0 corresponding to compactons. Using the change of variable
equation (4.5) is mapped to form (2.3) which possesses compacton solutions given by (2.5). Consequently, equation (4.5) admits stationary compacton solutions
4.6where
In addition, we have as above for any fixed , i.e. the compacton approximation (4.6) converges towards the Gaussian solitary wave approximation (3.16) when .

### Remark 4.1

The compacton approximation defined in (2.5) and obtained using Boussinesq equation (2.1) differs from compacton (4.6) deduced from the H-KdV equation. However, both approximations are equivalent when as they converge towards Gaussian profile (3.16). In fact, an infinity of compacton approximations could be constructed, depending on the approximation of the logarithmic nonlinearity introduced in (4.3).

Using the Galilean invariance of (4.4), the compacton (4.6) yields two (symmetric) families of compactly supported solitary waves of H-KdV equation (4.4)
4.7parametrized by the wave velocity . From expression (3.3), these profiles yield the approximate compacton solutions of FPU lattice (1.8)
4.8where we have set , *v*_{s}=1+*c*(*α*−1), , introduced an additional phase shift and used the fact that . For any fixed value of *c* and , we have . In this limit, the compacton (4.8) converges towards Gaussian approximation (3.20), therefore it is also close to exact FPU solitary wave (3.23).

## 5. Discussion

We have obtained two generalized KdV equations (log-KdV equation (3.8) and H-KdV equation (4.4)) as formal asymptotic limits of FPU lattices with Hertzian-type potentials, when the nonlinearity exponent *α*>1 goes to unity and slowly varying profiles are considered. Using numerical computations, we have checked that FPU solitary waves converge towards Gaussian solitary waves and compacton solutions to these KdV equations when and for near-sonic wave speeds. In addition, we have illustrated numerically the formation of stable solitary waves after a localized velocity perturbation in Hertzian FPU system (1.1) and (1.6) when *α*≈1, a limit in which the propagating pulse becomes nearly Gaussian. The linearized log-KdV equation preserves the spectral stability of solitary waves, which is lost when using other formal (Boussinesq-type) continuum models. While our study does not yield a complete proof of the asymptotic behaviour of exact FPU solitary waves when , it provides nevertheless an asymptotic framework to explain classical formal compacton approximations [1,28], whose justification remained unclear up to now.

It would be interesting to examine the dynamical properties of the log-KdV and H-KdV equations for different classes of initial conditions. Relevant questions include local well-posedness (or ill-posedness), derivation of *a priori* bounds, global well-posedness (or blow-up), scattering of some initial data and nonlinear stability of solitary waves. In our context, the study of the nonlinear orbital stability [45] or asymptotic stability [47] of solitary waves rises new difficulties, linked with the lack of smoothness of the energy functional (3.13).

In addition, if the well-posedness of the Cauchy problem for the log-KdV or H-KdV equation could be established for appropriate initial data, one could then study analytically or numerically their connection with FPU system (1.1) with homogeneous potentials (1.7). An open question is to check that well-prepared initial data evolve (up to higher order terms and on long finite times) according to the log-KdV or H-KdV equation (in the same spirit as the justification of the classical KdV equation for FPU chains [15–17]). This problem may be extended to Hertzian FPU system (1.1) and (1.6), at least close to a solitary wave solution and for a suitable topology (i.e. using a weighted norm flattening perturbations behind the propagating wave [10–13,47]). The construction of appropriate numerical methods for the time-integration of the log-KdV or H-KdV equations is of course also fundamental in this context. Another open problem concerns the dynamical stability of the solitary wave solutions of the FPU lattice (1.8). Our proof of the linear orbital stability of solitary waves for the log-KdV equation could be useful in this context, following the lines of [10–13] (using the linearized log-KdV equation instead of the linearized KdV equation). These open questions will be explored in forthcoming works.

## Funding statement

G.J. acknowledges financial support from the French Embassy in Canada and the Rhône-Alpes Complex Systems Institute (IXXI). D.P. acknowledges financial support from the NSERC.

## Acknowledgements

G.J. is grateful to V. Acary, B. Brogliato, W. Craig and Y. Starosvetsky for stimulating discussions on this topic. Part of this work was carried out during a visit of G.J. to the Department of Mathematics at McMaster University, to which G.J. is grateful for hospitality. G.J. introduced the time-dependent logarithmic and Hölderian KdV equations and their derivations from the Hertzian chain model. Numerical computations were performed by D.P. (figures 1, 3 and 4) and G.J. (figure 5). All authors contributed to the linear stability analysis of Gaussian solitary waves. D.P. showed the convergence of the Ahnert-Pikovsky compacton towards a Gaussian for Hertz force exponents close to unity. All authors contributed to the writing and editing of the manuscript.

- Received July 12, 2013.
- Accepted January 24, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.