Asymptotic approximations to travelling waves in the diatomic Fermi-Pasta-Ulam lattice

Abstract: We construct high-order approximate travelling waves solutions of the diatomic FermiPasta-Ulam lattice using asymptotic techniques which are valid for arbitrary mass ratios. Separately small amplitude ansatzs are made for the motion of the lighter and heavier particles, which are coupled The Fredholm alternative is used to derive consistency conditions, whose solution generates small amplitude expansions for both sets of particles.


Introduction
In this paper, we consider the effects that variations in the masses of nodes in the FPU lattice has on travelling waves by considering a diatomic FPU system in which lighter and heavier masses alternate. The original monatomic version of this lattice was simulated numerically simulated by Fermi, Pasta and Ulam [7], who discovered recurrence of the initial data. Similar recurrence behaviour was found in the Korteweg-de-Vries system and explained, by Zabusky and Kruskal [28], using solitons, which prompted investigations into solitary travelling waves in FPU [8].
The Klein-Gordon lattice (KG), nearest-neighbour interactions are linear and the nonlinearity is due introduced via an onsite potential, for example H = j 1 2 m j dq j dt 2 + ε(q j+1 − q j ) 2 + W(q j ). (1.1) Travelling waves solutions have been studied in diatomic KG lattices for some time, often in the presence of second neighbour interactions [22,26]. In contrast, the FPU lattice has no onsite potential, instead, each node experiences nonlinear interactions with its nearest neighbours via the Hamiltonian which gives rise to the equations of motion The standard monatomic case is given by m j = m for all j; however, here we consider the case where m j = m for odd j and m j = M for even j with M > m, giving a mass ratio µ = m/M < 1. The distinction between the two masses is indicated in the notation of Q j for the larger particles' displacements, and q j for the smaller, as shown in Figure 1. Here, we retain the restriction to nearest-neighbour interactions, acknowledging that many diatomic systems will have second-neighbour interactions, which significantly complicate the permitted dynamics [23]. Figure 1. Illustration of the one-dimensional diatomic FPU lattice considered here. The Q j (t) quantities represent the displacement of the heavier particles from their equilibrium positions for even values of j, and q j (t) denotes displacement of the lighter particles, at odd j.
Early interest in the behaviour of diatomic lattices was initiated by the work of Collins [5], who extended his earlier improved quasi-continuum expansions from standard FPU lattice [3,4] to the diatomic FPU case, again finding solitary waves which have the characteristic sech 2 profile at leading order. Pnevmatikos et al. [15] studied a generalised FPU-type lattice which had second-neighbour (SNI) as well as nearest-neighbour interactions (NNI). They used an ad hoc decoupling ansatz to reduce the system of coupled equations to a single generalised Boussinesq equation for the leading-order behaviour of the travelling solitary wave solution. One effect of SNI is that for certain parameter values, subsonic solitary travelling waves are predicted [25]. Pnevmatikos et al. [15] also obtain breather solutions of the FPU system with SNI.
Huang [11] investigated the diatomic FPU-KG system with symmetric nonlinearities, that is, only cubic nonlinearities in both the onsite potential and the NNI. They found breather modes in both optical and acoustic cases with forms governed by the NLS equation and travelling waves whose shape was determined by mKdV. A diatomic FPU system with no quadratic nonlinearity has been studied by Qin [18], who sought periodic travelling waves, and obtained the form of the wave for the lighter and heavier masses. Periodic travelling wave solutions of the diatomic FPU have been constructed by Betti and Pelinovsky [1] using analytic continuation from the limit of zero mass ratio. They also show that the larger wavelength solutions are stable. Vainchtein et al. [23] consider the case of highly different masses in the diatomic FPU finding a bifurcation from the simple single-humped sech 2 -like profile seen in monatomic FPU to a broader wave which supports multiple oscillations on its peak, which may be of small amplitude, and may have an amplitude similar to that of the wave itself. Faver and Wright [6] investigate the form of travelling wave solutions of a diatomic FPU lattice, and find that the underlying problem is singularly perturbed, leading to the presence of background small amplitude periodic wave upon which the solitary wave is superimposed. This system is studied further by Hoffman and Wright [10], in the limit of small mass ratios. They only find solutions when the mass ratio lies within certain ranges. Ponson et al. [16] analyse the propagation and scattering of waves in disordered diatomic FPU systems, as a model of granular crystals. They find that in disordered systems, the wave amplitude undergoes power law decay. Porubov and Andrianov [17] consider diatomic FPU lattice with both quadratic and cubic nonlinearities and derive generalised coupled Boussinesq equations for the mean displacements and differences in displacements, which describe a single system in which waves can travel in either direction.
Rigorous results on the FPU system with (periodically-arranged) polyatomic masses have been derived by Chirilus-Bruckner et al. [2] and Gaison et al. [9]. In the former work, the authors establish that the leading order asymptotic solution is valid for extremely long time intervals; the leading order solutions being those generated by the KdV approximation for travelling waves and for breather solutions, the NLS. Similar results are derived by Gaison et al. [9] who note that for arbitrary initial data, the solution decomposes into waves travelling in both directions. In both papers, corrections terms are shown to be small but are not determined.
To approximate travelling waves in this diatomic FPU system we assume separate profiles for the displacements of the heavy and light particles. The amplitude of the displacement is used as a small parameter; at each order of magnitude we find a pair of governing equations which are simplified using asymptotic expansions and solved at each order of magnitude in sequence using the Fredholm alternative.
The travelling waves discussed in this work are approximations and the system may not support exact moving solutions. Previous work, such as [2,5,9,11,15] has focused primarily on determining the form of the leading-order term; here, we aim to find the first correction terms as well, with the aim of explaining the oscillations observed in the numerical work of Vainchtein et al. [23]. A beyond-allorders asymptotic expansion, akin to the analysis of Segur and Kruskal [21], may show that no such exact solutions exist. However, numerical solutions of discrete systems suggest that these solutions decay extremely slowly and hence are relevant in understanding energy transport in applications. An analysis of this type has been carried out in the case of a diatomic Toda lattice with an asymptotically small mass ratio. Lustri and Porter [13] find 'nanoptera', that is oscillations in the tail away of the solitary wave, which are finite in amplitude but exponentially small.
In Section 2 we formulate the problem, defining variables for the masses, their displacements and interaction potentials, before reformulating the system in terms of differences in displacements, in which the system of differential equations has a form more closer to the classical travelling wave pde. Further transformations of variables are introduced in Section 3, where the problem of seeking travelling waves is posed as a system of coupled differential-delay equations. The two cases of quadratic and cubic nonlinear forces are analysed in detail in Sections 4 and 5. Finally, in Section 6, we compare the results of the two cases with each other and draw comparisons with the results of other authors.

