Improved models of dense anharmonic lattices

Physics Letters A 381 (2017) 87–93

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Improved models of dense anharmonic lattices P. Rosenau ∗ , A. Zilburg School of Mathematical Sciences, Tel-Aviv University, Tel Aviv 69978, Israel

a r t i c l e

i n f o

Article history: Received 12 July 2016 Received in revised form 26 October 2016 Accepted 27 October 2016 Available online 2 November 2016 Communicated by C.R. Doering Keywords: Anharmonic lattices Compactons Improved PDE description Quartic and Hertz potentials

a b s t r a c t We present two improved quasi-continuous models of dense, strictly anharmonic chains. The direct expansion which includes the leading effect due to lattice dispersion, results in a Boussinesq-type PDE with a compacton as its basic solitary mode. Without increasing its complexity we improve the model by including additional terms in the expanded interparticle potential with the resulting compacton having a milder singularity at its edges. A particular care is applied to the Hertz potential due to its nonanalyticity. Since, however, the PDEs of both the basic and the improved model are ill posed, they are unsuitable for a study of chains dynamics. Using the bond length as a state variable we manipulate its dispersion and derive a well posed fourth order PDE. © 2016 Elsevier B.V. All rights reserved.

1. Introduction A wide variety of physical systems are described by anharmonic lattices that on the macroscopic length-scale involve many unit cells [1] and thus in a first approximation are modeled using the continuum limit given in terms of partial differential equation(s) (PDE), with the discrete effects being completely washed away. Since, however, with few notable exceptions [7], the continuum limit leads to gradients catastrophe, the leading effect due to the discreteness was retained to yield a quasi-continuum PDE description of the lattice. The resulting solitary mode – the compacton [2] – is an example of energy-storing, essentially nonlinear entity and is both of fundamental importance and of longstanding interest [1, 3–5,11]. Note that although, in principle, lattices governed by Newton’s laws immediately spread any information posed over a finite domain, the purely anharmonic interparticle interaction creates a genuine screening effect beyond which there is no measurable excitation. In fact the quasi-continuum limit correctly predicts the main span of the discreton – the solitary solution of the anharmonic chain, which decays at a doubly exponential rate. In spite of the desired features of the compacton, the underlying PDE model hardly solves our quest for an analytically accessible model to study interactions on the lattice. For not only is the expansion not asymptotic and its utility unknown a priori, but the PDE itself is ill posed (occasionally referred to as a bad-Boussinesq

*

Corresponding author. E-mail address: [email protected] (P. Rosenau).

http://dx.doi.org/10.1016/j.physleta.2016.10.051 0375-9601/© 2016 Elsevier B.V. All rights reserved.

Equation). Its main use was to extract an approximate shape of the underlying compacton to be used as an input in the original lattice. Yet the fundamental importance of the dynamics on the lattice and, apart from few exceptions, our inability to analyze it, makes it highly desirable to derive a usable model. Note that the asymptotic procedure devised by Kruskal which resulted in the rebirth of the KdV, is applicable only to weakly anharmonic lattices. With these issues in mind we reexamine a class of purely and almost purely anharmonic lattices and demonstrate how a judicious use of the expanded interparticle potential and the discreteness begets equations with compactons which have an improved level of regularity and provide a better rendition of their discrete antecedents. However, to overcome the ill-posedness of the PDE models based on a direct expansion, we adopt another maneuver which begets a well posed fourth order PDE (capable of handling two sided propagation and head on collisions) which is then optimized for the best L 1 fit of its solitary waves with their discrete antecedents. This approach has two main drawbacks, it is limited to 1-D chains and the tails decay at the conventional exponential rate missing the screening effect of their discrete antecedents. In passing we stress again that our goal is not to duplicate the discrete profile of traveling waves (an exact solution profile can be computed following [13], see also the Appendix) but to derive a better PDE model of lattices dynamics which hopefully will enable a better handle on lattices dynamics. In fact with few notable exceptions the dynamics of nonlinear dense chains is beyond our ability to analyze it and apart from brute force simulation of the lattice, modelization via a PDE seems at this time to be the only tool available to us.

