Recent studies have presented first-order multiple time scale approaches for exploring amplitude-dependent plane-wave dispersion in weakly nonlinear chains and lattices characterized by cubic stiffness. These analyses have yet to assess solution stability, which requires an analysis incorporating damping. Furthermore, due to their first-order dependence, they make an implicit assumption that the cubic stiffness influences dispersion shifts to a greater degree than the quadratic stiffness, and they thus ignore quadratic shifts. This paper addresses these limitations by carrying-out higher-order, multiple scales perturbation analyses of linearly damped nonlinear monoatomic and diatomic chains. The study derives higher-order dispersion corrections informed by both quadratic and cubic stiffness and quantifies plane wave stability using evolution equations resulting from the multiple scales analysis and numerical experiments. Additionally, by reconstructing plane waves using both homogeneous and particular solutions at multiple orders, the study introduces a new interpretation of multiple scales results in which predicted waveforms are seen to exist over all space and time, constituting an invariant, multiharmonic wave of infinite extent analogous to cnoidal waves in continuous systems. Using example chains characterized by dimensionless parameters, numerical studies confirm that the spectral content of the predicted waveforms exhibits less growth/decay over time as higher-order approximations are used in defining the simulations' initial conditions. Thus, the study results suggest that the higher-order multiple scales perturbation analysis captures long-term, nonlocalized invariant plane waves, which have the potential for propagating coherent information over long distances.

Introduction

Wave propagation in linear monoatomic, diatomic, and other chains is of contemporary interest due to their nontrivial dispersion, filtering, and frequency bandgap behavior. These systems may arise in the modeling of one-dimensional (1D) waveguides, or wave propagation in three-dimensional crystals along preferred directions, such as the directions along the unit cell edges of zincblende crystals (e.g., GaAs and NaCl). Nonlinearities in these systems introduce behavior of interest to the engineering community such as amplitude-dependent dispersion and group velocity associated with weak nonlinearities and solitary wave propagation (e.g., solitons) associated with strong nonlinearities. Specific engineering technology that has been proposed from studying wave propagation in periodic structures are switches [1], filters [2], and diodes [3].

Several recent studies have analyzed amplitude-dependent dispersion (nontrivial temporal/spatial frequency relationships) and other nonlinear behavior in discrete chains. Many of such studies have applied perturbation theory to uncover amplitude-dependent phenomena. Vakakis and King [4] employed perturbation techniques to study the monocoupled nonlinear media. They identified amplitude-dependent attenuation and propagation zones, which were associated with synchronous motion of standing waves and nonsynchronous motion of traveling waves, respectively. Chakraborty and Mallik [5] obtained first-order perturbation results for a cubically nonlinear monoatomic chain, with a discussion on the bounding frequencies of the propagation zone. The interaction of oppositely traveling waves and the effect of boundary conditions for semi-infinite and finite chains were also examined. In Ref. [6], Dubus et al. applied higher-order multiple scales to a linear monoatomic chain coupled (linearly) to another monoatomic chain with quadratic nonlinearities. Tuning the linear interchain coupling stiffness shifted the magnitude of a critical wavenumber associated with resonant energy transfer between the chains. Narisetti et al. [7] and Manktelow et al. [8,9] investigated wave propagation in cubically nonlinear monoatomic and diatomic chains using Lindstedt–Poincaré and multiple scales analyses, respectively, with an emphasis on amplitude-dependent dispersion shifts. Further, they identified wave-based devices which exploit bandgap shifting to enable tunable filtering and wave-guiding. In Ref. [10], dispersion shifts of a monoatomic chain governed by quadratic and cubic interactions are obtained using a Lindstedt–Poincaré approach. Furthermore, particular solutions are obtained at the first-order only and waveform invariance is not considered. For a more comprehensive review of tunable dispersion behavior in discrete periodic systems, see Ref. [11].

The study of higher harmonics arising in these discrete systems has centered on their generation, energy exchange, and/or decay. Cabaret et al. [12] analytically and experimentally studied second harmonic generation in a precompressed weakly nonlinear diatomic chain, reporting amplitude-dependent self-induced attenuation and self-induced transparency at the band edges. Sánchez-Morcillo et al. [13] conducted a study on the propagation of spectral content of harmonically driven compressional waves in quadratically nonlinear 1D granular chains subject to a static preload. An interplay between the fundamental and second harmonic, governed by dispersion, was observed for nonevanescent modes. Swinteck et al. [14] employed second-order perturbation studies to examine the self-interaction and wave–wave interactions in a monoatomic chain with quadratically interacting masses. Molecular dynamics simulations affirmed that vibrational modes exchange energy and have lifetimes that are inherently long. These studies discuss the generation, energy exchange, and decay/lifetime of the spectral content in the chains, respectively, different from the invariant spectral content studied herein.

In the long wavelength limit, discrete media behave as if continuous, and within such an approximation, studies have identified the presence of invariant waveforms. Invariant waveforms, or specific harmonic content that remains unaltered over the course of propagation, are of interest to the engineering community because of their ability to maintain information over long distances and times. They arise due to an interplay between dispersion and nonlinearity, the former of which assigns different wave speeds for different harmonics, and the latter of which can shift such wave speeds to maintain the wave's spatial waveform [15]. Higher-order asymptotic analysis is often employed when analyzing these waveforms. Laitone [16] extended an expansion by Friedrichs [17] to higher-orders to capture the effect of vertical motions and pressure variation of finite-amplitude, long water waves. At the first-order, the solution by Friedrichs is shown to be a cnoidal wave. Cnoidal waves appear to be less-studied by the engineering community as compared to other nonlinear invariant solutions. Figure 1 provides a classical example of cnoidal water waves captured in 1933 [18]. These waves are characterized by their sharp crests formed by nontrivial spectral content that propagates without changing. It is known that spatially extended, long-crested cnoidal waves satisfy the nonlinear Kortegweg de-Vries (KdV) equation [19] and travel unaltered, much like solitons but with infinite, periodic extent (vis-à-vis compact support). Kraenkel and Manna [20] applied higher-order multiple scales perturbation theory to the KdV equations to yield solutions for solitary waves. Some studies have presented continuum approximations to nonlinear lattices and recovered the KdV equation as a result. Friesecke and Mikikits-Leitner [21] demonstrated how the nonlinear interaction of particles in a discrete chain uncovered, in the continuum limit, continuous partial differential equations, such as the KdV equation and thus cnoidal waves. The dynamics of discrete particles are not considered for the duration of their proofs. Gaison et al. [22] applied homogenization techniques to the long wavelength limit of a nonlinear monoatomic lattice with periodic material properties. They uncovered solutions that approximately satisfy a KdV equation. Friesecke and Pego [23] presented that sonic speed wave pulses in nonlinear lattices are governed by the continuum limit regardless of the length scale. The studies in Refs. [2123] provide an insight into how invariant waveforms may arise in these lattices in the long wavelength limit. However, no application of multiple scales to demonstrate invariance of waves in discrete lattices, for all wavelengths, has been found in the prior literature.