Formulation of problem
The Hamiltonian of our diatomic FPU chain is given by where V(·) is the potential due to interactions between nearest neighbour particles via nonlinear springs.
The light particles all have mass m, whilst the heavy particles have mass M; we assume that m < M and define µ = m/M. The quantities P 2 j and p 2 j+1 denote the momenta conjugate to the displacements Q 2 j and q 2 j+1 respectively. Using Hamilton's equations, we obtain Here, we assume the potential function and interparticle forces are given by where in this case the argument φ will be the difference of the displacements q j+1 and Q j (or vice versa). Analysis of the monatomic system (1.3) is simplified by using the difference in displacements φ j (t) = q j+1 (t) − q j (t) rather than the displacements (q j ) themselves, since this leads to a governing equation which is more clearly a discrete nonlinear wave equation, namely . We follow this procedure in the diatomic system, introducing where φ j and ψ j now describe the difference displacements between neighbouring light and heavy particles, (and heavy and light particles, respectively). We assume that even nodes are occupied by the heavier masses, so that Q j and φ j are defined for even j, and the odd nodes host the lighter masses (that is, q j and ψ j correspond to odd j). Under the reformation (2.5), the governing equations (2.2)-(2.3) become with the interaction potential V(·) being given by (2.4).

Travelling solitary waves
To investigate travelling waves, we introduce the spatial shift operator e ±∂ , where e ±∂ f ( j) = f ( j±1). We also define the mass ratio µ = m/M with 0 < µ < 1, rescale time via t = t √ µM and drop the tilde, so the system of equations (2.6)-(2.7) are rewritten as so that ξ = 1 2 (φ + ψ) and ζ = 1 2 (φ − ψ). Since both φ and ψ are differences in displacement, ξ is formally a mean difference in displacements and ζ is a second difference. By adding and subtracting (3.1)-(3.2) we obtain To proceed further, we need to scale both independent ( j, t) and dependent (ξ, ζ) variables, and these scalings differ, depending on whether a = 0 or a 0, so we consider each case separately. We start with the quadratic case a 0 and, for simplicity, assume b = 0. Later we consider the cubic case a = 0, b 0.
Following the standard method for deriving a KdV equation from the FPU lattice equations, we introduce h 1 and write At leading order, we expect φ = ψ = O(h 2 ), and hence the scaling for the second difference in displacement, ζ, needs to be smaller than that of the mean displacement, ξ. To ensure as many terms as possibly balance in the equation for v, ζ has to be just one power of h smaller than ξ.
We assume that the functions ξ and ζ are slowly-varying in j, and noting that ∂ = ∂/∂ j and ∂ y = ∂/∂y, we replace the discrete differences by their continuum limits The governing equations (3.4)-(3.5) thus become From the leading order terms in (4.5) we have v = 1 2 u y (1 − µ)/(1 + µ), and from (4.4) we have C 1 u ττ = C 2 u yy + C 3 v y for suitable constants C j . Thus the system has the form of a travelling wave equation.
Transforming to a moving wave frame, via z = y − c 0 τ, and T = h 2 τ, (4.6) implies ∂ y = ∂ z and ∂ τ = h 2 ∂ T − c 0 ∂ z . After expanding terms, integrating the equation for u (4.4) once with respect to z, moving all O(1) terms to the lhs, and O(h 2 ) to the right, and neglecting higher order terms, we find where Equations (4.7)-(4.8) have the form Aw = h 2 g where h 1 and we are seeking an O(1) solution for w. At leading order, the problem reduces to Aw = 0, which only has an O(1) solution for w if A is singular. Hence a condition for a solution to exist is that we choose the currently undetermined speed c 0 so that det(A) = 0. This requirement means that is the speed of sound in the (y, τ) coordinate system, which corresponds to the speed of sound in the original lattice being in ( j, t) coordinates. To obtain expressions for u(z, T ) and v(z, T ) we need to consider the next order terms, which requires more careful analysis, including an application of the Fredholm alternative [14].