88

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

2. Equations of motion and PDE approximation Consider the Hamiltonian



H disc = h

N  ρ y˙ n2

2

n=− N

+P





y n +1 − y n h



(1)

describing a chain of 2N + 1 particles of equal density ρ = m/h, equal spatial separation h and P being the interaction potential. The equations of motion are

y¨ n =

1 h



P

or for un =

u¨ n =



y n +1 − y n h



− P



y n − y n −1 h

 (2)

y n − y n −1 h

1  P (un+1 ) − 2P  (un ) + P  (un−1 ) . 2 h

(3)

To approximate by a continuum we fix the total length L = 2Nh while letting h ↓ 0 (and hence N ↑ ∞) which yields the continuum limit

 

ρ yt2

H cont =

2

We eliminate tion



+ P ( y x ) dx.

(4)

ρ by rescaling t yielding a nonlinear wave equa-

ytt = ∂∂x P  ( y x ) ,

P  ( S ) ≥ 0.

(5)

In the strict continuum limit not only do the solitary waves observed on a lattice disappear but, as is well known, for any reasonable initial excitation, the second derivatives of solutions may become infinite in a finite time for it is the dispersion due to the discrete lattice which prevents this blow-up, hence its singular effect on the dynamics. The nonlinear spatial gradients are a major change from the linear wave equation; for now we are in a realm typical for quasilinear wave equations. So the issue facing us is to construct an equation which incorporates the discrete effects in an adequate manner. However, since the problem is essentially nonlinear which amounts to saying that it has no phonons, there is no weakly nonlinear regime and an asymptotic expansion associated  with it. With that in mind, we expand the potential P

yn +1− yn h

in h (see Eq. (10)). Keeping the

leading order correction and assuming P  (s) ≥ 0 and smooth we get

ytt = ∂∂x P  ( y x ) +

h2 ∂ 12∂ x



P  ( y x ) ∂∂x





P  ( y x ) y xx + O (h4 ).

(6)

Define u = y x , β = h2 /12, differentiate once and rescale

(x, t ) = ( β z, β τ ) to obtain

2 2 u τ τ = ∂∂z2 P  (u ) + ∂∂z2

P  (u ) ∂∂z





P  (u )u z .

 P

y n +1 − y n h



⇒ P ( yx) +

+ a2 h2 P  ( y xx )2 + a4 h4 P (4) ( y xx )4 + a6 h6 P (6) ( y xx )6 + ... + b2 h4 P  ( y 3x )2 + b4 h8 P (4) ( y 3x )4 + b6 h12 P (6) ( y 3x )6 + ... + c 2 h6 P  ( y 4x )2 + ... + ...

(10)

where a2 = −1/24, a4 = −1/(10 · 242 ), b2 = 1/6!, etc. Clearly this expansion assumes an analytic P . In the following we will present two concrete examples: one where this assumption is satisfied and another where it is not [6] and further measures should be taken. If P is a polynomial of order 2K, P (2K +1) = 0 and each row terminates after K terms. First row terms render a fourth order term in the PDE, the second row begets a sixth order term, and so on. To keep complexity in check we shall not venture beyond a fourth order PDE, which amounts to using only the first row in the expansion of P . First, let us define

.

P 2m = P ( y x )(2m)

(8)

 P

y n +1 − y n h



∼ = P ( yx) +

2m−1

2m P 2m2m ( y xx )2m−2 ∂∂x

=

1) It is ill posed. 2) The expansion leading to Eq. (8) is not asymptotic. The first feature is well known and can be deduced either by inspecting the extended Hamiltonian

H qc =

ρ

yt2

2

+ P ( yx) −

h

2

24

 P  ( y x )( y xx )2

 a˜ 2m h2m P 2m ( y xx )2m . 2m − 1