Waveform invariance has received attention in recent years, primarily focused on the propagation of solitary waves and solitons. This is in contrast to single-frequency plane waves in nonlinear media, which are known to induce higher harmonics as they propagate. Sen et al. [24] provide an overview on the phenomenon of solitary waves in discrete media, highlighting differences from their characteristics in continuous media. They also discuss ways to manipulate solitary waves in granular chains for shock mitigation. Another review of solitons in nonlinear lattices is provided by Kartashov et al. [25] with an emphasis on theoretical implications. Their existence and stability in one, two, and three dimensions are considered. Experimental studies have confirmed the existence of solitary waves and solitons in discrete chains. Molerón et al. [26] investigated a 1D chain of repelling magnets analytically and experimentally. The analysis discussed a solitary wave whose profile and speed depends on amplitude. It was shown that in the low energy limit these wave profiles take the form of a KdV soliton with a sech2 shape, whereas in the high energy limit, energy becomes further localized, and the waveform exhibits a hat function with width on the order of the lattice spacing. Daraio et al. [27] experimentally observed solitary waves propagating down a strongly nonlinear chain of Teflon spheres, obtaining a good agreement between theoretically estimated values (based on Hertz interaction: 3/2 force–displacement characteristic) of the wave speed and experimental measurements.

The stability of wave propagation in nonlinear continuous and discrete systems has also been analyzed in a number of studies. Such studies are of practical interest because they define how systems may deviate from their expected behavior or define critical thresholds for system parameters. Newton and Keller [28] analytically studied the stability of continuous periodic plane waves by examining wave amplitude. They arrive at an eigenvalue problem in which positive eigenvalues correspond to exponential growth in amplitude (unstable), and negative eigenvalues correspond to decay in amplitude (stable). Bickham et al. [29] studied the stability of localized modes in a monoatomic chain with quadratic and cubic force interactions. Numerical simulations introduced local static distortions into chains. Instabilities of the localized vibration modes began to arise in simulations in which the cubic anharmonicity increased to the point where the lowest possible plane wave frequency was reached. Further increases in the cubic anharmonicity led to behavior reminiscent of a triatomic molecule until further increases created instabilities yet again. Friesecke and Pego [30] identified exponential growth in amplitude of waves in discrete lattices by applying Floquet theory and analyzing evolution equations. In Ref. [31], Flach and Gorbach uncovered a threshold for tangent bifurcations in discrete breathers. Huang and Hu [32] investigated the stability of diatomic lattices using a quasi-discreteness approach. They derived evolution equations and fixed points for acoustic and optical modes and identified the presence of asymmetrical gap solitons. Gorbach and Johansson [33] presented findings on the linear stability of discrete gap breathers in a 1D diatomic chain with harmonic intersite potential and nonlinear external potential. Six types of instabilities were identified, which took on either oscillatory or nonoscillatory forms. An “inversion of stability” regime was found characterized by practically radiationless mobility. The authors would like to emphasize that discrete breathers studied in Refs. [3133] are spatially localized and thus in contrast to the plane waves of infinite extent studied herein.

This paper builds upon recent studies of plane waves in nonlinear monoatomic and diatomic chains and other discrete systems by developing a higher-order multiple scales procedure to inform dispersion, stability, and waveform invariance. Previous efforts have focused primarily on first-order multiple scales analysis and have not yet considered both quadratic and cubic nonlinearities, which are of interest since they both necessarily arise in Taylor Series expansions of nonlinear interactions. With the higher-order analysis applied to these discrete systems, former assumptions neglecting quadratic nonlinearities on their amplitude-dependent dispersion shifts can now be rigorously justified. The reviewed literature also has not presented an analytical approach based on multiple scales to describe the loss of plane wave stability in nonlinear monoatomic and diatomic chains. This paper formulates a local fixed point analysis based on the multiple scales-derived evolution equations, and thus provides stability information in analytical form. It also numerically simulates nondimensionalized chains to identify and characterize their amplitude-dependent stability. Waveform invariance in discrete media has generally focused on localized waveforms such as solitary waves and solitons. To the best of our knowledge, no analysis has been put forth showing that multiple scales theory, applied to nonlinear monoatomic and diatomic chains, can capture solutions which asymptotically approach periodic, invariant waveforms akin to cnoidal waves. This paper provides both analytical and numerical evidences that a higher-order multiple scales approach yields such invariant waveforms. This contrasts, for example, with other equally valid interpretations [14] in which similar higher-order perturbation analyses are used instead to provide information about multiple wave scattering and phonon lifetimes.

System Descriptions

This section introduces the two systems studied herein, namely nonlinear monoatomic and diatomic chains, defining their equations of motion to which the asymptotic analyses are applied and presenting schematics for the ease of system visualization. For the monoatomic chain, the smallest repeatable subsystem, or unit cell, consists of a single mass and its connecting springs and dampers, as presented in Fig. 2.

A force balance on the jth mass yields its equation of motion 
mẍj+k1(2xjxj+1xj1)εk2(xj+1xj)2+εk2(xj1xj)2εk3(xj+1xj)3εk3(xj1xj)3+εc(2ẋjẋj+1ẋj1)=0
(1)

where xj=x(j,t) denotes the displacement from equilibrium of the jth mass, m its mass, k1, k2, and k3 are the linear, quadratic, and cubic stiffnesses, respectively, and c denotes the linear damping coefficient. Time derivatives of xjẋj and ẍj —denote the velocity and acceleration of the jth mass, respectively. The parameter ε is assumed to be small and introduced herein as a bookkeeping device [34].

By contrast, the unit cell of the diatomic chain, depicted in Fig. 3, contains two masses coupled with identical springs and dampers.

Summing the forces on each degree-of-freedom in the jth unit cell leads to its equation of motion 
Mx¯̈j+l=j1j+1Klx¯l+l=j1j+1εClẋ¯εf¯quadraticεf¯cubic=0¯
(2)
where 
x¯j=[xa(j,t)xb(j,t)]
(3)
 
M=(ma00mb)
(4)
 
Kj1=(0k100), Kj=(2k1k1k12k1), Kj+1=(00k10)
(5)
 
Cj1=(0c00), Cj=(2ccc2c), Cj+1=(00c0)
(6)
 
f¯quadratic=[k2(xb(j1,t)xa(j,t))2+k2(xb(j,t)xa(j,t))2k2(xa(j,t)xb(j,t))2+k2(xa(j+1,t)xb(j,t))2]
(7)
 
f¯cubic=[k3(xb(j1,t)xa(j,t))3+k3(xb(j,t)xa(j,t))3k3(xa(j,t)xb(j,t))3+k3(xa(j+1,t)xb(j,t))3]
(8)

It is noted here that the presence of linear and small cubic and quadratic force–displacement relationships describes many systems whose nonlinear interactions are expanded in a Taylor series.

Analysis Approach

Method of Multiple Scales: Monoatomic Chain.

In this section, a higher-order method of multiple scales (MMS) is developed for the monoatomic chain presented in Sec. 2, extending the analysis through the second-order to capture harmonic content and time evolution of amplitude and phase. It is shown in Sec. 5.3, that the inclusion of these higher-order features increases the invariance of the cnoidal-like wave. Higher-order results also facilitate the stability and dispersion analyses, providing insight to behavior not evident at solely the first-order.

Multiple time scales are assumed to exist such that 
T0=t, T1=εt, T2=ε2t, , Tn=εnt
(9)

where t denotes the original time scale, and Tn represents the nth time scale. Since ε1, each time scale advances more slowly than its predecessor. Introducing these times scales plays an important role in detecting amplitude-dependent dispersion shifts since phase change at each slow time scale can be interpreted as dispersion corrections.

In accordance with these time scales, differentiation with respect to time contains multiple orders 
()̇=D0()+εD1()+ε2D2()++εnDn()
(10)

where Dn() denotes an operator defined to represent differentiation with respect to Tn.