Leading-order analysis-the quadratic case
The singular matrix, A, has a kernel, span{k}, range, span{r}, and normal to the range, n, given by Hence the solution of Aw = 0 is w = w(z, T )k for any function, w(z, T ).
In order for there to be a solution of the perturbed problem Aw = h 2 g 0 we have an additional consistency condition, this time on the rhs vector, g. Since A maps all R 2 onto the line λr the rhs vector, g must lie on this line, hence in (4.8) we require g 2 = 2g 1 (1 + µ)/(1 − µ). As the vector n (4.11) is normal to the range, the condition can be written as n · g = 0 or, from (4.8) Since the problem (4.7) has two components, the O(h 2 ) correction term should also have two components, which we write as coefficients of the two independent vectors r and k. For 0 < µ < 1, the vectors r and k are never parallel (or antiparallel) so, in general, we can write the solution of and since w = (u z , v) T , we have (4.14) A component of the correction term in the k direction can simply be absorbed into the leading order solution, which we have already discussed. So, here, we only consider the rôle of the component w z in the r direction. Using the leading order expressions u = 2(1 + µ)w, v = (1 − µ)w z and (4.9), simplifying (4.12) yields This gives the leading order description of u, and hence also ξ, and φ, ψ.

Higher-order analysis-the quadratic case
Inserting (4.13) into Aw = h 2 g, we now aim to solve w z Ar = g, and so find the higher-order correction terms w. We have so substituting the solution (4.14) into (4.7)-(4.8), yields w(z, T ) = (1 − µ)(5 + 6µ + 5µ 2 )w zz 6(1 + µ)(3 + 10µ + 3µ 2 ) . . The ξ 0 , ζ 0 curves (thicker dashed lines) correspond to the leading order terms, that is, neglecting the w terms. The higher order approximations, ξ, ζ are shown in narrower solid curves; the larger amplitude curves, which have the typical sech 2 profile correspond to ξ 0 , ξ, whilst ζ 0 , ζ are given by the smaller amplitude curves which are more oscillatory and have zero mean. Hardly any difference between ξ and ξ 0 or between ζ and ζ 0 can be seen in the case h = 0.1 (upper panel). However, when h = 0.25 (lower panel), the development of more complex structure can be seen, which includes more oscillatory behaviour in ζ than in ζ 0 . All plots correspond to a mass ratio of µ = 0.02.
Inverting the transformations (4.14), (4.1), and (2.5), we obtain The functions ξ, ζ are illustrated in Figure 2, which shows both the leading-order solution (ξ 0 , ζ 0 ) arising from w = w z k in (4.13), and the solution w = w z k + h 2 w z r which includes the first corrections terms, leading to (ξ, ζ). We note that when h = 0.1, the curves for ζ 0 and ζ are almost identical, as are the profiles of ξ and ξ 0 . Furthermore, ζ is small and has the same shape as ξ , thus the combinations ξ(y, τ) ± ζ(y, τ) will be well-approximated by ξ(y ± s, τ) for some shift, s. However, at larger values of h, as illustrated in the lower panel of Figure 2, the more accurate approximation introduces a more complicated oscillation into the shape of ζ as shown by the lower solid line. The generation of such internal oscillations in the peak of the travelling wave has been observed in the numerical work of Vainchtein et al. [23], who consider the Toda lattice, V(φ) = (α/β)(e −βφ − 1) + αφ which corresponds to V (φ) = α(1−e −βφ ). Expanding this for small φ leads to V (φ) ≈ αβ(φ− 1 2 βφ 2 ), which gives the quadratic FPU system, and we note that no choice of α, β is approximated by the cubic FPU system. Small amplitude travelling waves in this diatomic Toda lattice are thus well-approximated by the quadratic diatomic FPU system.
The amplitude of the main pulse ξ 0 , given by the leading order sech 2 solution (4.16), and the perturbations due to ζ and w waves all depend on the mass ratio, µ. By introducing a shift (z = z − hδ From the solution (4.19)-(4.20) we determine a speed-amplitude relationship for the kink: In terms of the original ( j, t) variables the speed is given by c = ( c 0 + h 2 C)/ √ m, and the amplitude of the kink is So, in this system, the only slowly-varying travelling waves are supersonic. Equation (4.10) again gives the relationship between c 0 and c (4.9). Note that in the limit µ → 0 + , (4.22) gives large amplitude kinks.
5. Cubic case: a = 0, b 0 In the case, where there is no quadratic nonlinearity in the system, but there is a cubic term, that is, a = 0 and b 0, we need to use a different scaling: In place of (4.1), we define We still have the second difference of displacements, ζ, being smaller than the mean difference of displacements, ξ, by one order of magnitude in h, since we expect φ = ψ to leading order, but now φ, ψ = O(h), in place of O(h 2 ) in the quadratic case, compare (5.1) with (4.1).

(5.5)
Here, we have integrated the upper equation with respect to z; and in both equations, collected the leading order terms onto the lhs are leading order, leaving the O(h 2 ) first correction terms on the right. The matrix A on the lhs of (5.4) is the same as that in (4.7), and we follow the same strategy as in §4.1, that is, we write the solution as (4.13).

Higher-order analysis-cubic case
We write (5.4)-(5.5) as Aw = h 2 g where A is singular and has kernel, range and normal to the range, span{k}, span{r}, n given by (4.11) respectively. We again write the solution as (4.14) and require that the rhs vector g satisfies n · g = 0. At leading order in h, this consistency condition, together with the substitutions (4.14) implies which is an mKdV equation, rather than the KdV equation (4.15), which occurs for the quadratic case. Equation (5.6) has a solitary wave solution given by Converting this back to the original variables ( j, t), we find the speed of the wave is c = (c 0 +h 2 C)/ √ m: since we require C > 0, these waves are supersonic. The amplitude of the resulting kink is which has no dependence on the speed c (at leading order). A similar cancellation occurs with the standard continuum approximation of the monatomic FPU problem, where the dependence of amplitude on speed can be recovered by use of quasicontinuum approximations [19,20,24]. We can also determine the correction term, w z : from the assumed form (4.14) we have Aw = h 2 w z λr; and from the last line of (5.4)-(5.5), we obtain 2(1+µ)λ w z = −(2c 2 +1+µ)v zz − 12b(1+µ)u 2 v + (1−µ)(6bu 2 u z + 1 3 u zzz ), hence the first correction term in (4.13) is given by where w is given by (5.7). Note that while the same ansatz (4.14) is used in both quadratic and cubic cases, the expressions for w, (4.16) and (5.7), and w, (4.18) and (5.10), differ. Inverting the changes of variables, the spring extensions φ, ψ are given by The ξ 0 , ζ 0 curves (thicker dash-dotted lines) correspond to the leading order terms, that is, neglecting the w terms. The more accurate solutions, ξ, ζ, are shown in narrower solid curves. In both panels, the largeramplitude curve corresponds to the sech-type solution ξ 0 , ξ, and the smaller amplitude, and more oscillatory curves illustrate the second differences given by ζ 0 , ζ, which have zero mean, as in Figure 2. All plots correspond to a mass ratio of µ = 0.02.
(ξ, ζ) which includes both terms on the right hand side of (4.13), that is, it includes the effect of w. This correction term is seen to increase slightly the amplitude of the oscillation in ζ, but not make a significant change in wave form as was observed in the quadratic case ( Figure 2).
The first correction term in (5.11) can be removed by applying a phase shift z → z + δ with δ = 1 2 h(1 − µ)/(1 + µ), following which, we obtain The quantity κ is the upper curve plotted in Figure 3, for both the quadratic (4.21) and cubic cases (5.13). We note that in both cases, κ takes the value zero at µ = 1 (the monatomic case), and increases away from µ = 1, and is symmetric under the transformation µ → 1/µ. The limiting values as µ → 0 + are κ 2lim = 19/72 ≈ 0.264 and κ 3lim = 5/18 ≈ 0.278. Note, however, that the definitions of φ in the two cases (4.21) and (5.13) have different signs in the coefficient of κw z z due to the differing signs in (4.18) and (5.10). In both panels, the thin solid line represents the leading order solution, and the thicker solid line shows the solution including the O(h 2 ) correction terms, the difference between the two solutions is displayed using thick dashed lines. Note the effect of the different signs in front of the w z z terms in (4.21) and (5.13). All plots correspond to a mass ratio of µ = 0.02.
In Figure 5 we show the effect of the higher-order perturbative terms on the travelling wave, after the effects of a shift in the Goldstone mode is accounted for. In both panels we plot φ, which is identical to ψ once the phase shift has been accounted for. We see that in the quadratic case (upper panel), there is the development of an oscillation on the peak of the travelling pulse, whereas in the cubic case, (lower panel) the effect of the perturbation is to 'sharpen' the pulse, increasing the amplitude and reducing the half-height width.

Conclusions
We have investigated the properties of travelling waves in a one-dimensional diatomic FPU lattice with a polynomial interaction potential. To do this, we have derived continuum limit equations in multi-component systems, and used small amplitude asymptotic expansions to convert the governing equations to a system of coupled PDEs. Whilst many authors consider the diatomic system in the limit of asymptotically small mass ratio, for example [1,10], our results, in contrast, are valid for any mass ratio. We solve the equations at each order of the asymptotic expansion by application of the Fredholm alternative. This is similar to the approach based on Bloch waves used by Chirilus-Bruckner et al. [2] in that it removes the requirement of an ad hoc ansatz that was used by Pnevmatikos et al. [15] to reduce the coupled two-component system to a single equation.
In our formulation, the leading order problem is singular, giving the wave speed-amplitude relationship. At next order we have a singular forced equation, which is subject to a consistency condition. The leading-order solution gives rise to the usual sech 2 or sech-type solutions, which have been found by many authors previously, for example, [5,11,15], and proven to be accurate approximations [2,9] for long periods of time. At next order, the correction terms that we have found here correspond to simple phase shifts and, at the following order, we find terms which give rise changes in the shape of the solitary waves. We believe these last results to be new.
Comparing the equations (4.21) for the quadratic nonlinearity with the cubic case (5.13), we note a difference in the sign of the correction term: φ = 2h(1 + µ)[w ± κh 2 w z z ], (both κ 2 and κ 3 are positive, as shown in Figure 3). In the case of a quadratic nonlinearity, the '+' sign causes the emergence of an oscillation in the peak of the solitary wave; whilst in the cubic system, the '−' sign causes a 'sharpening' of the peak and possibly the development of an oscillation in the tail of the solitary wave. All the calculations presented herein are based on small amplitude asymptotic approximations. At larger amplitudes, we expect the oscillations to grow, both in amplitude and in number-as shown in the numerical results of Vainchtein et al. [23] who investigate the diatomic Toda chain. With the Toda potential, there is typically a combination of quadratic, cubic and higher-order nonlinearities involved. Also, we note from Figure 3 that extremely small or large mass ratios will exhibit a greater effect than mass ratios close to unity-the almost monatomic case.
In future work, we hope to apply these methods to find properties of travelling waves in mass-inmass FPU systems [12], which is an alternative generalisation of the FPU system with two components in each unit cell.