m =1

Unfortunately, Eq. (8) has two flaws

 

.

and a˜ 2m = (2m − 1)a2m .

Thus

(7)



which unlike its discrete antecedent is no longer bounded from below, or by inspection of the linearized dispersion. Actually this fact is by now well known, and for P  (s) ∼ const . the resulting equation of motion is referred to as the “bad” Boussinesq equation [14,15]. Originally it has already emerged in the very first analytical treatment of the FPU problem by Kruskal wherein P was a cubic polynomial and the nonlinearity was assumed to be weak. This resulted in P  ∼ const . rendering the leading fourth order term linear. Kruskal circumvented the ill-posedness by deriving the one sided KdV equation as a leading asymptotic approximation of the weakly nonlinear regime. However, surprisingly enough, in the present genuinely nonlinear case there is no small parameter. In fact, the scaling which eliminates h in (8) eliminates h in all higher order terms of the expansion as well! The size of h is thus irrelevant and the problem has no genuinely small parameter. The termination of the expansion at any level beyond the strict continuum, as done for example in Eq. (6), can be judged only a posteriori according to its utility in describing the sought after phenomena. That being the case, one may turn things around and rather than to terminate the expansion at a given power of h, to terminate it at a given level of complexity, which in the present context will mean that no derivatives of order higher than the second will be kept in the Hamiltonian. Consequently, the resulting equation of motion shall have no derivatives of order higher than fourth. To carry out this program we expand the potential. Integration by parts and elimination of null divergences begets



1



2m y xx P 2m

 1 ∂  ( y xx )2m P 2m ∂ x y xx

(9)

(12)

which enables to represent the resulting equation of motion in the form

  1 ∂  ytt = ∂∂x P  ( y x ) − ∂∂x a˜ 2m h2m ( y xx )2m P 2m . ∂ x y xx 2m

dx

(11)

Rescaling à la (7), we have in terms of u = y x