It follows that 
()̈=(D0()+εD1()+ε2D2()++εnDn())2=D02()+2εD0D1()+ε2D12()+2ε2D0D2()+O(ε3)
(11)
The MMS seeks a series expansion of the solution of the form 
xj=xj(0)+εxj(1)+ε2xj(2)++εnxj(n)
(12)
where xj(n) denotes the nth-order solution. Thus, the asymptotic solution is expanded about a known zeroth-order solution xj(0). Substituting Eqs. (10)(12) into Eq. (1) and collecting matching orders of ε yields n+1 equations. The first three are presented here 
ε0: mD02xj(0)+k1(2xj(0)xj+1(0)xj1(0))=0
(13)
 
ε1: mD02xj(1)+k1(2xj(1)xj+1(1)xj1(1))=2mD0D1xj(0)cD0(2xj(0)xj+1(0)xj1(0))+k2(xj+1(0)xj(0))2k2(xj1(0)xj(0))2+k3(xj+1(0)xj(0))3+k3(xj1(0)xj(0))3
(14)
 
ε2: mD02xj(2)+k1(2xj(2)xj+1(2)xj1(2))=2mD0D2xj(0)mD12xj(0)2mD0D1xj(1)cD1(2xj(0)xj+1(0)xj1(0))cD0(2xj(1)xj+1(1)xj1(1))+fquadratic(xj1,xj,xj+1)+fcubic(xj1,xj,xj1)
(15)
where 
fquadratic(xj1,xj,xj+1)=2k2xj1(0)xj1(1)+2k2xj1(0)xj(1)+2k2xj1(1)xj(0)2k2xj+1(1)xj(0)2k2xj+1(0)xj(1)+2k2xj+1(0)xj+1(1)
(16)
 
fcubic(xj1,xj,xj+1)=3k3xj1(0)2xj1(1)3k3xj1(0)2xj(1)6k3xj1(0)xj1(1)xj(0)+6k3xj1(0)xj(0)xj(1)+3k3xj1(1)xj(0)26k3xj(0)2xj(1)+3k3xj+1(1)xj(0)2+6k3xj+1(0)xj(0)xj(1)6k3xj+1(0)xj+1(1)xj(0)3k3xj+1(0)2xj(1)+3k3xj+1(0)2xj+1(1)
(17)

Note that the left-hand sides of Eq. (13) through (15) are similar, i.e., the linear kernels of the collected equations take the same form. Also note that each order is effectively forced by the solution obtained from the previous order.

By examination of Eq. (13), it is apparent that xj(0) takes the form of a Bloch wave since its governing equation is linear with periodic coefficients 
xj(0)=12Aeiω0T0eiμj+c.c.
(18)

where A denotes the wave's amplitude, ω0 and μ are the wave's fundamental temporal and spatial frequencies, respectively, and c.c. denotes the complex conjugate of all preceding terms.

Since A is generally a complex quantity, it is useful to express it in polar form 
A=αeiβ
(19)

where α=α(T1,T2,,Tn) and β=β(T1,T2,,Tn) are real quantities. Wave amplitude is a critical parameter as it determines the strength of the nonlinearities and consequently the amount of dispersion shifting and the relevancy of the MMS results.

The zeroth-order dispersion relationship follows from the Bloch solution form Eq. (18) and the governing equation (13) 
ω0=2k1m(1cosμ)
(20)

This relationship, which stems from the linear, undamped chain (ε=0), shifts for waves at different amplitudes in the presence of nonlinearities.

With xj(0) known, xj(1) is found by updating Eq. (14) with Eq. (18). Secular terms, i.e., those that contain eiω0T0eiμj, or its complex conjugate, induce resonant responses and thus cause the series solution to become nonuniform over the time scales of interest. Consequently, they are set to zero 
0=k3α3(34cos2μ3cosμ+94)+ciω0α(1cosμ)+imω0D1αmω0αD1β
(21)
which, after separating real and imaginary parts and setting to zero, restricts the evolution of α and β with respect to T1 
D1α=γα
(22)
 
D1β=δα2
(23)
where 
γ=cm(1cosμ)
(24)
 
δ=3k34mω0(cos2μ4cosμ+3)
(25)

Equations (22) and (23) represent evolution equations describing the dynamics of the wave on a slower time scale. Phase evolution governs dispersion shifts, whereas amplitude evolution affects wave stability as will be presented in Secs. 5.1 and 5.2.1, respectively.

With secular terms removed, Eq. (14) reduces to a linear, first-order ordinary differential equation, with constant coefficients and harmonic forcing at 2ω0 and 3ω0. The homogeneous solution is ignored since it contains no information than that already given by xj(0), so a particular solution is sought of the form 
xj(1)=12B1e2iω0T0e2iμj+12C1e3iω0T0e3iμj+c.c.
(26)
where B1 and C1 can be found using the method of undetermined coefficients 
B1=ik2A2(sin2μ2sinμ)4mω02+2k1cos2μ2k1=b1A2
(27)
 
C1=k3A3(3cos2μcos3μ3cosμ+1)18mω02+4k1cos3μ4k1=c1A3
(28)

The amplitudes in Eqs. (27) and (28) describe the amount of second- and third-harmonic content that make-up the first-order plane wave solution.

Turning attention to the second-order problem, Eq. (15), with known solutions for xj(0) and xj(1), the same procedure is carried out. Removal of secular terms now determines the evolution of α and β with respect to T2 
D2α=γδω0 α3
(29)
 
D2β=Xα4+Yα2+Z
(30)
where 
X=14mω0(k3c1(9cosμ9cos2μ+3cos3μ3)+2δ2m)
(31)
 
Y=k2b1imω0(2sinμsin2μ)
(32)
 
Z=γ22ω0
(33)
The particular solution for xj(2) takes the form 
xj(2)=12B2e2iω0T0e2iμj+12C2e3iω0T0e3iμj+12E2e4iω0T0e4iμj
 
+12F2e5iω0T0e5iμj+c.c.
(34)
where 
B2=b21α4e2iβ+b22α2e2iβ
(35)
 
C2=c21α5e3iβ+c22α3e3iβ
(36)
 
E2=e2A4
(37)
 
F2=f2A5
(38)
and 
b21=12eiμ(4mω02e2iμ+k1e4iμ2k1e2iμ+k1)(3k3b1e6iμ2k2c1e6iμ6k3b1e5iμ+2k2c1e5iμ16mb1ω0δe3iμ3k3b1e4iμ+2k2c1e4iμ+12k3b1e3iμ3k3b1e2iμ2k2c1e2iμ6k3b1eiμ2k2c1eiμ+3k3b1+2k2c1)
(39)
 
b22=4icb1ω0(1cos2μ)8ib1γmω04(mω02+k1cos2μk1)
(40)
 
c21=32(k1e6iμ+9mω02e3iμ2k1e3iμ+k1)c1(16k3e3iμcos4μ+2k3e6iμk3e5iμ+16k3e3iμcos2μ+12e3iμδmω0+2k3e4iμ6k3e3iμ+2k3e2iμk3eiμ+2k3)
(41)
 
c22=36iγc1mω0+icc1ω0(12cos3μ12)+ib1k2(4sin3μ+4sin2μ+4sinμ)2(9mω02+2k1cos3μ2k1)
(42)
 
e2=b1k3(12cosμ+6cos4μ12cos3μ6)+ic1k2(8sinμ8sin4μ+8sin3μ)8(8mω02+k1cos4μk1)
(43)
 
f2=c1k3(3cos2μ6cosμ+6cos4μ3cos5μ3cos3μ+3)50mω02+4k1cos5μ4k1
(44)
These amplitudes build upon the higher harmonic content introduced at the first-order to describe more completely the multiharmonic solutions that exist in this nonlinear media. As will be discussed in Sec. 5.3, this higher-order solution, when added to the first- and zeroth-order solutions, propagates more invariantly, i.e., with less growth and decay in its spectral content.