(13)

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

  2 2 1 m 2m ∂ ˜ u τ τ = ∂∂z2 P  (u ) − ∂∂z2 12 ) a ( u ) P ( . 2m z 2m uz ∂ z

89

(14)

2m

Since in practice we rarely, if ever, use more than the first two terms in the sum, we write those explicitly. P 2 and P 4 ≥ 0 are assumed to be smooth





u τ τ = ∂∂z2 P  (u ) + ∂∂z2 P 2 ∂∂z P 2 u z   3 ∂2  1/4 2 3/4 ∂ u ) P P u ( . + z z 2 4 ∂z 4 10 ∂ z 2

2



(15)

Of the two aforementioned flaws we have so far focused on the second. Yet the ill-posedness remains. Since a4 < 0 the new term adds to the underlying Hamiltonian an additional negative part. This guides us to look at the one sided rendition of (15) which, however, unlike its weakly nonlinear counter part, is not asymptotic. We have

u τ = ∂∂z P  (u ) + ∂∂z



P 2 ∂∂z





P 2uz

   3/4 1/4 + γ 2 ∂∂z (u z )2 P 4 ∂∂z P 4 u z .

(16)

To quantify the impact of the newly added term besides its strict application to the lattice dynamics, we have replaced the numerical factor 3/10 in (16) with a free parameter γ 2 . The Hamiltonian structure apart from conservation of u, u 2 and energy, assures an easy access to traveling waves. In fact, if r = z − λt, then after one integration

−λu + P  (u ) + 2

2

+ γ (u r )



P 2 ∂∂r

3/4 P 4 ∂∂r





P 2 ur

1/4 P 4 ur



= P  (0).

(17)

An additional integration yields a biquadratic entity

γ 2 P 4 (ur )4 + 2P 2 (ur )2 −  = 0,

(18)

where  = 2λu 2 − 4 P˜ , P˜ = P (u ) − P  (0)u. Its solution is thus

u r2

=

P2 +





√

(19)

.

P 22 + γ 2 P 4 

2.1. Two examples of lattice potentials We evaluate the impact of the improved PDE approximation on two canonical lattice potentials:

/4. Example 1. P ( S ) = S 4√ √ Setting ( z, τ ) → ( z 3/2, τ 3/2 ) we have





2 u τ = ∂∂z u 3 + ∂∂z u ∂∂z2 u 2 + γ12 (u z )2 u zz ,

(20)

where γ12 = 8γ 2 /3 and γ12 = 4/5 for the lattice. Note that in the quartic case only two √ terms survive in the first row of the expansion (10). Let u = 2λ v. Then the traveling waves are given via

v r2 =

v (1 − v 2 ) 2v +

√ 1

where 1 = γ12 + (4 − γ12 ) v 2 .

(21)

The critical points ofthe system are at v = 0 and v = 1. The radical vanishes at v = γ1 / 4 − γ12 > 1 thus its zero need not concern us.

γ1 = 0 and γ1 = 2 we have explicit compacton solutions r  √ γ1 = 0, u = 2λ cos H (π − |r |) For

2

and

Fig. 1. λ = 1 Traveling waves. The main contribution of the added force is to further smooth the tails; now at the tails u ∼ x2 .

 2

γ1 = 2, u = 2λ cos

r

√

2 2

 H (π

√

2 − |r |)

(where the first represents the compacton via the conventional approach) which, respectively, exhibit quite clearly the smoothing effect of the improved approximation. Note that √ whereas com2, mainly due to pacton’s width has increased by a factor of extension of the tails there is a mass increase of only ∼11 percent. For other values of γ1 , compacton’s profile has to be drawn numerically. At compacton’s edge where v γ1 the tail decays as v ∼ r 2 /4γ1 leading to one additional degree of smoothness. Apart from a more faithful rendition of wave’s profile, the increased smoothness at compacton’s edge makes also their numerical treatment much easier. Since u ∼ s2 at the edge of the compacton all terms in the Hamiltonian are now continuous as opposed to u ∼ s at the edge of the original compacton. This provides us additional tools of analysis not available for the original case. In Fig. 1 we present the standard compacton solution vs. its improved approximation rendition compared with the exact solution of the lattice. The profiles of discretons were obtained via an iterative algorithm presented in [13], see the Appendix for details. It is clearly seen that the main impact of the extension is to bring the continuum rendition into a better proximity with its discrete antecedent. We demonstrate the dynamical properties of the PDE approximation by presenting several numerical simulations of Eq. (20) in Figs. 2–4. In Fig. 2 we display the emergence of improved ap-

90

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

Fig. 2. Eq. (20), γ12 = 4/5. Simulation of evolution of compactons which emerge out of the initial datum u ( z, 0) = cos2 ( 4z ) H (2π − | z|).

Fig. 5. Eq. (3). Numerical simulation of interaction of lattice discretons with a quartic potential. Discretons profiles were calculated via the iterative algorithm, [13], presented in the Appendix.

Example 2. The Hertz potential. Assume a chain of beads with their interaction potential being

P ( yn+1 − yn ) = (δ − yn+1 + yn )5/2

(22)

where δ is the pre-loading depth [6]. In a strongly nonlinear regime wherein δ | yn+1 − yn |, δ has been, as a rule, neglected. Following that route one assumes P (s) = s5/2 but then the expansion cannot be extended beyond second derivative, as with third derivative and on, all derivatives become unbounded at the origin. The difficulty may be resolved by employing the original potential P (s) = P 0 (δ + s)5/2 which removes the spurious singularity at the origin and leaves all derivatives intact. Adopting this route we assume

Fig. 3. Eq. (20), γ = 4/5. Simulation of interaction of improved approximation compactons, due to a quartic potential. 2 1

P (u ) =

2 5

(δ + u )5/2 − δ 5/2



where u = − y x ,

(23)

to obtain

 6 2 2 2 u τ τ = ∂∂z2 (δ + u )3/2 + ∂∂z2 (δ + u )1/4 ∂∂z2 (δ + u )5/4 5  9 ∂2  2 −9/8 ∂ 2 5/8 + + − u ) u u ( . (δ ) (δ ) z 2 2 ∂z 50 ∂ z

(24)

As in the quartic example, the last term is new. Traveling waves −3 are given via (19). (Since now P 4 = 8(δ+ < 0, we shall have to u )3/2 assure that the radical remains positive. Use of the simplified Hertz potential ∼u 5/2 would have caused the radical to turn negative for small u.) It is convenient to define

u = λ2 v

δ = λ2 δ1 ,

and

then

Fig. 4. Eq. (20), γ = 0. Simulation of interaction of conventionally derived compactons. Note the much larger residuum left after the interaction as compared with Fig. 3. 2 1

proximation compactons out of an initial excitation and in Fig. 3 interaction of two such compactons. To compare with Fig. 3, in Fig. 4 we present interaction of two standard compactons. Note that whereas Fig. 5 describes interaction on the original lattice (3), Figs. 2–4 refer to the dynamics on its one sided quasi-continuous rendition.

v r2 =

∗ =

4[ v 2 − 2 P˜ ( v , δ1 )] 3( v + δ1 )1/2 [1 + 27 25

5/2

−

4δ1

√

∗ ]

where

3/2

+ 10δ1 v + 5v 2 . 50( v + δ1 )5/2

(25)

Using (25) for v δ1 , ∗ → 1 thus at the tails we have

 v r2

1/2

2 − 3δ1

1/2

3δ1



v2

 ⇒

2

u λ exp −r



2λ 3δ 1/2

 −1 ,

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

Fig. 6. λ = 1 Traveling waves on the lattice corresponding to the Hertz potential (23) and varying values of δ .

implying that for every δ there is a minimal λ for which there is a traveling wave solution on a zero background. Though this speed limitation may seem as an artifact due to the additional term in the expansion, we note that a similar lower bound on λ was also observed for the original lattice: the iterative algorithm failed to converge in this regime into a traveling wave profile. Note that unlike Example 1 where monomial potential was used, there is no symmetry which allows to scale away the speed and thus render a universal profile as done for instance in [13]. In fact, for a given 0 < δ , the solitary patterns that emerge in simulations of (3) with small generic initial data, exhibit non-stationary tails and thus are not exact traveling waves. This regime warrants a separate treatment which will not be undertaken here. Though formally now the tails decay exponentially, the decay rate depends crucially on the pre-loading depth which effectively cuts the tails off at the rate of δ 1/4 . Figs. 6 and 7 display both exact lattice solutions and improved approximation profiles. 3. Strategies to restore well-posedness Trying to keep complexity in check has its price; our improved equation rendered a bit ‘nicer’ compacton profiles, but as its predecessor being ill posed it cannot be used for direct simulations. In fact, in earlier studies, cf., [8–10] the derived PDE served as a vehicle to derive an approximate form of the discreton and was then used as an input in the study of the original lattice. In what follows we discuss three strategies to induce well posedness in the PDE models of the lattice.

91

Fig. 7. λ = 1 Traveling waves due to the Hertz potential (23) and δ = 0.01. The discrete lattice solution is compared with the solutions of Eq. (24) and the solitary solution of the RB Eq. (33).

3.1. A sixth order PDE The simplest path to well posedness is to pull the first term in the second row of the expanded potential (10). Thus, ρ = 1,

H cont =

  2 yt 2

+ P ( y x ) − |a2 |h2 P  ( y xx )2



 + h4 − |a4 | P  ( y xx )4 + b2 P  ( y 3x )2 dx

(26)

where b2 = 61! . Now the Hamiltonian is bounded from below and the leading order term has the ‘right’ sign and begets a sixth order PDE as it adds to the RHS of Eq. (15)

+

 1 ∂2  2 P 3 (u zz )2 + 2 ∂∂z2 ( P 2 u zz ) . 2 ∂ z 5 Specifically, for the quartic potential we have

 2  4 2 2 u τ τ = ∂∂z2 u 3 + ∂∂z2 u ∂∂z2 u 2 + (u z )2 u zz + 5   8 ∂2  2 2 ∂2 + u u zz . 2 u (u zz ) + ∂ z2 ∂ z 15

(27)

The new equation has all the desired properties but one: it is essentially intractable neither as a PDE nor on the TW level rendering it more or less useless. Clearly, simplicity won’t do here.

92

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

3.2. A positive definite Hamiltonian approximation For a simple potential like the quartic one may attempt to overcome the well posedness issue by approximating the original Hamiltonian by a non-negative quadratic like entity

H cont =

  2 yt 2



1

+

y 2x −

4

h2 4

y 2xx

2  dx.

(28)

However, not only the leading neglected term is larger than the last one retained, it also has an opposite sign rendering the approximation useless (a similar phenomenon happens if one uses a three terms approximation). In fact the resulting entity does not support compact solitary waves but a compact kink. 3.3. Alternative representation In 1-D when no site forces are involved, there is another way to overcome the ill-posedness, [12]. To this end we use the bond length U (x, t ): U (x, t ) = Y x (x, t ) = (un+1 (t ) − un (t ))/h as the dependent variable. In terms of Y the original Hamiltonian may be rewritten as, [12]

H=

  1 2

Y t LR Y t



 + P (Y x ) dx

(29)

with its Fourier symbol Lˆ R (k), ∂ → ik,

.

1 ˆ Lˆ − R (k) = L D (k) =

4 sin2 (kh/2) h 2 k2

= 1 − β k2 +

2β 2 5

k4 + ...,

(30)

and the resulting equation of motion, at least formally, preserves the potential’s original form and looks far less daunting

2β 2 5

∂x4 )( P  (U ))xx + O (h6 ).

(31)

As before, retaining only the leading term beyond continuum begets an ill posed entity whereas engaging the O (h4 ) term begets an insurmountable sixth order PDE (though in a more appealing form than Eq. (27)) which even when slightly simplified via the β approximation Lˆ D (k) ∼ (1 − 2 k2 )2 begets







U τ τ = (1 + ∂z2 )2 P  (U )

zz

where (x, t )/ β/2 → ( z, τ ),

(32)

which in spite of its appealing well posed form is analytically intractable. Again, we are trapped within opposite constraints, tractability vs. well posedness. However, now we may escape the trap by shifting our efforts into the left part of the PDE: we approximate Lˆ D (k) ∼ 1/(1 + β k2 ) and invert the resulting bounded operator to obtain a Regularized Boussinesq equation

RB :





U tt = P  (U )

xx

Though the focus in this work has been on strict anharmonic lattices, in general a strongly nonlinear regime is accompanied by a small harmonic part. This makes the RB model even more relevant. We can now examine how our two study cases are resolved by Eq. (33). The solitary waves are given via

−λ2 U 2 + 2 Pˆ (U ) + λ2 β(U  )2 = 0 where s = x ± λt and Pˆ = P (U ) − − P (0). Assuming Pˆ (U ) = solution of (34) is given by P  (0)U

 U=