After substituting the results from Eqs. (18), (26), and (34) into Eq. (12), the solution for xj is now known up to O(ε3).

Method of Multiple Scales: Diatomic Chain.

Next, the MMS is developed for the diatomic chain. As with the monoatomic chain, higher-order results reveal the magnitudes of the spatially and temporally invariant harmonics induced by the quadratic and cubic nonlinearities as well as amplitude and phase evolution equations utilized in the stability and dispersion shift analyses. It is again assumed that there are multiple time scales at which the system evolves, and a series expansion of the solution is assumed 
x¯j=x¯j(0)+εx¯j(1)+ε2x¯j(2)++εnx¯j(n)
(45)
where x¯j(n) denotes the nth-order solution vector 
x¯j(n)=[xa(n)(j,t)xb(n)(j,t)]
(46)
Collecting matching orders of ε in Eq. (2) yields n+1 vector equations. The first two are presented here 
ε0:MD02x¯j(0)+l=j1j+1Klx¯l(0)=0¯
(47)
 
ε1:MD02x¯j(1)+l=j1j+1Klx¯l(1)=2MD0D1x¯j(0)+f¯damping,1+f¯quadratic,1+f¯cubic,1
(48)
 
f¯damping,1=[cD0(xb(0)(j,t)xa(0)(j,t))+cD0(xb(0)(j1 ,t)xa(0)(j,t))cD0(xa(0)(j+1,t)xb(0)(j,t))+cD0(xa(0)(j,t)xb(0)(j,t))]
(49)
 
f¯quadratic,1=[k2(xb(0)(j,t)xa(0)(j,t))2k2(xb(0)(j1,t)xa(0)(j,t))2k2(xa(0)(j+1,t)xb(0)(j,t))2k2(xa(0)(j,t)xb(0)(j,t))2]
(50)
 
f¯cubic,1=[k3(xb(0)(j,t)xa(0)(j,t))3+k3(xb(0)(j1,t)xa(0)(j,t))3k3(xa(0)(j+1,t)xb(0)(j,t))3+k3(xa(0)(j,t)xb(0)(j,t))3]
(51)
By inspection of Eq. (47), the zeroth-order solution takes the Bloch form 
x¯j(0)=12A ϕ¯eiω0T0eiμj+c.c.
(52)

where A denotes the complex amplitude and ϕ¯ a wave mode shape.

Substituting Eq. (52) into Eq. (47) leads to an eigenvalue problem for two wave modes and their dispersion relationship, ω0(μ). There are two independent solutions, corresponding to the acoustic and optical branches 
ϕac¯=[(1+eiμ)mbmbma+ma2+mb2+2mambcosμ1]
(53)
 
ϕopt¯=[1mamb+ma2+mb2+2mambcosμ(1+eiμ)mb]
(54)
 
ω0,ac=k1(mb+ma)mambk124mb2+4ma2+8mambcosμma2mb2
(55)
 
ω0,opt=k1(mb+ma)mamb+k124mb2+4ma2+8mambcosμma2mb2
(56)

As with the monoatomic chain, these zeroth-order dispersion relationships completely characterize the linear, undamped chain (ε=0). Higher-order results add amplitude-dependent corrections to these relationships due to the nonlinearities.

It is worth noting that at μ=π, when mb>ma, the optical mode shape simplifies from an indeterminate form to [10] through application of L'Hopital's rule. Also, the mode shapes presented here are the complex conjugates of those in Ref. [8] due to a difference in sign convention for the zeroth-order solution in Eq. (52).

The known form of x¯j(0) is then substituted into Eq. (48). Projecting the secular terms against the Hermitian conjugate of the corresponding mode shape leads to evolution expressions for D1α and D1β 
D1α=γdα
(57)
 
D1β=δdα2
(58)
where 
γd=γd(μ,c,ϕa¯,ma,mb)
(59)
 
δd=δd(μ, k3, ϕa¯, ϕb¯, ma,mb, ω0)
(60)

The functional dependence only is provided although complete expressions can be found in a straight-forward manner using a symbolic manipulator. Such dependence is significant because of arguments made about the influence of different chain parameters on dispersion shifts and invariance in Secs. 5.1 and 5.3, respectively.

A particular solution for x¯j(1) can be found using the form 
x¯j(1)=12B1¯e2iω0T0e2iμj+12C1¯e3iω0T0e3iμj+c.c.
(61)
where the vectors B1¯ and C1¯ are found using the method of undetermined coefficients 
B1¯=[b1,aA2b1,bA2]
(62)
 
C1¯=[c1,a A3c1,bA3]
(63)
where 
b1,a=b1,a(μ, ma, mb, k1, k2,ϕa¯, ϕb¯, ω0)
(64)
 
b1,b=b1,b(μ, ma, mb, k1, k2,ϕa¯, ϕb¯, ω0)
(65)
 
c1,a=c1,a(μ,ma, mb, k1,k3, ϕa¯,ϕb¯, ω0)
(66)
 
c1,b=c1,b(μ,ma, mb, k1,k3, ϕa¯, ϕb¯, ω0)
(67)

These amplitudes detail how much second and third harmonics arise due to the quadratic and cubic nonlinearities, respectively, and will be shown to make up an invariant waveform.

Note that this solution may not exist very close to μ values associated with internal resonance—i.e., where (2ω0, 2μ) or (3ω0, 3μ) satisfy the zeroth-order dispersion relationship in Eq. (55) or (56), inducing a resonant response. For these cases, the particular solution becomes unbounded and must be removed using an additional (e.g., detuning) parameter. Chains exhibiting this behavior are omitted in this work, but this phenomenon can be explored in greater detail in future studies.

Reapplying the procedure on the second-order yields evolution equations for α and β of the following form: 
D2α=D2α(α,μ,ma,mb,k2,k3,c,γd,δd,ϕa¯,b1,a,b1,b,c1,a,c1,b,ω0)
(68)
 
D2β=D2β(α,μ,ma,mb,k2,k3,c,γd,δd,ϕa¯,b1,a,b1,b,c1,a,c1,b,ω0)
(69)
The second-order particular solution then takes the form 
x¯j(2)=12B2¯e2iω0T0e2iμj+12C2¯e3iω0T0e3iμj+12E2¯e4iω0T0e4iμj+12 F2¯e5iω0T0e5iμj+c.c
(70)

where the vectors B2¯, C2¯, E2¯, and F2¯ can be found using the method of undetermined coefficients.

As with the monoatomic chain, second-order analysis adds to the spectral content at the second and third harmonics and introduces new spectral content at the fourth and fifth harmonics further defining the multiharmonic invariant wave. It is also noted here that the diatomic results can recover the monoatomic results by setting ma=mb=m and replacing μ with μ/2.

Numerical Simulation

This section briefly presents the procedure for numerically injecting perturbation-based waveforms into numerical simulations of Eqs. (1) and (2), forming the foundation to the numerical results presented in the remaining sections. They rely on initial conditions associated with the multiple scales solutions, assigning the order of terms at which to truncate the multiharmonic wave solution. In doing so, this procedure exhibits more control over the ensuing wave propagation than the approach used in Refs. [7] and [9], whereby the chain is excited with an arbitrary forcing until a measurable response is generated. It also allows for measurement of the validity of the MMS solutions since waveforms existing in the medium should persist after they are initially set, as detailed in Sec. 5.3. This section also introduces dimensionless parameters that combine the effects of nonlinear stiffness and wave amplitude into convenient expressions useful for simulating chains and interpreting results.

Displacement and velocity initial conditions required to numerically integrate Eqs. (1) and (2) are found using the series solutions in Eqs. (12) and (45), respectively. Thus, x(j,0)=xj|t=0 and ẋ(j,0)=(d/dt)(xj)|t=0, where xj|t=0 is truncated to a desired order n, and (d/dt)() must correspond to the total time derivative defined in Eq. (10), truncated at the same order n.

By carrying the expansions out and truncating, first- and second-order initial conditions can be determined. To emulate an infinite domain, the number of unit cells N should be sufficiently large—generally on the order of 1800–2000—such that μN1. Within this domain, analysis of a chain's behavior remains at a central region—usually the middle 800–1000 unit cells—such that reflection from the boundaries does not propagate to the region within the time scale of the simulation. In these studies, either fixed–fixed or free–free boundary conditions are used.

A given simulation is characterized by the following dimensionless parameters introduced here: 
τ=ω0t,Π1=εcω0k1, Π2=εk2α0k1, Π3=εk3α02k1, Π4=mbma  
(71)

where Π1, Π2,  and Π3 can be interpreted as relative strengths of the damping, quadratic nonlinearity, and cubic nonlinearity, respectively, and Π4, a mass ratio, applies only to the diatomic chain. The parameters Π2 and Π3 are particularly useful since they combine the effects of the nonlinear stiffness and wave amplitude into single parameters, both of which affect the magnitude of dispersion shifts. Note, however, that Π2 and Π3 are proportional to—not exact measurements of—the strength of the quadratic and cubic nonlinearities, respectively. That is to say, these parameters are not an exact ratio of spring forces because α0 is the zeroth-order wave amplitude and not necessarily the spring stretch, xj+1xj. Thus, Π2,Π3<0.1 are conservative bounds on the MMS applicability (weak nonlinearity assumption). In general, the parameters in Eq. (71) facilitate parameter definition as outlined in the following procedure outline.

These steps were carried out for each simulation:

  • Specify Π1, Π2, and Π3 values. Set m and k1 equal to 1, without loss of generality. For diatomic simulations, prescribe ma, mb values (based on a desired Π4 value) as well as a branch of excitation (acoustic or optical).

  • Specify a propagation constant μ restricted to the first irreducible Brillouin zone. A corresponding ω0 value can then be assigned to it according to Eqs. (20), (55), or (56). Note that this is not guaranteed to be the ensuing temporal frequency of the chain but merely a parameter used to evaluate initial conditions. For diatomic simulations, ϕac¯ or ϕopt¯ can now be evaluated, depending on the branch of excitation that was selected.

  • Set k3 equal to 1 (hardening) or −1 (softening) and then solve for c, k2,  and α0 values which satisfy the desired Π1,Π2, and Π3 values, respectively.

  • Define the number of unit cells to emulate an infinite system (Usually N=1800 or 2000 such that μN1)

  • Assign an initial displacement and velocity distribution to the chain according to Eqs. (12) or (45) evaluated at time t=0 and truncated at a specified order n (n=0,1, or 2). The ε parameter is chosen to be small (e.g., ε=0.1).

  • Set two boundary conditions. For each end, use either fixed—zero displacement—or free—zero force.

  • Simulate the equations of motion by numerically integrating Eq. (1) or (2) using matlab's ODE45 routine.

  • For analyzing the results, study only a fixed middle portion of the chain (e.g., the middle 800–1000 unit cells) such that boundary effects can be ignored. Visually inspect the results to make sure that the boundary effects have not propagated into this central region during the time range of interest.

Results

Dispersion Shifts.

This section details amplitude-dependent dispersion shifts arising from both cubic and quadratic nonlinearities, a finding useful for device design [7]. It will be shown that the MMS analysis uncovers dispersion shifts to the zeroth-order (ε=0, linear undamped) relationship. The measurement of these predicted shifts in numerical simulations affirms the validity of the MMS results at the first and second-orders and demonstrates that quadratic nonlinearities govern shifts at higher-orders only, i.e., an effect not captured by a first-order analysis.

At the zeroth-order, the temporal frequency ω0 relates to the spatial frequency μ through dispersion relationships represented by Eqs. (20), (55), and (56). Figure 4 displays these zeroth-order expressions for the monoatomic and diatomic chains. Note that the monoatomic chain has a cutoff for ω0>2(k1/m), and the diatomic chain exhibits a bandgap.

Closed-form, amplitude-dependent dispersion shifts will now be derived from the MMS results. For undamped chains, D1α=D2α=0. Thus, the wave amplitude is constant (α=α0). For the monoatomic chain, the evolution equations for phase, Eqs. (23) and (30), reduce to 
D1β=δα02,     c=0     
(72)
 
D2β=Xα04+Yα02+Z,    c=0     
(73)

Similar expressions result for the diatomic chain by replacing δ in Eq. (72) with δd from Eq. (60), and replacing the left side of Eq. (73) with the expression in Eq. (69) evaluated at α=α0.

Reconstituting the phase expression using D1β and D2β in Eq. (10) and returning to the original time scale, leads to the wave's temporal phase behavior 
β̇=εD1β+ε2D2β=εδα02+ε2(Xα04+Yα02+Z)
(74)

for the monoatomic chain; a similar expression occurs for the diatomic chain.

Since D1β and D2β do not explicitly depend on time, they can be interpreted as corrections to the frequency at their respective order, yielding first- and second-order-accurate corrections 
ωO(ε1)=ω0+εD1β
(75)
 
ωO(ε2)=ω0+εD1β+ε2D2β
(76)

Inspection of Eqs. (72), (73), (75), and (76) shows that the quadratic nonlinearity does not shift the frequency until the second-order correction is included, implying a small influence as compared to the cubic nonlinearity, which appears at both the first- and second-orders. For positive coefficients k2 and k3, the shifts further increase the frequency (for a given μ) with increasing wave amplitude α0. Such effects can inform the design of tunable devices (filters, waveguides, acoustic logic devices, etc.) which selectively include or exclude waves based on wave gain and subsequent bandgap shifting [7,8].

For both monoatomic and diatomic chains, numerical studies of undamped systems confirm the theoretical shifts as documented in Fig. 5. For different μ values prescribed in simulation initial conditions, the first- and second-order MMS-predicted temporal frequencies are compared to the measured dominant temporal frequency of a central unit cell in the chain. Recall that assigning a μ value does not guarantee the ensuing temporal frequency but rather is used to define all initial parameters. Thus, the simulations given first- and second-order initial conditions (n=1, 2) were measured to have dominant temporal frequencies close to the first- and second-order perturbation corrections (ωO(ε1), ωO(ε2)), respectively. Theoretical expressions for the zeroth-order, unshifted temporal frequency, ω0 (Eq. (20)—monoatomic and Eqs. (55) and (56)—diatomic) are also plotted as references to visualize the shifts that occur. For each simulation, the dominant temporal frequency is computed by taking an fast Fourier transform (FFT) of the time history of a central mass in the chain. To collect accurate data, temporal sampling frequencies are many integer multiple times higher than the expected shifted frequency, e.g., fsampling=600 ((ωO(ε1))/2π) or fsampling=600((ωO(ε2))/2π).

Inspection of the bottom two subplots in Fig. 5 reveals that k2 bears little to no effect on the first-order dispersion correction. As Π2 changes from 0.04 to 0, which corresponds to the no quadratic nonlinearity present, the first-order perturbation and simulation results change by an average amount of 0.043%. For lightly damped chains D1α, D2α0, corrections can be interpreted as instantaneous dispersion shifts since the amplitude decays slowly. The nonconstant reconstituted expressions for β̇ are then interpreted similar to Eqs. (75) and (76) at each instance of time.