U tt = L D ( D )( P  (U ))xx

= (1 + β∂x2 +

Fig. 8. The discrete solitary wave due to the quartic potential vs. the solitary profiles due to Eq. (35) (compare with Fig. 1). β = (m)β0 and β = βm yield essentially the same profiles.

+ β U xxtt

√

mλ



2 cosh(s/ β α

which is both simple, tractable and well posed. The RB equation has two major advantages over the other approaches:

Yet the present form comes with its own price; As now the discrete effects enter via a linear force, consequently, rather than a compact support they yield the usual exponential tails which overshoot the doubly exponential decay of their discrete antecedents.

2)

, α=

2 m−2

(34) then the

(35)

.

For the Hertz potential, when 0 < δ it follows from (34) that as in (25) where finite δ implied an exponential tail, solitary waves are possible provided that 3δ 1/2 /2 < λ2 . Luckily we still have ‘one rabbit in the hat’. In Eq. (34) we have used the β = β0 h2 = h2 /12 value which assures a precise low-k approximation, corresponding to the first two terms of expansion. This means that for larger k’s our simple Pade will deviate from the exact expression, missing by a factor of 3 as k → ∞. Though not concerned with that limit, it may be beneficial to choose β in such a way that it better follows the exact dispersion for a wider range of k’s. For instance, to force an agreement with the exact dispersion at kh = π /2 or kh = π (which would call for β = 1.1366β0 h2 and β = 1.7841β0 h2 , respectively). Keeping the solitons affairs in mind we determine an effective βm minimizing in L 1 sense the difference between the solitary solution (35) and its discrete antecedent. This was done for m = 5/2, 3, 4, 6 and the respective values are displayed in the table

(33)

1) It can address both two-sided propagation and head on collisions of solitary waves. 2) On a finite interval unlike equations (8) and (15), it calls for one rather than two boundary conditions on each side, which is more appropriate physically.