Stability Considerations.

The stability of wave propagation is of primary concern for designing devices informed by these discrete systems. Herein, instability is defined as any qualitative shift from the multiharmonic plane wave solution predicted by the MMS. Thus, it is important to characterize the types of instabilities that take place in these systems and the conditions at which they occur. An analytical stability analysis is first presented to see if the evolution equations for the amplitude reveal any information on stability loss. Informed by this analysis, numerical studies are carried out to determine the onset of instabilities.

Local Stability Analysis.

Because it is a critical parameter that determines the strength of the nonlinearities in the system, amplitude is chosen to be analyzed to assess the stability of the wave. It will be shown that the first-order results reveal only a damping-induced decay of amplitude (stable solution), whereas the inclusion of second-order results reveals a possible growth in amplitude (unstable solution). Though it is outside the conservative bounds of the MMS solution, the unstable solution provides a foundation to the numerical stability studies presented in Sec. 5.2.2.

Applying Eq. (10) to wave amplitude α and recalling that D0α=0 leads to its reconstituted evolution equation 
α̇=εD1α+ε2D2α
(77)
which suggests the existence of fixed points, or α* such that α̇=0. These fixed points are found from 
α̇|α*=0
(78)
A local stability analysis can then be performed near each fixed point by introducing small perturbations from α*. The result of doing so is that the stability of each fixed point depends on the sign of λ arising from 
ddαα̇|α*λ
(79)
where λ>0 yields an unstable fixed point, λ<0 a stable fixed point, and λ=0 is neutrally stable. Thus, the fixed points for the monoatomic chain can be readily found by substituting in the values from Eqs. (22) and (29) into Eqs. (77) and (78) 
α*=0,±ω0εδ
(80)
The fixed point α*=±(ω0/εδ) exists only when δ>0 (cubic hardening) since α is assumed to be strictly real. Application of Eq. (79) yields 
λ|α*=0 =εγ
(81)
 
λ|α*=±ω0εδ=2εγ
(82)

The fixed point α*=0 is of the highest relevancy since (1) amplitudes near it are small-enough to satisfy the conservative weakly nonlinear ansatz, Π2,Π3<0.1; (2) this is the fixed point to which asymptotically predicted plane waves will be attracted to in the presence of damping, therefore, its stability establishes plane wave stability of interest herein. Since, λ|α*=0 <0 for damped chains (c>0), the fixed point α*=0 is stable. Note that α*=0 is the only fixed point informed by a first-order analysis as α*=±(ω0/εδ) originates from second-order terms. Also note that the presence of damping enables the stability of the fixed points to be classified: undamped chains yield λ|α*=0 =λ|α*=±(ω0/εδ)=0, which is neutrally stable.

By similar arguments, since Eq. (82) is positive for damped chains (c>0), the fixed point α*=±(ω0/εδ) is unstable. However, closer inspection of this fixed point reveals that it is at a large-enough amplitude that its accuracy cannot be expected. After substituting Eqs. (20) and (25) into the expression for α*, a critical Π3 value results dependent on only μ 
Π3,crit=εk3α*2k1=8(1cosμ)3(cos2μ4cosμ+3)
(83)

Examining Eq. (83), the lowest value of Π3,crit in the first Brillouin zone, 0μπ, is approximately 0.66, which is beyond the conservative weak nonlinearity ansatz used for MMS analysis. However, as evident in the results presented in Sec. 5.2.2, the existence of such a fixed point is likely since large-amplitude plane waves are observed to grow unboundedly in numerical simulations.

The stability of the monoatomic chain is summarized in Fig. 6. As indicated, only initial wave amplitude, and not initial phase, determines the basins of attraction.

Stability of the diatomic chain can also be studied analytically. Its reconstituted evolution equation for α can be written as 
α̇=P1α5+P2α3+P3α
(84)
where P1,P2,andP3 denote constants known from Eqs. (57) and (68). Fixed-points can be solved for by setting α̇=0. Thus, a given system admits at most five fixed points; however, because α is assumed to be real, imaginary fixed points are ignored. Inspection of Eq. (84) indicates that 
α*=0,±[P22P1±12P22P124P3P1]12
(85)
where it is again emphasized that α* must be strictly real. The stability of a given fixed point can be studied by examining the sign of each value of λ 
λ=5P1α*4+3P2α*2+P3
(86)

Results are similar to the monoatomic chain. For α*0, many values of λ can be computed using the P1, P2, and P3 values associated with a prescribed set of Π1, Π2,Π3, and Π4 values and branch. The fixed point α*=0 is always stable for c>0. Other fixed points are unstable but have corresponding Π2 and Π3 values never simultaneously less than 0.1. Therefore, they are of large-enough amplitude to not be of interest to plane waves predicted herein.

Numerical Stability Studies.

This section details numerical studies that investigate stability beyond the applicable range of the MMS. It also confirms the general trends predicted by the local stability analysis, which is that at low amplitudes, multiharmonic plane waves are stable, and at high amplitudes (or, equivalently, large nonlinearities) the plane wave solutions predicted may lose stability. Note that such a conclusion cannot be drawn from only first-order information. These studies quantify the stability of the predicted waveforms over a wide range of Π2 and Π3 values for a given value of Π1 and, if diatomic, Π4. Simulations begin with the initial conditions corresponding to the series solutions in Eqs. (12) and (45) for the monoatomic and diatomic chain, respectively. An unstable response is defined to be any deviation from the nature of the plane wave propagation predicted by the perturbation analysis. Thus, such studies aim to identify the stability limits of the MMS-predicted multiharmonic solutions. Although different types of instability are observed, it is not the primary concern of this work to explore further the ensuing solutions.

One such type of instability that occurred in numerical simulations is characterized by a spreading of the wave's spectral content. Chains with large Π3 values exhibit this behavior. The spectral content of these waves initially contains integer multiples of ω0 and μ, but following longer simulation times, the plane wave nature of the response breaks down into an assortment of incommensurate frequencies. Figure 7 depicts this phenomenon.

The other type of instability observed in simulations is characterized by an unbounded growth in the wave's amplitude. Chains with large positive or large negative Π2 values, as well as chains with large negative Π3 values, exhibit this behavior. The initial amplitude is unstable such that, as time progresses, the wave propagates with increasingly larger amplitudes. Note that this is similar to the exponential growth in amplitude instability presented in Ref. [28]. Figure 8 depicts this phenomenon. It is noted that all oscillations cease when the amplitude starts to grow.

For automating instability detection in the simulated results, criteria are developed which determine when one of the two types occur, thus signifying a loss of stability, or deviation from the multiharmonic plane wave solution. When either one of the two instability types occur, the parameter set is considered “unstable,” without distinguishing exactly which type of instability developed. Thus, an instability occurs when a set number of peaks in a spatial FFT exceeds a certain fraction of the maximum peak at that instant or when the instantaneous amplitude grows to a certain multiple of the simulation's original amplitude. A time limit is prescribed for simulations: if an instability does not arise within the prescribed time, or τ=12(2π) for these simulations, the parameter set is considered stable. For the diatomic chain, instability criteria are applied to xa and xb separately.

Figure 9 presents the stability results for the monoatomic chain as a function of Π2 and Π3. The parameters Π2,crit and Π3,crit are the values of Π2 and Π3, accurate to the nearest tenth, at which an unstable response is first observed when the other nonlinear Π-value is set to zero (i.e., Π2,crit occurs when Π3=0 and Π3,crit occurs when Π2=0). Each plot verifies the local stability analysis; namely, that plane waves characterized by small values of Π2/Π2,crit and Π3/Π3,crit are stable. Figure 10 presents similar stability results for the diatomic chain applied to xa. Plots for the stability study of xb depict similar trends and are shown in Figs. 1315 in the Appendix.