α

√

1 m U m

m

5/2

3

4

6

βm /β0 (m)

1.3488 1.3540

1.6308 1.5811

2.0268 1.9149

2.3856 2.4083

Perhaps surprisingly, we have found that the simple interpolation relation



(m) =

m−

1 m−1

 βm /β0

is a remarkably good approximation. Figs. 8 and 9 display the solitary profiles based on both the ‘raw’ and the modified βm . Since the profiles based on β = (m)β0 are essentially the same as with β = βm only the later are shown. As larger m represents stronger nonlinearity, it is not surprising that βm grows with m. It is thus

P. Rosenau, A. Zilburg / Physics Letters A 381 (2017) 87–93

93

Setting h = 1, in the monomial case, the solitary waves of (3) are solutions of

u  (x) =

1

λ2



u p (x − 1) − 2u p (x) + u p (x + 1) .

(A.1)

For sufficiently quickly decaying solutions this is equivalent by the Fourier transform and its inverse to

u (x) =

1

λ2

( ∗ u p )(x)

(A.2)

where (x) is the tent function. Starting with u 0 (x) = (x), the algorithm consists of carrying the following iterations, m = 1, 2, . . . : p

u˜ m (x) = ( ∗ um−1 )(x),

⎛

Cm = ⎝ Fig. 9. The discrete solitary wave due to the ‘raw’ δ = 0-Hertz potential vs. the solitary profiles due to Eq. (35). β = (m)β0 and β = βm yield essentially the same profiles.

natural that in the Hertz case wherein the nonlinearity is relatively weak, the modification of the ‘raw’ β is quite small. 4. Closing comments We have reconsidered the Newtonian lattice, aiming to improve its modeling via a PDE in a strongly nonlinear regime. Two paths were chosen. In the first, the standard expansion was extended by two rather than one term beyond the strict continuum without an increase in the complexity of the resulting PDE as would have been the case had we followed the usual path which would have resulted in a beastly sixth order PDE. The added term improves significantly the correspondence between the resulting compacton and its discrete antecedent at the tails and turns it into a classical solution amenable, for instance, to stability analysis via its Hamiltonian. The extended model was exemplified via two canonical cases: the quartic and the Hertzian potentials. The second case being of both practical and fundamental interest, as it is a non-analytical function of the stress. To bypass this difficulty we have retained the pre-loading depth into the potential which typically is ignored. The resulting potential being analytical admits expansion in gradients. It reveals a small amplitude regime which does not support solitary waves. The second path chosen resulted in a tractable well posed model. Its main drawback is its limitation to 1-D lattices and ‘fat tails’; whereas the compactons approximated the doubleexponential decay of their discrete antecedents via compactness of their support, the solitary waves decay at the usual exponential rate and thus overshoot their antecedents. In conclusion, though we made another step forward, the issue of modeling the Newtonian lattices by PDEs is still far from coming to its end. Acknowledgements The work of the first author was supported in part by the Bauer-Neuman Chair in Applied Mathematics and Theoretical Mechanics. Appendix A. Direct computation of traveling waves 1 For a monomial P (s) = p + s p +1 lattice potential, the solitary 1 profiles may be computed via an iterative algorithm, [13]. We restate the derivation with a minor but a necessary generalization.