For some chains simulated using free–free boundary conditions, a rigid body (DC) harmonic appears that causes the two masses in each unit cell to separate by a constant amount, in addition to their expected oscillation. The location of the center of mass remains unchanged, as expected to satisfy the conservation of momentum. However, since the magnitude of this separation is generally small for the parameters being studied, such responses are considered stable.

Although the amplitudes at which point stability is lost in the numerical simulations do not align with the nonzero fixed points predicted by multiple scales approach, the results of these numerical studies qualitatively confirm the findings. Wave amplitudes below a critical value decay over time toward zero due to the presence of damping. Above this critical amplitude value, the wave will become unstable by either growing unboundedly or decomposing its original multiharmonic content into incommensurate frequencies.

Invariance.

This section describes how higher-order MMS results converge to an invariant waveform, a finding that offers new insight for predicting and modeling the wave propagation in these weakly nonlinear discrete systems and further validation of the MMS approach. In the absence of damping, the presented perturbation approach does not distinguish between any two points in space, nor any two points in time. Thus, in the limit as the asymptotic expansions are taken to higher orders, and assuming convergent vis-à-vis divergent behavior, it can be hypothesized that the waveforms should propagate for all space and time without change—i.e., they should be an invariant form akin to cnoidal waves studied in other contexts. It is worth noting that such an argument does not rely on a continuum approximation as in other studies [2123]. Note that this is a very different perspective taken from other works carrying-out perturbation approaches in wave propagation problems [14,35], which more conventionally analyze higher harmonic generation and/or the closely related concept of phonon generation and lifetimes.

To test the invariant hypothesis, numerical integration of Eqs. (1) and (2) is carried out using the procedure presented in Sec. 4. Thus, plane waves corresponding to the asymptotic expansions in Eqs. (12) and (45) truncated at different orders are injected into the simulated chains. During simulation, the spatial harmonic content of the waves, is quantified using FFTs at different instances in time to track the generation (or loss) of the initially prescribed harmonic content. As mentioned in Sec. 4, to emulate an infinite chain, the analysis was restricted to a central region of the chain so that boundary effects could be ignored. According to the hypothesis, waves incorporating higher-order terms should exhibit measurably less variation, and hence propagate more invariantly. Only Π3>0 chains were sampled since it was found while conducting the numerical stability studies in Sec. 5.2.2, that the negative Π3 values become unstable at much lower magnitudes than positive ones.

For monoatomic chains, a parameter that measures the spatial variance of a simulated waveform is proposed as 
σp(n)1Ll=1L|1al,p|
(87)

where al,p denotes the amplitude of the pth harmonic at the lth instant in time normalized by its initially prescribed value, L is the number of time instances that make up the numerical solution, and n denotes the order of perturbation terms retained in Eq. (12). Since σp(n) measures variance, lower values of σp(n) indicate more invariant waveforms, and thus, it is expected (if the hypothesis holds) that σp(n) decreases as n increases. Figure 11 verifies the expected reduction in variance for monoatomic chains given second versus first-order initial conditions—all values of (Π2, Π3) simulated experience decreases in variance. Because Fig. 11 plots the normalized reduction in variance— σp(1)σp(2)/σp,max(1) (p=2,3) —a value of 1.0 indicates elimination of all variance using second-order initial conditions, while a value of 0.0 indicates no variance reduction.

For diatomic chains, the variance parameter specializes to each degree-of-freedom: σp,a(n) and σp,b(n) for xa and xb, respectively. Similar numerical simulations are carried out and the second and third spatial harmonics are tracked over time for waves on the acoustic and optical branches. Figure 12 presents the variance results for xa; the Appendix provides the results for xb. The trend in variance reduction is again present.

Figures 11 and 12 also show that an increase in the magnitude of Π2 leads to a large variance reduction in the third harmonic, and conversely, an increase in Π3 leads to a large variance reduction of the second harmonic. This counter-intuitive crossing effect can be understood by close examination of the second-order results in Eqs. (39)(42) for the monoatomic chain (the same argument can be made for the diatomic chain). There, k2 influences the particular solution governing the third harmonic, and k3 influences the particular solution governing the second harmonic. Note that, the first-order solution does not contain such a crossing effect and thus injected waveforms based solely on it will not be sufficient for preventing an exchange of energy between different harmonics. In contrast, because the presence of the crossing terms in second-order injected waves keeps spectral energy, where it is initially assigned, the crossing effect becomes necessary for increasing invariance. This is yet another example of the richness in system behavior uncovered by higher-order perturbation analysis.

Conclusions

Higher-order multiple scales solutions have been detailed for predicting multiharmonic plane wave propagation in cubically and quadratically nonlinear monoatomic and diatomic chains. The method yields higher-order amplitude-dependent dispersion relationships, the result of which clearly show the importance of cubic terms vis-à-vis quadratic terms in frequency/wavenumber shifts. Wave stability has also been assessed through introduction of linear damping and a local stability analysis, and the notion of wave invariance has been explored and demonstrated for waveforms composed using higher-order terms in the asymptotic expansions. Finally, trends in dispersion shifts, wave stability, and wave invariance have been verified using dimensionless parameters and direct numerical integration.

Acknowledgment

The authors would like to thank the National Science Foundation for support of this research under Grant No. CMMI 1332862.

Appendix

References

References
1.
Li
,
F.
,
Anzel
,
P.
,
Yang
,
J.
,
Kevrekidis
,
P. G.
, and
Daraio
,
C.
,
2014
, “
Granular Acoustic Switches and Logic Elements
,”
Nat. Commun.
,
5
, pp.
1
6
.
2.
Jensen
,
J. S.
,
2003
, “
Phononic Band Gaps and Vibrations in One- and Two-Dimensional Mass–Spring Structures
,”
J. Sound Vib.
,
266
(
5
), pp.
1053
1078
.
3.
Sun
,
H.-X.
,
Zhang
,
S.-Y.
, and
Shui
,
X.-J.
,
2012
, “
A Tunable Acoustic Diode Made by a Metal Plate With Periodical Structure
,”
Appl. Phys. Lett.
,
100
(
10
), p.
103507
.
4.
Vakakis
,
A. F.
, and
King
,
M. E.
,
1995
, “
Nonlinear Wave Transmission in a Monocoupled Elastic Periodic System
,”
J. Acoust. Soc. Am.
,
98
(
3
), pp.
1534
1546
.
5.
Chakraborty
,
G.
, and
Mallik
,
A. K.
,
2001
, “
Dynamics of a Weakly Non-Linear Periodic Chain
,”
Int. J. Nonlinear Mech.
,
36
(
2
), pp.
375
389
.
6.
Dubus
,
B.
,
Swinteck
,
N.
,
Muralidharan
,
K.
,
Vasseur
,
J. O.
, and
Deymier
,
P. A.
,
2016
, “
Nonlinear Phonon Modes in Second-Order Anharmonic Coupled Monoatomic Chains
,”
ASME J. Vib. Acoust.
,
138
(
4
), p.
041016
.
7.
Narisetti
,
R. K.
,
Leamy
,
M. J.
, and
Ruzzene
,
M.
,
2010
, “
A Perturbation Approach for Predicting Wave Propagation in One-Dimensional Nonlinear Periodic Structures
,”
ASME J. Vib. Acoust.
,
132
(
3
), p.
031001
.
8.
Manktelow
,
K. L.
,
Leamy
,
M. J.
, and
Ruzzene
,
M.
,
2014
, “
Weakly Nonlinear Wave Interactions in Multi-Degree of Freedom Periodic Structures
,”
Wave Motion
,
51
(
6
), pp.
886
904
.
9.
Manktelow
,
K.
,
Leamy
,
M. J.
, and
Ruzzene
,
M.
,
2010
, “
Multiple Scales Analysis of Wave-Wave Interactions in a Cubically Nonlinear Monoatomic Chain
,”
Nonlinear Dyn.
,
63
(
1
), pp.
193
203
.
10.
Narisetti
,
R. K.
,
2010
, “
Wave Propagation in Nonlinear Periodic Structures
,”
Ph.D. thesis
, Georgia Institute of Technology, Atlanta, GA.
11.
Hussein
,
M. I.
,
Leamy
,
M. J.
, and
Ruzzene
,
M.
,
2014
, “
Dynamics of Phononic Materials and Structures: Historical Origins, Recent Progress, and Future Outlook
,”
ASME Appl. Mech. Rev.
,
66
(
4
), p.
040802
.
12.
Cabaret
,
J.
,
Tournat
,
V.
, and
Béquin
,
P.
,
2012
, “
Amplitude-Dependent Phononic Processes in a Diatomic Granular Chain in the Weakly Nonlinear Regime
,”
Phys. Rev. E
,
86
(
4
), p.
041305
.
13.
Sánchez-Morcillo
,
V. J.
,
Pérez-Arjona
,
I.
,
Romero-García
,
V.
,
Tournat
,
V.
, and
Gusev
,
V. E.
,
2013
, “
Second-Harmonic Generation for Dispersive Elastic Waves in a Discrete Granular Chain
,”
Phys. Rev. E
,
88
(
4
), p.
043203
.
14.
Swinteck
,
N. Z.
,
Muralidharan
,
K.
, and
Deymier
,
P. A.
,
2013
, “
Phonon Scattering in One-Dimensional Anharmonic Crystals and Superlattices: Analytical and Numerical Study
,”
ASME J. Vib. Acoust.
,
135
(
4
), p.
041016
.
15.
Hasegawa
,
A.
, and
Tappert
,
F.
,
1973
, “
Transmission of Stationary Nonlinear Optical Pulses in Dispersive Dielectric Fibers—I: Anomalous Dispersion
,”
Appl. Phys. Lett.
,
23
(
3
), pp.
142
144
.
16.
Laitone
,
E.
,
1960
, “
The Second Approximation to Cnoidal and Solitary Waves
,”
J. Fluid Mech.
,
9
(
03
), pp.
430
444
.
17.
Friedrichs
,
K. O.
,
1948
, “
On the Derivation of the Shallow Water Theory
,”
Commun. Pure Appl. Math.
,
1
, pp.
81
85
.
18.
Grosvenor
,
G.
,
1933
, “
Flying
,”
Natl. Geogr. Mag.
,
63
(
5
), pp.
585
630
.
19.
Fenton
,
J. D.
,
1999
, “
The Cnoidal Theory of Water Waves
,”
Developments in Offshore Engineering
,
J. B. Herbich and K. A. Ansari K.
,
eds.
, Gulf Publisher Company, Houston, TX, pp. 2–32.
20.
Kraenkel
,
R. A.
,
Manna
,
M. A.
,
Merle
,
V.
,
Montero
,
J. C.
, and
Pereira
,
J. G.
,
1996
, “
Multiple-Time Higher-Order Perturbation Analysis of the Regularized Long-Wavelength Equation
,”
Phys. Rev. E.
,
54
(3), pp.
2976
2981
.
21.
Friesecke
,
G.
, and
Mikikits-Leitner
,
A.
,
2015
, “
Cnoidal Waves on Fermi-Pasta-Ulam Lattices
,”
J. Dyn. Differ. Equations
,
27
(
3
), pp.
627
652
.
22.
Gaison
,
J.
,
Moskow
,
S.
,
Wright
,
J. D.
, and
Zhang
,
Q.
,
2014
, “
Approximation of Polyatomic FPU Lattices by KdV Equations
,”
Multiscale Model. Simul.
,
12
(
3
), pp.
953
995
.
23.
Friesecke
,
G.
, and
Pego
,
R. L.
,
1999
, “
Solitary Waves on FPU Lattices: I. Qualitative Properties, Renormalization and Continuum Limit
,”
Nonlinearity
,
12
(
6
), pp.
1601
1627
.
24.
Sen
,
S.
,
Hong
,
J.
,
Bang
,
J.
,
Avalos
,
E.
, and
Doney
,
R.
,
2008
, “
Solitary Waves in the Granular Chain
,”
Phys. Rep.
,
462
(
2
), pp.
21
66
.
25.
Kartashov
,
Y. V.
,
Malomed
,
B. A.
, and
Torner
,
L.
,
2011
, “
Solitons in Nonlinear Lattices
,”
Rev. Mod. Phys.
,
83
(
1
), pp.
247
305
.
26.
Molerón
,
M.
,
Leonard
,
A.
, and
Daraio
,
C.
,
2014
, “
Solitary Waves in a Chain of Repelling Magnets
,”
J. Appl. Phys.
,
115
(
18
), p.
184901
.
27.
Daraio
,
C.
,
Nesterenko
,
V. F.
,
Herbold
,
E. B.
, and
Jin
,
S.
,
2005
, “
Strongly Nonlinear Waves in a Chain of Teflon Beads
,”
Phys. Rev. E
,
72
(
1
), p.
016603
.
28.
Newton
,
P. K.
, and
Keller
,
J. B.
,
1987
, “
Stability of Periodic Plane Waves
,”
SIAM J. Appl. Math.
,
47
(
5
), pp.
959
964
.
29.
Bickham
,
S. R.
,
Kiselev
,
S. A.
, and
Sievers
,
A. J.
,
1993
, “
Stationary and Moving Intrinsic Localized Modes in One-Dimensional Monatomic Lattices With Cubic and Quartic Anharmonicity
,”
Phys. Rev. B
,
47
(
21
), pp.
14206
14211
.
30.
Friesecke
,
G.
, and
Pego
,
R. L.
,
2004
, “
Solitary Waves on Fermi–Pasta–Ulam Lattices: III. Howland-Type Floquet Theory
,”
Nonlinearity
,
17
(
1
), pp.
207
227
.
31.
Flach
,
S.
, and
Gorbach
,
A. V.
,
2008
, “
Discrete Breathers—Advances in Theory and Applications
,”
Phys. Rep.
,
467
(1–3), pp.
1
116
.
32.
Huang
,
G.
, and
Hu
,
B.
,
1998
, “
Asymmetric Gap Soliton Modes in Diatomic Lattices With Cubic and Quartic Nonlinearity
,”
Phys. Rev. B
,
57
(
10
), pp.
5746
5757
.
33.
Gorbach
,
A. V.
, and
Johansson
,
M.
,
2003
, “
Discrete Gap Breathers in a Diatomic Klein-Gordon Chain: Stability and Mobility
,”
Phys. Rev. E
,
67
(
6
), p.
066608
.
34.
Nayfeh
,
A. H.
, and
Mook
,
D. T.
,
1979
,
Nonlinear Oscillations
,
Wiley
, New York.
35.
Hamilton
,
M. F.
, and
Blackstock
,
D. T.
,
1998
,
Nonlinear Acoustics
,
Academic Press
,
San Diego, CA
.