⎞⎛

∞

um−1 (x)dx⎠ ⎝

−∞

∞

⎞ −1 u˜ m (x)dx⎠

,

−∞

um (x) = C m u˜ m (x).

(A.3)

Now assume that 0 < u ∗ (x) = lim um (x) is a uniform limit and

!∞

−∞ u ∗ (x)dx = 1. Then C m

m→∞

→ C ∗ and u ∗ is a solution of

p

u ∗ (x) = C ∗ ( ∗ u ∗ )(x).

(A.4)

By scaling

u (x) = (C ∗ λ2 )1/( p −1) u ∗ (x)

(A.5)

we recover the solution of (A.2). Let this algorithm to potentials of the form P (s) = us extend

1 (s + δ) p +1 − δ p +1 , with δ 1 a constant. The traveling wave p +1 profile is then given by

u (x) =

1 

λ2



  ∗ (u + δ) p − δ p (x)

(A.6)

where in order to apply the Fourier transform δ p has been twice added and subtracted. As there is no scaling symmetry we replace the first step of the algorithm by

dm = (C m−1 λ2 )−1/( p −1) ,





u˜ m (x) =  ∗ [um−1 + dm δ] p − [dm δ] p



(x)

(A.7)

and set C 0 = 1. As in the monomial case, assuming that the algorithm converges, we recover via the scaling (A.5) a solution of (A.6). References [1] D. Campbell, S. Flachand, Y. Kivshar, Phys. Today 57 (1) (2004) 43. [2] P. Rosenau, Not. Am. Math. Soc. 52 (2005) 738; P. Rosenau, J.M. Hyman, Phys. Rev. Lett. 70 (1993) 564; P. Rosenau, J.M. Hyman, M. Staley, Phys. Rev. Lett. 98 (2007) 024101; P. Rosenau, Phys. Rev. Lett. 73 (1994) 1737. [3] S. Flach, Phys. Rev. E 50 (1994) 3134. [4] S. Flach, Phys. Rev. E 51 (1995) 1503. [5] S. Flach, C.R. Willis, Phys. Rep. 295 (1998) 181. [6] V. Nesterenko, J. Appl. Mech. Tech. Phys. 5 (1983) 733; V. Nesterenko, Dynamics of Heterogeneous Materials, Springer, 2001. [7] P. Rosenau, S. Schochet, Phys. Rev. Lett. 94 (2005) 045503. [8] R. Rosenau, A. Pikovsky, Phys. Rev. E 89 (2014) 022924. [9] K. Ahnert, A. Pikovsky, Chaos 18 (3) (2008) 037118. [10] K. Ahnert, A. Pikovsky, Phys. Rev. E 79 (2009) 026209. [11] G. James, C. R. Acad. Sci. Paris Sér. I Math. 332 (2001) 581. [12] P. Rosenau, Phys. Lett. A 118 (1986) 222; P. Rosenau, Phys. Lett. A 311 (2003) 39. [13] J.M. English, R.L. Pego, Proc. Am. Math. Soc. 133 (6) (2005) 1763–1768. [14] J. Boussinesq, J. Math. Pures Appl. (2) 17 (1872) 55–108. [15] S.K. Turitsyn, Phys. Rev. E 47 (2) (1993) 796–799.