Nanoptera in weakly nonlinear woodpile and diatomic granular chains
NNanoptera in Weakly Nonlinear Woodpile and Diatomic GranularChains
G. Deng , C. J. Lustri ∗ , and Mason A. Porter Department of Mathematics and Statistics, 12 Wally’s Walk, Macquarie University, NewSouth Wales 2109, Australia Department of Mathematics, 520 Portola Plaza, University of California, Los Angeles,California 90024, USA
Abstract
We study “nanoptera”, which are non-localized solitary waves with exponentially small but non-decayingoscillations, in two singularly-perturbed Hertzian chains with precompression. These two systems are woodpilechains (which we model as systems of Hertzian particles and springs) and diatomic Hertzian chains withalternating masses. We demonstrate that nanoptera arise from Stokes phenomena and appear as special curves,called Stokes curves, are crossed in the complex plane. We use techniques from exponential asymptotics toobtain approximations of the oscillation amplitudes. Our analysis demonstrates that traveling waves in asingularly perturbed woodpile chain have a single Stokes curve, across which oscillations appear. Comparingthese asymptotic predictions with numerical simulations reveals that this accurately describes the non-decayingoscillatory behavior in a woodpile chain. We perform a similar analysis of a diatomic Hertzian chain, and wethat the nanpteron solution has two distinct exponentially small oscillatory contributions. We demonstratethat there exists a set of mass ratios for which these two contributions cancel to produce localized solitarywaves. This result builds on prior experimental and numerical observations that there exist mass ratios thatsupport localized solitary waves in diatomic Hertzian chains without precompression. Comparing asymptoticand numerical results in a diatomic Hertzian chain with precompression reveals that our exponential asymptoticapproach accurately predicts the oscillation amplitude for a wide range of system parameters, but it fails toidentify several values of the mass ratio that correspond to localized solitary-wave solutions.
The behavior of a chain of frictionless spherical particles under compression is often characterized in terms ofthe interaction effects between particles. In particular, the Hertzian interaction [35] describes the relationshipbetween the repelling force from the contact area of adjacent frictionless spherical particles and the compressionbetween them.A chain of spheres with Hertzian interactions is called a “Hertzian chain” and has been the subject ofnumerous investigations in the last few decades [18, 66, 77]. Hertzian chains are experimentally realizable andhighly configurable, so there has been much research both on their fundamental physical properties and ontheir potential uses for practical engineering applications [73]. ∗ Corresponding Author. Electronic address: [email protected] a r X i v : . [ n li n . PS ] F e b R − δ F . . . . . . F
Figure 1: Schematic illustration of a monoatomic Hertzian chain with precompression. The chain consists of N identical, alignedspheres with radius R . The interaction potential between adjacent spheres is Hertzian and is governed by (2) with α = 3 / F to the chain compresses each sphere by a uniform distance δ . A general equation that governs particle interactions in a chain is the following system of differential–difference equations: m ( n )¨ x ( n, t ) = φ (cid:48) ( x ( n + 1 , t ) − x ( n, t )) − φ (cid:48) ( x ( n, t ) − x ( n − , t )) , (1)where n ∈ Z , the quantity m ( n ) is the mass of the n -th particle, x ( n, t ) is the position of the n -th particleat time t , a dot denotes differentiation with respect to time, a prime denotes differentiation with respect tospace, and the interaction potential between adjacent particles is φ ( r ) = c ( δ − r ) α +1 , r ≤ δ , r > , c = constant , (2)where α > δ is the equilibrium overlap of adjacent particles due to the precompression that is inducedby an external force. We obtain the Hertzian interaction potential with the choice α = 3 /
2. In Figure 1, weillustrate a Hertzian chain with precompression. The Hertzian potential (2) is 0 when the precompression is0, and it cannot take negative values. For algebraic convenience, we set c = 1 / ( α + 1) in all of our examples.A particle chain supports the propagation of non-dispersive solitary waves [30] if the exponent α + 1 inthe interaction potential (2) is larger than 2 (i.e., if α > | x ( n − , t ) − x ( n, t ) | /δ (cid:28)
1, one can expand the right-hand side of (1) asa Taylor series in the normalized deformation size. Neglecting terms of order O (( | x ( n − , t ) − x ( n, t ) | /δ ) )produces the frequently studied Fermi–Pasta–Ulam–Tsingou (FPUT) chain. In the long-wavelength limit, inwhich the characteristic spatial size of a solution is much larger than the lattice spacing, one can approximatethe discrete lattice equation by a continuous system. In the limit of both small deformations and long wave-length in comparison to particle size, the monoatomic Hertzian system with precompression reduces to theKortweg–de Vries (KdV) equation [66]. In this long-wavelength limit, one can approximate the solitary-wavebehavior in a Hertzian chain as a soliton solution of the KdV equation [50]. Building on this idea, the interac-tion of solitary waves in a Hertzian chain in the long-wavelength limit was described in [80] using a two-solitonsolution of the KdV equation.A remarkable feature of Hertzian chains (and higher-dimensional generalizations of them) is their tunability.It is possible to produce a wide range of dynamical behaviors through simple modifications of a homogeneousHertzian system, such as by incorporating different types of particle heterogeneities. Such modifications canproduce behavior that is very different to that of homogeneous Hertzian chains. For example, researchershave explored the behavior of Hertzian chains with impurities [25, 63, 79], disordered and random particlearrangements [33,48,58,59,61], quasiperiodicity [62], and particles that are composed of segments with differentmaterials (i.e., so-called “compound chains”) [20,67,85,86]. It has been demonstrated that impurities can cause .....Identical rods (a) A physical woodpile system . . . . . . : Particle of mass m : Particle of mass m : Spring with: spring constant k (b) A woodpile-chain model . . . . . . : Particle of mass m : Particle of mass m (c) A diatomic chain of particlesFigure 2: The chains of particles that we examine in the present study. The schematic illustration in panel (a) shows the physicalconfiguration of orthogonal rods that is known as a “woodpile chain”. The schematic in panel (b) shows an idealizedmathematical model of the physical configuration in (a). This model consists of heavy spherical particles in physicalcontact via Hertzian interactions and light spherical particles (which are sometimes called “resonators”) that are attachedto each heavy particle by a spring. The schematic illustration in panel (c) shows a diatomic chain of alternating heavy andlight spheres in physical contact via Hertzian interactions. We analyze the models in (b) and (c). solitary waves to scatter [25, 63, 79], and the transport and localization of energy in Hertzian systems withdisorder and quasiperiodicity were studied in [48, 61, 62]. In Hertzian chains with random masses, solitarywaves can delocalize, such that the energy that is carried by a solitary wave spreads among many particlesin a chain [33, 58, 59]. The scattering of solitary waves at an interface of a compound Hertzian system wasstudied in [67, 85]. Other phenomena that have been studied in numerical and experimental investigations ofheterogeneous Hertzian chains include energy trapping [20, 86], shock disintegration [20], and the generationof secondary solitary waves [85].It is also interesting to consider other types of heterogeneous chains, including (1) particle–rod systemsthat model so-called “woodpile chains” [18, 44, 47, 49, 53, 54, 88, 89] and (2) diatomic chains [11, 18, 34, 39, 42,43,46,64,70–72,74,83]. Both the woodpile and diatomic chains are feasible to obtain experimentally. One canuse the electromagnetic counterpart of the woodpile chain, which is known as a “woodpile photonic crystal”,to controllably manipulate electromagnetic waves [28, 52].A woodpile chain consists of orthogonally stacked slender rigid cylinders with some mass m [see Fig-ure 2(a)]. The interaction force along the direction of the stack is described by a Hertzian force, and onecan model the elastic deformation along the direction perpendicular to the stack direction using internal “res-onators”, where one determines the mass and coupling constant of these resonators based on the material andshape of the cylinders. When considering a stack of identical cylinders, each resonator has the same mass m and elastic constant k . We model such a woodpile chain as a monoatomic Hertzian chain in which each particleis connected to an external particle by a rod, which we treat as a linear spring. Each particle in the chain hasmass m , each external particle has mass m , and each linear spring has spring constant k . We illustrate thismodel woodpile chain in Figure 2(b). For the remainder of this study, we use the term “woodpile chain” torefer to this idealized system.Much of the existing work on woodpile chains has concentrated on their linear elastic responses [44,49,88].Some recent work studied the behavior of traveling waves in woodpile chains with no precompression [47, 89],which is a strongly nonlinear regime. Other studies have investigated the existence of discrete breathers inthe woodpile chain both without [54] and with precompression [53].A diatomic chain consists of particles with some characteristic, such as particle mass, that alternates inadjacent particles. We consider chains [see Figure 2(c)] in which even particles have mass m and odd particleshave mass m . Diatomic Hertzian chains have rich dynamics [18]. They can have both discrete breathersolutions [11, 83] and stationary shock waves [64], and a variety of other phenomena have also been studiedin them. This includes the propagation and scattering of nonlinear waves in both ordered and disordereddiatomic Hertzian chains [34, 70–72] and chaotic dynamics in diatomic Hertzian chains that are both dampedand driven [39]. = ct (a) A solitary wave x = ct (b) A one-sided nanopteron x = ct (c) A two-sided nanopteronFigure 3: Comparison of the profiles of (a) a standard solitary wave, (b) a one-sided nanopteron, and (c) a two-sided nanopteron thateach propagate at speed c . The solitary wave is localized spatially, whereas the nanoptera have non-decaying oscillatorytails on (b) one side or (c) both sides of the wave front. The waves in panels (a) and (c) propagate without decaying,but the wave in panel (b) cannot propagate indefinitely. In (b), the one-sided radiation draws energy from the wave front,leading to eventual decay of the wave front. Because this decay occurs over a long time scale, the one-sided nanopteron in(b) is said to be “metastable”. It has also been shown that diatomic Hertzian chain models can, in principle, support solitary wavesthat propagate indefinitely without decaying [42], and several prior studies [42, 46, 74] have examined sets ofparameters for Hertzian chains that produce non-decaying solitary wave solutions. These sets of parametersarise from satisfying a condition known as an “anti-resonance condition”. Subsequently, [60] examined non-decaying solitary waves in two-dimensional Hertzian systems. Hertzian chains can also have related wavesolutions that travel with particularly strong dispersion and attenuation. Such behavior occurs if a chainssatisfies a ca “resonance condition”, which one obtaines in a similar fashion to an anti-resonance condition [43,46, 74].
Typically, traveling-wave solutions in woodpile and diatomic Hertzian systems are not truly localized [42, 47,89]; instead, they also include a small non-localized traveling and thus take the form of a nanopteron [12].A nanopteron is the superposition of a central solitary wave and a persistent oscillation in the “far field”(which refers to the region in which the distance from the central solitary wave tends to infinity) on oneor both sides of the central wave . These oscillations typically have an exponentially small amplitude withrespect to some small parameter. We show examples of nanoptera in Figure 3. A genuine solitary wave isexponentially localized, but a nanopteron is localized only up to algebraic orders in the small parameter.References [43, 46, 74, 89] demonstrated in uncompressed woodpile chains and diatomic Hertzian chains thatnon-decaying oscillations are absent for a set of parameters that satisfy anti-resonance conditions. When suchanti-resonance conditions are satisfied, both types of chains admit genuinely localized solitary-wave solutionsthat propagate indefinitely.Nanoptera have been studied in many other particle chains, such as diatomic Toda chains [56, 68, 82, 84],diatomic FPUT chains [26, 38, 55, 84], and chains with an on-site nonlinear potential in which each particleis linearly coupled to its neighbors [41]. The existence of nanopteron solutions in a diatomic FPUT chainwas proven rigorously in [26, 38], and nanoptera in periodic Toda chains were studied using a multiple-scaleapproach in [84]. Recently, [27] studied diatomic FPUT chains using numerical continuation. They exploredhow nanopteron solutions for chains in which one mass is much larger than the other (i.e., the so-called “small-mass-ratio” regime) are related to traveling-wave solutions, including ones in regions of parameter space inwhich the mass ratio is not small.Anti-resonance conditions, at which the small oscillations disappear entirely and a nanopteron becomes agenuine solitary wave, have been identified in diatomic Toda chains [56, 82, 84], a diatomic FPUT chain [56],and a discrete nonlinear Schr¨odinger equation [1]. It was demonstrated in [82] that diatomic Toda chains have Several studies use the term “nanopteron” to refer only to the additional small oscillations. We do not follow this convention;instead, we use the term “nanopteron” to refer to a solution that includes both a central wave and such small oscillations. mass ratio at which the small oscillatory tail is absent. In [84], it was shown that that there exists a discreteset of mass ratios that cause the oscillations to vanish. Fixed-point methods were used in Xu et al. [89] toderive numerical schemes for studying nanopteron solutions for the Hertzian chain without precompression;in the process of this derivation, they also obtained anti-resonance conditions.More recently, an exponential asymptotic method was applied in [55,56] to study anti-resonance conditionsin diatomic Toda and FPUT chains. These studies derived asymptotic anti-resonance conditions for massratios for these systems. They also showed that the exponentially small oscillations in these diatomic systemsare examples of the “Stokes phenomenon”, which refers to behavior that is switched on when special curves(called “Stokes curves”) are crossed in the complex plane. Traveling-wave solutions in diatomic Toda andFPUT chains possess two Stokes curves, which generate two distinct oscillations with the same amplitudebut different phases. For particular values of the mass ratio, the oscillations are precisely out of phase. Theoscillations therefore cancel, producing localized solitary waves. In the present study, we use exponentialasymptotic methods to study nanoptera in singularly perturbed woodpile and diatomic Hertzian chains.As we illustrated in Fig. 3, nanoptera can have oscillations on one or both sides of a central traveling wave.The existence of non-decaying symmetric two-sided nanopteron solutions, where the amplitudes of the oscilla-tions on both sides is identical, has been proved for the fifth-order KdV equation [40]. Subsequently, symmetricnanopteron solutions for the fifth-order KdV equation were constructed numerically using perturbation series.However, attempts to numerically compute one-sided nanoptera to all orders were unsuccessful [13]. It wasthen proved in [6] that one-sided nanopteron solutions for the fifth-order KdV equation cannot propagate in-definitely without changing form. The one-sided oscillation draws energy from the central wave, which causesthe amplitude of the central core to decay. The eventual decay of nanopteron solutions in the fifth-order KdVequation is a result of energy conservation, and one can make a similar argument against the existence ofone-sided nanoptera in any energy-conserving system. The instability of the one-sided nanopteron solutionsin the fifth-order KdV equation was discussed further in [32].Giardetti et al. [31] simulated one-sided nanoptera in the diatomic FPUT model with small mass ratios fora long enough time to obtain a solution that appeared to be a stable traveling wave. In their simulations, boththe central solitary wave and the oscillatory tail traveled without any apparent decay. However, for diatomicchains without driving or damping forces, the total energy of the system must be conserved. Consequently, itis impossible for both the central core and the one-sided oscillation to have constant amplitude indefinitely,because the one-sided oscillation constantly draws energy from the central core [31, 42, 84]. Therefore, it wasconjectured in [31] that the decay of one-sided nanoptera in diatomic FPUT chains occurs on a time scale that isexponentially large in the limit that the mass ratio parameter tends to 0. Solutions that decay on exponentiallyslow time scales are sometimes known as “metastable” solutions (which are also called “quasistable” solutions),which describe asymptotic solutions that appear to be stable over a duration that is large in comparison to anyinverse power of the small asymptotic parameter before eventually decaying. The nanoptera that we considerare all metastable; however, as we will note at an appropriate point in our discussion, it is straightforward toextend our analysis to the stable case of two-sided nanoptera.In the present study, we compute the behavior of nanopteron solutions in precompressed woodpile anddiatomic Hertzian chains in the asymptotic limit that m /m →
0. In this limit, both woodpile and diatomicHertzian chains are singularly perturbed around a monoatomic Hertzian chain. The study in [42] investigatedthis limit for diatomic Hertzian chains, but existing studies on woodpile chains have not considered the small-mass-ratio regime.We apply an exponential asymptotic method, which was developed in [17, 69] and was used in [55, 56] tostudy traveling waves in diatomic Toda and FPUT chains with small mass ratios. In a typical asymptoticpower-series analysis, one expands a solution as a series in some small parameter. One then computes theterms in the series using a recursion relation that one obtains by asymptotic matching. This process can nevercapture exponentially small behavior, which is smaller than any algebraic series term in the asymptotic limit.By contrast, the exponential asymptotic method that we use in the present paper is capable of describingasymptotic behavior on this exponentially small scale. Additionally, it only requires the direct computationof the leading-order series behavior. We describe this exponential asymptotic approach in detail in Section 2. .3 Paper Outline We first investigate traveling-wave solutions in a woodpile chain with precompression in the limit that thewavelength of the central traveling wave is large in comparison to the particle size; this is a weakly nonlinearregime. We show that the traveling waves in this system are nanoptera and that the appearance of non-decaying exponentially small oscillations behind the wave front is an example of the Stokes phenomenon.Using an exponential asymptotic approach, we obtain simple asymptotic expressions for the exponentiallysmall, constant-amplitude waves in the wake of the leading-order solitary wave. We demonstrate that thesolution has a single oscillatory contribution in the wake of the leading wave and thus that it can never becanceled out. Consequently, the woodpile chains that we examine do not have an anti-resonance condition.This contrasts with the results of Xu et al. [89], who identified anti-resonance conditions for a woodpile systemwith no precompression.We then study traveling waves in a diatomic Hertzian chain using the same exponential asymptotic ap-proach. The key difference between Hertzian chains and woodpile chains is that traveling waves in diatomicHertzian chains have two Stokes curves, which yields solutions with two oscillatory contributions. For cer-tain mass ratios, the two oscillatory contributions cancel precisely and the traveling-wave solution becomes agenuine localized solitary wave.For both woodpile chains and diatomic Hertzian chains, we use the long-wavelength approximation from [66]to describe the leading-order solitary wave behavior (as in [55]). This differs from a diatomic Toda system [56],for which there exists an analytical expression for the leading-order wave behavior. Therefore, our analysisof woodpile and diatomic Hertzian chains have two distinct small parameters, which correspond to the smallmass ratio and the large traveling-wave length scale.Our paper proceeds as follows. In Section 2, we discuss the exponential asymptotic method that we employfor our analysis. In Sections 3.1–3.3, we use this exponential asymptotic method to obtain an asymptoticapproximation for nanopteron solutions in a singularly perturbed woodpile chain. In Section 3.4, we comparethe results of this approximation with numerical computations. We obtain strong agreement between theasymptotic and numerical results. In Sections 4.1–4.3, we perform a similar exponential asymptotic analysison a singularly perturbed diatomic Hertzian chain. In Section 4.4, we compare our asymptotic and numericalresults for the Hertzian chain. From this comparison, we see that the employed exponential asymptotic methodis useful for approximating the solution behavior for a wide range of mass ratios, but it fails to detect someimportant features that arise from the interaction between oscillatory contributions. In Section 5, we concludeand further discuss our results.
We aim to determine the asymptotic behavior of the exponentially small oscillations in the wake of the solitary-wave front in singularly perturbed variants of a Hertzian chain. We use η to denote the small parameter,which is related to small mass ratio and satisfies η = m /m . It is impossible to determine the form of theseoscillations using only classical asymptotic power series, because the oscillation amplitude is exponentiallysmall in the asymptotic limit η →
0. In particular, the amplitude is smaller than any algebraic power of η .Therefore, we use an exponential asymptotic approach to study the exponentially small oscillations.Consider a singularly perturbed differential equation of the form F ( x, g ( x ) , g (cid:48) ( x ) , g (cid:48)(cid:48) ( x ) , . . . ; η ) = 0 , (3)where η is a small parameter. To find an asymptotic solution to (3), one typically expands a solution aboutsome leading-order approximation as a power series in the parameter η . Once one has obtained this leading-order solution, one analytically continues it into the complex plane. Typically, this solution is singular atsome set of points, which are the end points of Stokes curves [81]. Stokes curves play an important role in thepresent study: as one crosses a Stokes curve, the amplitude of the exponentially small contribution undergoesa smooth but rapid change in value. In the present work, we require the exponentially small contributions tovanish on one side of the Stokes curve, so they “switch on” as the Stokes curves are crossed. n general, an asymptotic power series of the solution to a singularly perturbed problem is divergent [24].Nevertheless, one can truncate a divergent asymptotic series to approximate the solution [16]. If one choosesthe truncation point to minimize the difference between the approximation and the exact solution, the approx-imation error is typically exponentially small in the asymptotic limit [15]. This is called “optimal truncation”.We use N opt to denote an optimal truncation point. We then express the solution as the sum of an optimallytruncated power series and an exponentially small contribution. Substituting this sum into (3) produces anequation that describes the behavior of the exponentially small contribution.Early work that applied this idea include [7, 8], which determined the Stokes switching behavior in severalimportant special functions. The analysis in [8] demonstrated that the switching behavior depends predictablyon the manner in which the asymptotic power series diverges. Subsequent work in [10] established techniquesknown as “hyperasymptotics” to further reduce the exponentially small error that is generated by truncatingthe series. See [9] for a summary and discussion of results in [8, 10].In the present study, we apply an exponential asymptotic method that was developed in [17, 69]. Weexpress the solution g of governing equations as an asymptotic power series g ∼ ∞ (cid:88) j =0 η rj g j as η → , (4)where r is the number of times that one needs to differentiate g j − to obtain g j .We first substitute the series (4) into the governing equations (3), and we then asymptotically match termsin the subsequent expression to obtain a recursion relation for g j . In a singularly perturbed problem, obtaining g j using the recursion relation typically requires differentiating earlier terms in the series each time one appliesthe recursion relation. If the series terms have singular points, this repeated differentiation ensures that theseries terms diverge in a predictable fashion called a factorial-over-power divergence [24]. As j becomes large,behavior of this form dominates the series terms.To capture this divergence, we write an ansatz for the behavior of g j in the limit that j → ∞ . This yieldsthe so-called “late-order terms” of an asymptotic series. As j → ∞ , factorial-over-power divergent behaviordominates the series terms. Motivated by this pattern, Chapman et. al [17] proposed applying a late-orderansatz to approximate the form of the late-order terms, even if computing earlier terms in a series is challengingor intractable. This ansatz is a sum of terms of the form g j ∼ G Γ( rj + γ ) χ rj + γ as j → ∞ , (5)where the parameter γ is constant and G and χ are functions of any variables but are independent of j . Thefunction χ is known as the “singulant” and it is equal to 0 at one or more singularities of the leading-ordersolution g . This ensures that the late-order ansatz for g j is singular at the same location as the leading order,with a singularity strength that grows as j increases. Substituting (5) into the differential equation (3) andmatching orders of j allows one to determine the functional forms of χ and G . By comparing the ansatz (5)with the local behavior of g in the neighborhood of singular points, one can calculate γ .As we noted earlier, one can optimally truncate the power series (4) using the additional information fromthe late-order ansatz (5). Optimal truncation points typically occur after an asymptotically large number ofterms, so one can apply the late-order ansatz (5) to obtain a simplified expression using the heuristic thatwas described in [15]. This heuristic requires truncating the series after the value of N opt for which the term η N opt g N opt has the smallest magnitude.One can use the truncated series to approximate the exact solution with exponentially small error. Onewrites g = N opt − (cid:88) j =0 η rj g j + g exp , (6)where g exp denotes the exponentially small error term and N opt is the optimal truncation point that onecalculates using the late-order ansatz (5).Substituting the truncated series expression (6) into the differential equation (3) produces an equation forthe exponentially small remainder term [69]. Away from Stokes curves, one can determine this remainder sing a straightforward application of the WKB method [36]. In the neighborhood of the Stokes curves, weapply an exponential ansatz for the exponentially small asymptotic terms g exp and write g exp ∼ S G e − χ/η as η → , (7)where S is a function that is known as the “Stokes multiplier”. The Stokes multiplier is approximately constantexcept in the neighborhood of Stokes curves, where it varies rapidly in a neighborhood of width O ( √ η ) as η →
0. This behavior is known as “Stokes switching”. The locations of Stokes curves are determined completelyby form of χ . As shown in [10], Stokes curves occur only in locations where χ is real and positive. One can useexponential asymptotics to directly calculate the exponentially small contributions to the asymptotic behaviorof solutions that appear as the Stokes curves are crossed. One cannot express these contributions using classicalasymptotic power series.The exponential asymptotic approach that we have described often requires the explicit calculation of onlythe leading-order solution in (4) to determine the exponentially small contributions. This is convenient, ascomputing series terms beyond a leading-order expression in nonlinear problems can be complicated or evenintractable. See [14,15] for more details on exponential asymptotics and their applications to nonlocal solitarywaves, [8, 9] for examples of other studies of exponential asymptotics, and [17, 69] for more details on theparticular methodology that we apply in the present paper. We consider a precompressed woodpile chain, such as the one in Figure 2(a). We model it by a particle chain;each particle has mass m and is connected to an outside particle of mass m by a linear spring with springconstant k . The governing equations of this model, which we illustrated in Figure 2(b), are m ¨ u ( n, t ) = [ δ + u ( n − , t ) − u ( n, t )] α + − [ δ + u ( n, t ) − u ( n + 1 , t )] α + − k [ u ( n, t ) − v ( n, t )] , (8) m ¨ v ( n, t ) = k [ u ( n, t ) − v ( n, t )] , (9)where δ is the precompression and u ( n, t ) and v ( n, t ), respectively, denote the displacements of the n -thparticle of mass m and m at time t . The choice α = 3 / α for which the leading-order solution can be approximatedby a soliton solution of the KdV equation. In practice, this corresponds to α > u = m ˆ u and v = m ˆ v , and we rewrite the governing equations in termsof ˆ δ = δ /m and ˆ k = k/m . Setting η = m /m gives the following scaled governing equations:¨ˆ u ( n, t ) = [ˆ δ + ˆ u ( n − , t ) − ˆ u ( n, t )] α + − [ δ + ˆ u ( n, t ) − ˆ u ( n + 1 , t )] α + − ˆ k [ˆ u ( n, t ) − ˆ v ( n, t )] , (10) η ¨ˆ v ( n, t ) = ˆ k [ˆ u ( n, t ) − ˆ v ( n, t )] . (11)We consider a precompressed chain in which the argument inside a bracket is always non-negative, as theparticles always maintain physical contact. Therefore, we omit the subscript + in all subsequent analysis. Forconvenience, we perform our analysis on the scaled system (10)–(11); in our subsequent notation, we omit thehats that indicate this scaling.We are interested in asymptotic behavior of traveling-wave solutions of (10)–(11) when 0 < η (cid:28)
1. Thesesolutions consist of a localized wave core, which we generate using an asymptotic power series in η , andexponentially small non-decaying oscillations that we compute using exponential asymptotic techniques. We expand u ( x, t ) and v ( x, t ) in an asymptotic power series in η in the limit η → u ( n, t ) ∼ ∞ (cid:88) j =0 η j u j ( n, t ) , v ( n, t ) ∼ ∞ (cid:88) j =0 η j v j ( n, t ) . (12) o find a leading-order solitary wave, which is the first step to construct a nanopteron solution of (10)–(11),we apply the series expression (12) to (11) and match at leading order in the limit η → u ( n, t ) = v ( n, t ) . (13)Inserting the relation (13) into (10) gives¨ u ( n, t ) = [ δ + u ( n − , t ) − u ( n, t )] α − [ δ + u ( n, t ) − u ( n + 1 , t )] α . (14)One can obtain an approximation to the leading-order behavior by following the analysis in [66]. We reportthe key steps of this analysis. Assume that the deformation in the chain is sufficiently small so that | u ( n − , t ) − u ( n, t ) | (cid:28) δ . We now write (14) as¨ u ( n, t ) = αδ α − [ u ( n − , t ) − u ( n, t ) + u ( n + 1 , t )]+ α ( α − δ α − [( u ( n − , t ) − u ( n, t )) − ( u ( n, t ) − u ( n + 1 , t )) ]+ O (( | u ( n − , t ) − u ( n, t ) | /δ ) ) . (15)In the long-wavelength limit, in which the characteristic size of a wave is large in comparison to the particleradius R , we obtain the KdV equation from (15), ∂w∂t + c ∂w∂x + γ ∂ w∂x + σ c w ∂w∂x = 0 , w = − ∂u ∂x , (16)where c = 2 R √ αδ α − , γ = c R , σ = 2( α − c Rδ . (17)Solutions to (16) approximate the leading-order behavior of u ( x, t ), and one can obtain the leading-orderbehavior of v ( x, t ) using (13). We are particularly interested in solutions that are perturbations of a leading-order traveling wave, and we therefore select w ( x, t ) in (16) to be the KdV traveling-wave solution. We obtainthe leading-order behavior u ( x, t ) by integration: u ( n, t ) = − δ(cid:15)α − (cid:15) ( n − c (cid:15) t )) + O ( (cid:15) / ) , (18)where (cid:15) = 2 RL , c (cid:15) = δ ( α − / √ α + √ α (cid:15) δ ( α − / . (19)The long-wavelength limit corresponds to (cid:15) →
0, and we will refer to (cid:15) as the “long-wavelength parameter”.Motivated by (18), we use the co-moving frame ξ = n − c (cid:15) t . In terms of ξ , the leading-order behavior is u ( ξ ) = − δ(cid:15)α − (cid:15)ξ ) + O ( (cid:15) / ) as (cid:15) → , (20) v ( ξ ) = u ( ξ ) . (21)The leading-order behavior has singularities at ξ N ± = ± i(2 N + 1) π/ (cid:15) for N ∈ Z . The singularities thatare closest to the real axis occur at ξ ± = ± i π/ (cid:15) and are associated with singulants with the smallest valuesof | χ | when evaluated at real ξ . Consequently, these terms dominate the late-order behavior as n → ∞ . Nearthese singularities, we calculate that u ( ξ ) ∼ − δ ( α − − ξ − ξ ± as ξ → ξ ± (22)and that v has the same behavior as u in the neighborhood of the singularity. Writing the governing equations (10)–(11) in terms of ξ and matching at each order of η gives c (cid:15) u (cid:48)(cid:48) j ( ξ ) = [ u j ( ξ − − u j ( ξ ))( δ + u ( ξ − − u ( ξ )] / − [ u j ( ξ ) − u j ( ξ + 1)][ δ + u ( ξ ) − u ( ξ + 1)] / − k [ u j ( ξ ) − v j ( ξ )] + . . . , (23) c (cid:15) v (cid:48)(cid:48) j − ( ξ ) = k [ u j ( ξ ) − v j ( ξ )] , (24) here we omit the terms that are products that include both u j − k and v j − k with k >
1. These terms aresubdominant with respect to the terms that we retain in the limit j → ∞ . We retain all terms that include u j , v j , and derivatives of v j − . Given the general form of the factorial-over-power ansatz (5), we conclude thatthe omitted terms do not contribute to the behavior of the asymptotic solution in our subsequent analysis.We will confirm this claim explicitly once we obtain the form of the late-order ansatz (25).In principle, one can apply (23) and (24) recursively to obtain terms in the series (12) up to arbitrarily largevalues of j . This process is challenging technically because it requires solving both a differential–differenceequation (23) and an algebraic equation (24) at each order. Fortunately, this is not necessary for our analysis;obtaining terms up to arbitrary order does not reveal the presence of oscillations in the far field, where ξ → −∞ , because the oscillations are exponentially small in the singularly perturbed limit η →
0. Instead,we obtain the asymptotic form of the late-order terms as part of the exponential asymptotic process that weuse to calculate the behavior of these oscillations.The late-order ansatz consists of a sum of terms of the forms u j ∼ U ( ξ )Γ(2 j + β ) χ ( ξ ) j + β , v j ∼ V ( ξ )Γ(2 j + β ) χ ( ξ ) j + β as j → ∞ . (25)To ensure that late-order terms have singularities at the same locations as the leading-order solution, we set χ = 0 at a particular choice of ξ = ξ ± . The full late-order expression is then the sum of the late-ordercontributions from each of the singularities.For sufficiently large j , the terms of the asymptotic series (25) diverge in a factorial-over-power fashion,which confirms that u j (cid:29) u j − k and v j (cid:29) v j − k as j → ∞ for k >
0. Additionally, inserting the late-orderansatz (25) into (23) shows that only β + 2 = β produces a nontrivial asymptotic balance, implying that u j = O ( v j − ) as j → ∞ and hence that v j (cid:29) u j as j → ∞ .Inserting the late-order ansatz (25) into (24) gives c (cid:15) ( χ (cid:48) ( ξ )) V ( ξ )Γ(2 j + β ) χ ( ξ ) j + β − c (cid:15) χ (cid:48) ( ξ ) V (cid:48) ( ξ )Γ(2 j + β − χ ( ξ ) j + β − − c (cid:15) χ (cid:48)(cid:48) ( ξ ) V ( ξ )Γ(2 j + β − χ ( ξ ) j + β − + · · · = − kV ( ξ )Γ(2 j + β ) χ ( ξ ) j + β + · · · , (26)where the omitted terms are no larger than O ( v j − ) in the j → ∞ limit.Matching terms at O ( v j ) in this limit gives the singulant equation c (cid:15) ( χ (cid:48) ( ξ )) = − k . This implies that χ (cid:48) ( ξ ) = ± i √ k/c (cid:15) , which we integrate to obtain χ ( ξ ) = ± i √ k ( ξ − ξ ± ) c (cid:15) . (27)Stokes switching can occur only if Re( χ ) >
0, which corresponds to the positive sign choice for ξ andthe negative sign choice for ξ − . Consequently, we retain these solutions and ignore the solutions that areassociated with the remaining sign choices, as those can never appear in the asymptotic solution.Matching terms at O ( v (cid:48) j − ) gives the prefactor equation 2 V (cid:48) ( ξ ) χ (cid:48) ( ξ ) = 0, so V is a constant, with a valuethat depends on the choice of singularity. For clarity, we subsequently use Λ ± to denote the constant prefactorthat is associated with a singularity at ξ = ξ ± . For the singular late-order behavior to be consistent with thelocal behavior of the leading-order solution (22), we calculate that β = 1. In Appendix A.1, we perform alocal expansion of the solutions u ( ξ ) and v ( ξ ) in a neighborhood of size O ( (cid:15) ) near the singularity, and we useasymptotic matching to determine that Λ + = − i δ √ k ( α − c (cid:15) . (28)We thereby fully determine the asymptotic behavior of u j in the j → ∞ limit.In Section 3.3, we give a detailed calculation of the exponentially small oscillation that is associated withthe singularity at ξ = ξ . At the conclusion of our analysis in Section 3.3, we state the correspondingcontribution that arises from the singularity at ξ = ξ − ; this contribution is the complex conjugate of thecontribution from the first singularity. .3 Stokes switching We truncate the asymptotic series after N terms to obtain u ( ξ ) = N − (cid:88) j =0 η j u j ( ξ ) + S N ( ξ ) , v ( ξ ) = N − (cid:88) j =0 η j v j ( ξ ) + R N ( ξ ) , (29)where S N and R N are the remainder terms that we obtain by truncating the series. The remainder termsare exponentially small if we optimally truncate the series; we again denote the optimal truncation point as N = N opt . To do this truncation, we apply the heuristic that was described in [15] and determine the optimaltruncation point by finding the point at which consecutive terms in the series are equal in size. This heuristicgives N opt = | χ | / η + ω , where we choose ω ∈ [0 ,
1) in a way that ensures that N opt is an integer. Note that N opt → ∞ as η → c (cid:15) S (cid:48)(cid:48) ( ξ ) ∼ − kR N ( ξ ) , (30) η c (cid:15) R (cid:48)(cid:48) ( ξ ) + η N c (cid:15) v (cid:48)(cid:48) N − ( ξ ) ∼ − kR N ( ξ ) as η → , (31)where omitted terms are smaller than those that we retain in the limit η →
0. Equation (31) decouples from(30), so we can study it independently. Applying the late-order ansatz and rearranging gives η c (cid:15) R (cid:48)(cid:48) N + kR N ∼ − Λ η N ( χ (cid:48) ) Γ(2 N + 1) χ N +1 as η → . (32)The right-hand side of (32) is exponentially small except in a neighborhood around the Stokes curve. Awayfrom the Stokes curve, we use the WKB method to obtain R N ∼ C e − χ/η as η → , (33)where C is a constant that we have not yet determined.To capture the variation in the neighborhood of a Stokes curve, we adapt the form of (33) to include aStokes-switching parameter A ( ξ ), which is constant except in a region in the neighborhood of the Stokes curve.In such a neighborhood, R N takes the form R N ( ξ ) ∼ A ( ξ )e − χ/η as η → . (34)Inserting (34) into (32) and rearranging yieldsd A d ξ ∼ Λ χ (cid:48) η N − Γ(2 N + 1)2 χ N +1 e χ/η as η → . (35)Writing N opt in terms of χ , expanding the gamma function using Stirling’s formula, and transforming to make χ the independent variable givesd A d χ ∼ Λ √ πη | χ | /η +2 ω − ( | χ | /η ) | χ | /η +2 ω − / e −| χ | /η √ χ | χ | /η +2 ω +1 e χ/η as η → . (36)We transform (36) into polar coordinates using χ = ρ e i θ and consider variations in the angular direction.After some simplification, we obtaind A d θ ∼ iΛ (cid:114) πρ η exp (cid:18) ρη (e i θ − − i θρη − ωθ (cid:19) as η → . (37)The right-hand side of (37) is exponentially small in η , except in the neighborhood of θ = 0. Defining an innerregion θ = η / ¯ θ , we find that d A d¯ θ ∼ iΛ η (cid:114) πρ − ρ ¯ θ / . (38)By integrating (38), we see that the behavior of A as the Stokes curve is crossed is A ∼ iΛ η (cid:114) π (cid:90) √ ρ ¯ θ −∞ e − s / d s as η → . (39) e evaluate the integral in (39) and find that the difference between the values of A on the two sides of theStokes curve is [ A ] + − ∼ i π Λ η as η → , (40)where [ A ] + − denotes the change in A as the Stokes curve is crossed from θ < θ >
0. Recalling (28) and (34),we find that the exponentially small contribution from ξ = ξ is[ R N ] + − ∼ δπ √ k ( α − c (cid:15) η e − i √ k ( ξ − i π/ (cid:15) ) /c (cid:15) η as η → . (41)The exponentially small contribution from ξ = ξ − is given by the complex conjugate of (41). Therefore, thetotal exponentially small contribution is[ R N ] + − ∼ δπ √ k ( α − c (cid:15) η e − i √ k ( ξ − i π/ (cid:15) ) /c (cid:15) η + c . c . as η → , (42)where c.c. denotes the complex conjugate. We express (42) in terms of trigonometric functions to obtain[ R N ] + − ∼ δπ √ k ( α − ηc (cid:15) exp (cid:32) − π √ k c (cid:15) (cid:15)η (cid:33) cos (cid:32) ξ √ kc (cid:15) η (cid:33) as η → . (43) We compare the amplitude of the oscillations that we predict using our asymptotic analysis to those fromnumerical simulations. We employ a symplectic integrator in the form of a velocity Verlet algorithm [2, 87],which is convenient for studying chains of particles; it has also been used to study molecular dynamics. Thisapproach is designed to conserve the energy of a system, which is not true of many common numerical methods,such as Runge–Kutta algorithms. The velocity Verlet algorithm uses the discretization x n +1 = x n + v n ∆ t + 12 a n (∆ t ) , v n +1 = v n + 12 ( a n + a n +1 )∆ t , (44)where ∆ t is the size of the time step and x n , v n , and a n are vector quantities that encode the displacements,velocities, and accelerations of the particles at time t = t n . In a woodpile chain, we obtain the acceleration a n as an algebraic function of x n using (10)–(11).We truncate the domain and impose periodic boundary conditions. Instead of directly computing theabsolute displacements and velocities of the particles, we perform our computations using the relative dis-placements and velocities, which we compute by calculating the differences in the positions and velocities ofadjacent particles. This is convenient for our computations because the far-field behavior of the leading-ordertraveling wave approaches 0 in both directions in these coordinates. The relative displacements of the heavyparticles and light particles are r and r , respectively, where r ( n, t ) = u ( n + 1 , t ) − u ( n, t ) , r ( n, t ) = v ( n + 1 , t ) − v ( n, t ) . (45)In relative coordinates, the governing equations are¨ r ( n, t ) = 2[ δ − r ( n, t )] α + − [ δ − r ( n − , t )] α + − [ δ − r ( n + 1 , t )] α + − k [ r ( n, t ) − r ( n, t )] , (46) η ¨ r ( n, t ) = k [ r ( n, t ) − r ( n, t )] . (47)The truncated domain has M = 2 particles, whose indices are n ∈ {− M/ , . . . , M/ } . The initialcondition is the leading-order traveling-wave solution (20)–(21). To avoid interactions between the far-fieldoscillations and the leading-order traveling wave in the periodic domain, we apply a window function to thesolution, as in [31,55], at each iteration. We obtain the window by multiplying r ( n, t ) and ˙ r ( n, t ) by a function W ( n − n max + M/ W ( k ) = , | K | ≤ M − M (cid:0) | K | − M (cid:1) , M < | K | ≤ M , M < | K | ≤ M (48) − − − t . . . r ( n, t )UndisturbedWave front Exponentially small oscillations Figure 4: Numerically-calculated relative displacement for a single particle at some index n in a woodpile chain with α = 1 .
5. Themass-ratio parameter of the depicted particle is η = 0 .
4, the precompression of the chain is δ = 5, and the long-wavelengthparameter is (cid:15) = 0 . and n max denotes the location of the maximum of the leading-order solution.In [31], it was argued that a window of this form cannot affect the behavior of the main wave or trailingoscillations in simulations of FPUT systems, as any disturbances that are associated with the windowing musttravel more slowly than the leading-order solitary wave and hence than the window itself (see [29]). We donot perform a comparable analysis for a woodpile chain, but we do not observe any discernible differences inthe amplitudes of the trailing oscillations as a result of our windowing.In Figure 4, we show a numerically-calculated profile of one particle in a woodpile chain. In this simulation,we use a mass-ratio parameter of η = 0 .
4, a precompression of δ = 5, and long-wavelength parameter of (cid:15) = 0 . t . The particleis initially at rest before the leading-order solitary wave reaches the particle, causing it to become displaced.Once this solitary wave has passed the particle, the particle continues to oscillate with an exponentially smallamplitude. In Figure 5, we measure this oscillation amplitude from numerical simulations and compare it toour asymptotic approximations for a range of parameter choices.In our numerical simulations, we calculate the amplitude of the far-field waves for δ = 5, δ = 7 . δ = 10,and δ = 15 with long-wavelength parameters of (cid:15) = 0 . (cid:15) = 0 .
25, and (cid:15) = 0 .
3. We choose the values of η sothat the far-field waves large enough to be detected numerically, but small enough that our small- η asymptoticapproximation can accurately describe the behavior of the solution. We run each of our simulations for enoughtime so that any time variation in the far-field oscillation amplitude is not visible in the solution. We showthe results of these simulations, which we compare to the asymptotic results from (43) in Figure 5.It is apparent from Figure 5 that our asymptotic approximation captures the far-field wave amplitudeeffectively for a wide range of the parameters and that the approximation error decreases as η →
0, aswe expect from the asymptotic nature of the approximation. For a fixed (cid:15) , our asymptotic approximationbecomes more accurate for progressively larger δ . For progressively smaller δ , the system approaches thestrongly nonlinear regime and the long-wavelength approximation becomes less accurate, which in turn causesthe late-order approximation to become less accurate. Our asymptotic approximation also appears to becomemore accurate for progressively larger α , but the reason for this improvement is not apparent from our analysis.By comparing our numerical computations and our asymptotic results, we conclude that we can explain theoscillations that appear in woodpile chains as a consequence of Stokes phenomena and that the exponentialasymptotic analysis that we described in this section is a useful technique for studying nanoptera in thesechains. = 5 − − − − − − − − − − .
15 0 . . l og ( A m p li t ud e ) α = 1 . η δ = 7 . − − − − − − − − − − .
15 0 . . η δ = 10 − − − − − − − − − − .
15 0 . . η δ = 15 − − − − − − − − − −
11 0 .
15 0 . η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . δ = 5 − − − − − − − − − −
11 0 . . . l og ( A m p li t ud e ) η α = 2 δ = 7 . − − − − − − − − − −
11 0 . . . η δ = 10 − − − − − − − − − − . . . η δ = 15 − − − − − − − − − − . . . η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . δ = 5 − − − − − − − − − − . . . l og ( A m p li t ud e ) η α = 2 . δ = 7 . − − − − − − − − − − . . η δ = 10 − − − − − − − − − − . . η δ = 15 − − − − − − − − − − . . η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . Figure 5: Comparison of our asymptotic approximations and numerical simulations of the amplitude of exponentially small oscillationsin the far field (i.e., for ξ → −∞ ) for woodpile chains with interaction exponents of (a) α = 1 .
5, (b) α = 2, and (c) α = 2 . α , we show simulations for δ = 5, δ = 7 . δ = 10, and δ = 15. For each choice of δ , we set the long-wavelengthparameter to be (cid:15) = 0 . (cid:15) = 0 .
25, and (cid:15) = 0 .
3. The accuracy of our asymptotic approximation improves for progressivelysmaller values of η , as expected for an η → δ . This is a consequence of the use of a leading-order approximation for the weakly nonlinear regime; decreasing δ approaches a configuration that needs to be treated as strongly nonlinear. Our approximation also improves as we increase α , although the reason for this improvement is not apparent. Singularly Perturbed Diatomic Hertzian Chain
We now consider a singularly perturbed diatomic Hertzian chain. This problem is related to the woodpilechains that we analyzed in Section 3. Once we derive a leading-order solution for a traveling wave in adiatomic Hertzian chain, we can perform a similar exponential asymptotic analysis. Obtaining the singulantequation for diatomic Hertzian chains requires solving a complicated integral expression; this is significantlymore complex than the corresponding step in our analysis of woodpile chains.We consider a diatomic Hertzian chain, which is governed by (49)–(50). We denote the masses in the chainby m ( j ) = m for even j and m ( j ) = m for odd j . We consider a small ratio between the masses of the twotypes of particles. This choice amounts to a precompressed Hertzian analog of the diatomic Toda and FPUTchains that were studied in [55, 56]. Those previous studies used far-field asymptotic behavior to identify anorthogonality condition, which yields values of the mass-ratio parameter η that cause the far-field waves tocancel and produce a genuine solitary wave.We derive an asymptotic approximation in a diatomic Hertzian chain to identify an orthogonality conditionthat finds some some (but not all) values of η that produce solitary waves. Specifically, our asymptotic approachsystematically fails to identify every second value of η that causes the far-field oscillations to cancel.Much of our analysis in this section is similar to that in Section 3 and in [55, 56]. Therefore, we presentonly an outline of our calculations.The governing equations for a diatomic chain of particles is m ¨ u (2 n, t ) = [ δ + v (2 n − , t ) − u (2 n, t )] α + − [ δ + u (2 n, t ) − v (2 n + 1 , t )] α + , (49) m ¨ v (2 n + 1 , t ) = [ δ + u (2 n, t ) − v (2 n + 1 , t )] α + − [ δ + v (2 n + 1 , t ) − u (2 n + 2 , t )] α + , (50)where u and v , respectively, represent the displacement of even and odd particles, the + subscript has thesame meaning as in (10)–(11), and we set α = 3 / We apply thesame scalings as for woodpile chains, so we let u = m ˆ u and v = m ˆ v . We rewrite (49)–(50) in terms of ourscaled variables and a scaled parameter ˆ δ = δ /m . Setting η = m /m yields the scaled system¨ˆ u (2 n, t ) = [ˆ δ + ˆ v (2 n − , t ) − ˆ u (2 n, t )] / − [ˆ δ + ˆ u (2 n, t ) − ˆ v (2 n + 1 , t )] / , (51) η ¨ˆ v (2 n + 1 , t ) = [ˆ δ + ˆ u (2 n, t ) − ˆ v (2 n + 1 , t )] / − [ˆ δ + ˆ v (2 n + 1 , t ) − ˆ u (2 n + 2 , t )] / , (52)As in the woodpile chains, the system (51)–(52) is sufficiently precompressed that the quantities in the squarebrackets are never negative. We therefore omit the + subscript in our subsequent notation. In our analysis ofthe scaled system (51)–(52), we also omit the hats from the variables and parameters for notational convenience.As in the woodpile chains of Section 3, we study traveling-wave solutions to (51)–(52) when the mass-ratioparameter η is small. Specifically, we consider 0 < η (cid:28)
1. As before, the solution consists of a localized wavecore that depends algebraically on η and has exponentially-small, non-decaying oscillations. We again useexponential asymptotic methods to study these oscillations. The first step to construct a nanopteron solution of (51)–(52) is to find a leading-order solitary wave. Therefore,we begin by expanding the solution u ( x, t ) and v ( x, t ) using an asymptotic power series in η in the limit η → u ( n, t ) ∼ ∞ (cid:88) j =0 η j u j ( n, t ) , v ( n, t ) ∼ ∞ (cid:88) j =0 η j v j ( n, t ) . (53)Inserting the series expansions (53) into (52) and matching at leading order in the η → v (2 n + 1 , t ) = [ u (2 n, t ) + u (2 n + 2 , t )] . (54) In our analysis of woodpile chains in Section 3, we permitted α to take any value that produces leading-order solitary waves, asthe choice of α did not significantly impact the analysis. For diatomic chains, the choice of α has a more significant effect on theexponential asymptotic analysis, so we restrict our attention to Hertzian interactions, which is the most physically significant case. e( ξ )Im( ξ ) ξ − = i π(cid:15) − ξ = i π(cid:15) + 1 ξ − = − i π(cid:15) − ξ = − i π(cid:15) + 1 Figure 6: Singularities of v ( ξ ) in (59) that contribute to the asymptotic form of the far-field oscillations of (84). We give a detaileddiscussion of the analysis for the oscillations that arise from the singularity at ξ = ξ − . We then state the contributionsfrom the other three singularities. Inserting (54) into (51) yields¨ u (2 n, t ) = (cid:18) δ + u (2 n − , t ) − u (2 n, t )2 (cid:19) / − (cid:18) δ + u (2 n, t ) − u (2 n + 2 , t )2 (cid:19) / . (55)Using a similar analysis as in Section 3.1, we obtain an approximation to the leading-order solution byassuming that the displacement of a particle is small in comparison to the precompression parameter δ andthat the characteristic length scale of the wave is large in comparison to the radius of the spherical particles.We again define a long-wavelength parameter (cid:15) , and we find that the leading-order solution is approximately u ( ξ ) = − δ(cid:15) tanh ( (cid:15)ξ/
2) + O ( (cid:15) / ) , v ( ξ ) = ( u ( ξ + 1) + u ( ξ − , (56)where the variable ξ defines the moving frame ξ = n − c (cid:15) t , c (cid:15) = δ / √ δ / (cid:15) √ . (57)We now extend the leading-order solution (56) into the complex plane. We see that v ( ξ ) is singular at ξ = ξ N ± in the complex plane, where ξ N ± = (2 N − π i (cid:15) ± , N ∈ Z . (58)The behavior near the singularities is v ( ξ ) ∼ − δ ( ξ − ξ N ± ) − as ξ → ξ N ± . (59)As with woodpile chains, we find that the dominant contribution to the late-order behavior arises from thesingularities that are closest to the real axis; this occurs when N = 0 and N = 1 in (58). We illustrate thepositions of these singularities in Figure 6. In our subsequent analysis, we show the details for calculating thecontributions from the singularity at ξ = ξ − . We then state the analogous results for the contributions fromthe singularities at ξ = ξ , ξ − , and ξ . .2 Terms in the late-order series Inserting the series (53) into the system (51)–(52) and matching at each order of η gives a recurrence relationfor j ≥
2. The recurrence relation is c (cid:15) u (cid:48)(cid:48) j ( ξ ) = ( v j ( ξ − − u j ( ξ ))( δ + v ( ξ − − u ( ξ )) / − ( u j ( ξ ) − v j ( ξ + 1))( δ + u ( ξ ) − v ( ξ + 1)) / + ( v j − ( ξ − − u j − ( ξ ))( v ( ξ − − u ( ξ ))[ δ + ( v ( ξ − − u ( ξ ))] − / − ( u j − ( ξ ) − v j − ( ξ + 1))( u ( ξ ) − v ( ξ + 1))[ δ + ( u ( ξ ) − v ( ξ + 1))] − / + · · · , (60) c (cid:15) v (cid:48)(cid:48) j − ( ξ ) = ( u j ( ξ − − v j ( ξ ))( δ + u ( ξ − − v ( ξ )) / − ( v j ( ξ ) − u j ( ξ + 1))( δ + v ( ξ ) − u ( ξ + 1)) / + ( u j − ( ξ − − v j − ( ξ ))( u ( ξ − − v ( ξ ))[ δ + ( u ( ξ − − v ( ξ ))] − / − ( v j − ( ξ ) − u j − ( ξ + 1))( v ( ξ ) − u ( ξ + 1))[ δ + ( v ( ξ ) − u ( ξ + 1))] − / + · · · , (61)where we omit the terms that are products that include both u j − k and v j − k with k >
1. The terms that weretain are those that contain u j − , v j − , u j , and v j . All omitted terms are subdominant in comparison to theretained terms as j → ∞ .We again use a factorial-over-power late-order ansatz to approximate the late-order terms, so we write u j ∼ U ( ξ )Γ(2 j + β ) χ ( ξ ) j + β , v j ∼ V ( ξ )Γ(2 j + β ) χ ( ξ ) j + β as j → ∞ , (62)where β and β are constants. Inserting the late-order ansatz (62) into (60) shows that only β = β − χ Inserting the ansatz (62) into the recurrence relation (60)–(61) and matching at O ( z j ) as j → ∞ gives c (cid:15) ( χ (cid:48) ) = − (cid:18) δ + u ( ξ − − u ( ξ + 1)2 (cid:19) / . (63)Integrating (63) and recalling that χ = 0 at the singularity location ξ = ξ − implies that χ = ± √ c (cid:15) (cid:90) C (cid:18) δ + ( u ( s − − u ( s + 1))2 (cid:19) / d s , (64)where C is a contour from ξ − to ξ . Although one can select any contour, it is convenient to divide the contourinto a vertical component C and a horizontal component C . We show our contour in Figure 7.The contribution to the integral in (64) from C has both real and imaginary components, but the contri-bution to the integral from C is purely imaginary. The real and imaginary parts of χ areRe( χ ) = ± Re (cid:26) √ c (cid:15) (cid:90) C (cid:18) δ + u ( s − − u ( s + 1)2 (cid:19) / d s (cid:27) , (65)Im( χ ) = ± Im (cid:26) √ c (cid:15) (cid:90) C (cid:18) δ + u ( s − − u ( s + 1)2 (cid:19) / d s (cid:27) ± √ c (cid:15) (cid:90) C (cid:18) δ + u ( s − − u ( s + 1)2 (cid:19) / d s , (66)where the sign choices are either all positive or all negative. We see that Re( χ ) is constant for real-valued ξ ,so the far-field oscillations have a constant amplitude. Because | u ( ξ − − u ( ξ + 1) | → | ξ | → ∞ , wealso see that χ (cid:48) ( ξ ) tends to a constant value and that Im( χ ) depends linearly on ξ in this limit. Consequently,the associated far-field oscillations must tend to a constant wavelength.Stokes switching can occur only if Re( χ ) >
0, which corresponds to the positive signs in both (65) and (66).Therefore, we restrict our analysis to this choice of sign. s )Im( s ) ξ − = i π(cid:15) − ξ C C Figure 7: The contour of integration for (64) that connects the singularity at s = ξ − with the point s = ξ . We indicate the locationof the singularity using a cross, and we show the contour using a thick black curve. We divide the contour into a verticalcomponent C and a horizontal component C . The contribution to the contour integral from C has both real and imaginarycomponents, whereas the contribution to the integral from C is purely imaginary. This implies that Re( χ ) is constant forreal-valued ξ . For our subsequent analysis, it is helpful to have an asymptotic expression for the local behavior of χ nearthe singular point ξ − . From direct calculation using local expressions for u ( ξ + 1) and u ( ξ − c (cid:15) ( χ (cid:48) ) ∼ − √ δ ( ξ − ξ − ) − / as ξ → ξ − . (67)We rearrange this local expression to give χ (cid:48) ∼ ± √ δ / i c (cid:15) ( ξ − ξ − ) − / as ξ → ξ − . (68)We choose the positive sign to be consistent with the full expression for χ , and we integrate to obtain χ ∼ √ δ / i3 c (cid:15) ( ξ − ξ − ) / as ξ → ξ − . (69) V and β Matching (60)–(61) at order O ( v (cid:48) j − ) as j → ∞ gives V = Λ − (cid:112) χ (cid:48) ( ξ ) , (70)where Λ − is a constant that we can determine by taking an inner expansion of the solution in the neighborhoodof ξ = ξ − and matching the outer limit of this expansion with the inner limit of the late-order ansatz (62).We need to determine the value of β in (62) to perform this inner analysis.To determine β , we combine (68) and (70) to obtain a local expression for V . This expression is V ∼ Λ − c / (cid:15) / δ / √ i ( ξ − ξ − ) / as ξ → ξ − . (71)Substituting (69) and (71) into (62) yields a local expression for the late-order ansatz: v j ( ξ ) ∼ Λ − c / (cid:15) ( ξ − ξ − ) / Γ(2 j + β )6 / δ / √ i(4 √ δ / i( ξ − ξ − ) / / c (cid:15) ) j + β . (72)We note that the leading-order solution v (56) has a singularity of order 1 at ξ = ξ − . For (72) to be consistentwith repeated differentiation of the leading-order behavior, we require that 3 β / − / β = 3 / − . The late-order expansion (62) breaks down in the neighborhood of ξ = ξ − in the region η χ − = O (1) as η →
0. By matching the inner limit of the late-order behavior [whichis given in (72)] with the outer limit of a local expansion near this singularity, it is possible to obtain Λ − . Weshow the computational portion of this procedure in Appendix A.2 and obtain c (cid:15) Λ − /δ / ≈ . .3 Calculations of remainders As in Section 3, we truncate the asymptotic series after N terms. This yields u ( ξ ) = N − (cid:88) j =0 η j u j ( ξ ) + S N ( ξ ) , v ( ξ ) = N − (cid:88) j =0 η j v j ( ξ ) + R N ( ξ ) , (73)where S N and R N are the remainder terms that we obtain by truncating the series. The optimal truncationpoint, which we find using the same method as in Section 3.3, is N opt = | χ | / η + ω , where we choose ω ∈ [0 , N opt is an integer. Note that N opt → ∞ as η → c (cid:15) η R (cid:48)(cid:48) N ( ξ ) − c (cid:15) χ (cid:48) ( ξ ) R N ( ξ ) ∼ − η N χ (cid:48) ( ξ ) V ( ξ )Γ(2 N + 3 / χ ( ξ ) N +3 / as η → χ ) = 0, which is the Stokes curve. To capture the behavior of the remainder near this Stokescurve, we write the remainder R N using the same adapted WKB ansatz as (34). That is, R N ( ξ ) ∼ A ( ξ )e − χ/η as η → , (75)where A ( ξ ) is a Stokes switching parameter that varies rapidly near the Stokes curve and is constant outsideof the rapidly varying neighborhood. Inserting (75) into (74) yields − A (cid:48) χ (cid:48) e − χ/η ∼ − η N − ( χ (cid:48) ) Γ(2 N + 3 / χ N +3 / . (76)We transform (76) to treat χ as the independent variable and write χ = r e iθ in polar coordinates. We fix ρ , and we consider only variations in the angular direction. Using the optimal truncation N = N opt andsimplifying using Stirling’s formula givesd A d θ ∼ i η (cid:114) πρ (cid:18) ρη (cid:16) e i θ − (cid:17) − i θrη − i θ (1 / ω ) (cid:19) . (77)The right-hand side of (77) is exponentially small in η , except in the neighborhood of θ = 0. Defining an innerregion θ = η / ¯ θ , we find that d A d¯ θ ∼ i (cid:114) πρ η e − r ¯ θ / . (78)By integrating (78), we see that the behavior of A as the Stokes curve is crossed is A ∼ i (cid:114) π η (cid:90) √ ρ ¯ θ −∞ e − s / d s . (79)We evaluate the integral in (79) and use the form of R N from (75) and find that the exponentially smallcontribution from ξ − is [ R N ] + − ∼ i π Λ η / √ χ (cid:48) e − χ/η . (80)To capture the change across the Stokes curve, we also need to include the contribution from ξ − ; this is givenby the complex conjugate of (80). Adding (80) and its complex conjugate, we find that the change in theexponentially small contribution as the Stokes curve is crossed from left to right is[ R N ] + − ∼ i π Λ η / √ χ (cid:48) e − χ/η + c . c . , (81)where “c.c.” denotes the contribution from singularity ξ − that is conjugated to ξ − . map: ”that is conjugated to”: what does this mean exactly? I am confused by this phrasing he overall exponentially small contribution to the asymptotic behavior of v ( ξ ) in the wake of the leading-order solitary wave is given by the sum of the contributions from each of the four singularities (see Fig. 6). Inthe limit that η →
0, the exponentially small terms have the asymptotic expression v exp ∼ i πη / e − Re( χ − ( ξ )) /η Λ − e − iIm( χ − ( ξ )) /η (cid:112) χ (cid:48) − ( ξ )+ i πη / e − Re( χ ( ξ )) /η Λ e − iIm( χ ( ξ )) /η (cid:112) χ (cid:48) ( ξ ) + c . c , (82)where χ − is the singulant that is associated with ξ = ξ − and χ is the singulant that is associated with ξ = ξ . Using the relation Re( χ ) = Re( χ − ) from (65), the asymptotic expression (82) becomes v exp ∼ i πη / e − Re( χ − ( ξ )) /η (cid:34) Λ − e − iIm( χ − ( ξ )) /η (cid:112) χ (cid:48) − ( ξ ) + Λ e − iIm( χ ( ξ )) /η (cid:112) χ (cid:48) ( ξ ) (cid:35) + c . c . (83)In Appendix A.2, we show that Λ = iΛ − . Far behind the central solitary wave, | u ( ξ − − u ( ξ + 1) | → χ (cid:48) , ± ( ξ ) ∼ i √ δ / /c (cid:15) in the ξ → −∞ limit, as indicated in (63). This yields a convenient trigonometric simplification: v exp ∼ π Λ − η / (cid:112) χ (cid:48) − ( ξ ) e − Re( χ − ( ξ )) /η cos (cid:18) Im( χ − ( ξ ) − χ ( ξ ))2 η + π (cid:19) × cos (cid:18) Im( χ ( ξ ) + χ − ( ξ ))2 η − π (cid:19) . (84)From the integral expression (64) for the singulant χ , we see that Re( χ , ± ) and Im( χ − ) − Im( χ ) do notdepend on ξ . Consequently, we can write the amplitude of the far-field waves as v amp ∼ πc / (cid:15) Λ − / δ / η / e − Re( χ − ( ξ )) /η cos (cid:18) Im( χ − ( ξ ) − χ ( ξ ))2 η + π (cid:19) as η → . (85) As in Section 3.4, we compare our asymptotic results for the far-field amplitude to numerical simulations. Weagain use the velocity Verlet algorithm and again simulate the equations of motion on a large periodic domainwith relative coordinates. For diatomic Hertzian chains, the relative coordinate system is given in terms of r ( n, t ), where r (2 n, t ) = v (2 n + 1 , t ) − u (2 n, t ) , r (2 n + 1 , t ) = u (2 n + 2 , t ) − v (2 n + 1 , t ) . (86)In this coordinate system, the equations of motion (51)–(52) become¨ r (2 n, t ) = (cid:18) η (cid:19) [ δ − r (2 n, t )] / − [ δ − r (2 n − , t )] / − η [ δ − r (2 n + 1 , t )] / , ¨ r (2 n + 1 , t ) = (cid:18) η (cid:19) [ δ − r (2 n + 1 , t )] / − η [ δ − r (2 n, t )] / − [ δ − r (2 n + 2 , t )] / . As with the woordpile chain, the domain includes M = 2 particles, with indices n ∈ {− M/ , . . . , M/ } .At each time step, we multiply the solution by the windowing function in (48). The initial condition is givenby the leading-order traveling-wave solution (56).We calculate the amplitude of the far-field waves for the precompressions δ = 5, δ = 7 . δ = 10, and δ = 15with long-wavelength parameters (cid:15) = 0 . (cid:15) = 0 .
5, and (cid:15) = 0 . η . As with woodpile chains, we run each simulation for a sufficiently long time so that the amplitude of thefar-field oscillations appear to reach a constant value. In Figure 8, we present the results of our computationsand compare them to our asymptotic results from (85).We make two observations from Figure 8. First, our asymptotic approximation accurately predicts theamplitude of the far-field waves for a large range of values of η . Second, our asymptotic calculation does notalways identify the values of η for which destructive interference causes the far-field oscillations to vanish. In − − − −
12 0 . . . .
15 0 .
35 0 . l og ( A m p li t ud e ) η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (a) δ = 3 − − − − −
12 0 . . . .
15 0 .
35 0 . l og ( A m p li t ud e ) η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (b) δ = 5 − − − − −
12 0 . . . .
15 0 .
35 0 . l og ( A m p li t ud e ) η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (c) δ = 10 − − − − −
12 0 . . . .
15 0 .
35 0 . l og ( A m p li t ud e ) η (cid:15) = 0 . (cid:15) = 0 . (cid:15) = 0 . (d) δ = 15Figure 8: Comparison of our asymptotic results and numerical simulations of the far-field wave amplitude of nanoptera in a diatomicHertzian chain with precompressions of (a) δ = 3, (b) δ = 5, (c) δ = 10, and (d) δ = 15. For each value of δ , we show ourcomputations for long-wavelength parameters of (cid:15) = 0 . (cid:15) = 0 .
5, and (cid:15) = 0 .
6. Our asymptotic prediction for the far-fieldamplitude is accurate for a wide range of values of the mass-ratio parameter η , with certain important exceptions. For fixedvalues of δ and (cid:15) , two oscillatory wave trains interact and interfere destructively for specific values of η . Our asymptoticapproximation is capable of predicting the destructive interference that occurs for η ≈ .
32 and η ≈ .
22, and the accuracyof the prediction becomes better for progressively smaller values of η . However, the approximation does not predict thedestructive interference that occurs for η ≈ .
25 and η ≈ .
19. Our approximation for the destructive interference becomesless accurate for progressively smaller (cid:15) . each panel of Figure 8, we observe four values of η (these are η ≈ . η ≈ . η ≈ .
22, and η ≈ . (cid:18) Im( χ − ( ξ ) − χ ( ξ ))2 η + π (cid:19) = 0 , (87)which simplifies to solving2 πη (cid:18) K + 14 (cid:19) = √ c (cid:15) (cid:90) ξ ξ − (cid:18) δ + ( u ( ξ − − u ( ξ + 1))2 (cid:19) / d s (88)for integer values of K . The expression (88) is an orthogonality condition [1, 42, 55, 56, 84]. Values of η thatsatisfy (88) correspond to dips in the asymptotically predicted amplitude in Figure 8. These values yieldlocalized solitary waves. We can see that the condition (88) is capable of predicting the localized solitarywaves that arise for η ≈ .
32 and η ≈ .
22. However, it fails to predict the localized solitary waves thatarise for η ≈ .
25 and η ≈ .
19. In fact, our exponential asymptotic analysis appears to systematically missevery second point (i.e., if it captures one point, then it misses the next point, and vice versa) at which theoscillations vanish. Additionally, for progressively smaller values of (cid:15) , the orthogonality condition (88) becomesless accurate at predicting the values of η that give localized solitary waves.Systematically missing half of these solitary-wave solutions appears to be a limitation of the method that wehave used to obtain the asymptotic contribution. There are several approximations in our analysis. Notably,we analytically continue the long-wavelength KdV approximation from the real axis into the complex plane.The orthogonality condition (88) depends on the behavior of u ( ξ ) on a contour that connects the singularities t ξ − and ξ (see Figure 6). These singularities are far from the real axis, and their distance from the realaxis increases as the long-wavelength parameter (cid:15) decreases. Consequently, it is possible that errors in thisapproximation — and, in particular, in the locations of the singular points — become larger for progressivelysmaller (cid:15) .Detecting destructive interference involves determining configurations in which two oscillations with veryshort wavelengths precisely cancel each other. The individual Stokes contributions have wavelengths that areof size O (1 /η ) as η →
0. Errors in the positions of the singularities on this very short scale are capable ofsignificantly changing the values of η for which cancellation occurs. It is possible that approximation errors inthe leading-order solution lead to shifts in the singularities in the analytically-continued solution. This maydisrupt the precise cancellation that is necessary for the manifestation of localized solitary-wave solutions,causing some of them to be missed by our asymptotic approximation.The issue of missing solitary-wave solutions did not arise in a recent study of a diatomic FPUT chain [55],despite the use of a similar long-wavelength approximation to continue a leading-order solution into thecomplex plane. It is not apparent why the method worked effectively for diatomic FPUT chains, but hasfailed to identify analogous behavior in our diatomic Hertzian chain. This is an interesting open question thatmerits further study. In the present study, we derived asymptotic approximations for traveling waves in two particle chains thatare singular perturbations of a monoatomic Hertzian particle chain with precompression. Specifically, weconsidered a woodpile chain and a diatomic Hertzian chain with a small ratio between the masses of itstwo types of particles. The traveling-wave solutions of these systems are nanoptera, which consist of anexponentially localized solitary wave along with periodic oscillations of exponentially small amplitude.A monoatomic Hertzian chain supports traveling solitary-wave solutions, which one can approximate assoliton solutions to the KdV equation under suitable assumptions. We used this monoatomic solitary waveas the leading-order behavior for both woodpile chains and diatomic Hertzian chains. In both systems, wefound that non-decaying oscillations appear in the wake of a primary traveling wave. These oscillations arisefrom the Stokes phenomenon, in which exponentially small contributions to an analytically-continued solutionswitch on or off as Stokes curves are crossed in the complex plane.Using exponential asymptotic analysis, we obtained an asymptotic form for these exponentially small con-tributions, and we used this form to determine the amplitude of the non-decaying far-field waves that appearbehind the traveling wave. For woodpile chains, the solution has one Stokes curve, which produces one expo-nentially small oscillation in the region behind the wave front. We determined an asymptotic approximationof this behavior, and we compared our asymptotic results with numerical computations. We saw that theasymptotic prediction provided an accurate approximation for the wave behavior in the far field. However,the asymptotic approximation of the far-field wave amplitude becomes less accurate for progressively smallervalues of the precompression parameter δ . It is possible that this is a consequence of the inadequacy of exam-ining a weakly nonlinear regime of the woodpile chain as nonlinear effects become more significant, leading toerrors in the description of exponentially small effects. A study of woodpile chains with zero precompressionwould reveal whether it is possible to analyze the behavior of the oscillations in the strongly nonlinear regimeusing exponential asymptotic methods. We discuss the possibility of such a study at the end of this section.We also used exponential asymptotic analysis to study nanoptera in a diatomic Hertzian chain in whichthe mass ratio between light and heavy particles is small. We found that there are two Stokes curves in theanalytically-continued traveling-wave solution. These produce exponentially small oscillations — which haveidentical amplitude but different phases — behind the wave front. The interaction between these oscillationsis particularly important, as there exists a set of values of the mass-ratio parameter η at which the oscillationsinterfere destructively to produce a localized solitary wave.Our comparison of the asymptotic results and numerical computations in diatomic Hertzian chains showsthat our exponential asymptotic analysis accurately predicts the amplitude of the far-field waves for many alues of η . However, our computations also reveal that our asymptotic analysis is unable to detect the wavecancellation that occurs for several values of η . Furthermore, the predicted mass ratio values that produce wavecancellation are less accurate than those in analyses of other particle chain behavior (such as the exponentialasymptotic calculations from [55, 56]), particularly for small values of the long-wavelength parameter (cid:15) . Itis possible that this inaccuracy arises from our use of a long-wavelength approximation for the leading-ordersolution, which may have limited accuracy away from the real axis. This is an important question, as theStokes phenomenon depends on the behavior of singularities in the analytically-continued solution.In addition to deriving an explanation for the above undetected features of nanoptera in diatomic Hertzianchains, our investigation opens several interesting questions. In the present study, we used a similar asymptoticapproach as the one in [55,56] and evaluated our asymptotic approximations by comparing them to numericalcomputations. It would be interesting to develop a rigorous existence proof of nanopteron solutions in bothdiatomic Hertzian chains and woodpile chains. It may be possible to adapt the Beale ansatz method of [26,38]for this purpose.As we discussed in Section 1, one-sided nanopteron solutions in both woodpile and diatomic Hertzianchains are metastable, because the oscillations slowly draw energy away from the wave front. This indicatesthat the wave cannot persist indefinitely; instead, it must eventually decay. This property is not present inthe leading-order approximation, nor is it apparent in the form of the exponentially small oscillations in theanalysis in either the present study or in [55, 56]. The long-time decay must arise from interactions betweenthe leading-order traveling wave and the exponentially small oscillations. It would be useful to determine (1)how the leading-order wave and the exponentially small oscillations interact and (2) whether one can detectsuch interactions by using an appropriate long-time rescaling or by calculating higher-order corrections to theoscillatory wave train.In our analysis, we considered only a weakly nonlinear approximation to the leading-order solitary wave forboth woodpile and diatomic Hertzian chains. In principle, one can apply the exponential asymptotic methodthat was developed in [17] with few additional complications to strongly nonlinear systems, such as a Hertziansystem without precompression. However, there are two significant challenges that need to be overcome toundertake such an analysis. The first challenge is that one no longer has access to a closed-form leading-ordersolution to analytically continue in a straightforward fashion. We hope that it may be possible to use theleading-order approximation that was derived in [77] using a rational-function approximation, but it is notapparent whether this approximation is valid if it is analytically continued away from the real axis. This isan important point to clarify, as the Stokes phenomenon is caused by singularities in the complex plane. Thesecond challenge relates to the form of the potential in (2). This potential is not smooth at any point atwhich the argument in the brackets is 0. In the present study, the precompression parameter δ is alwaysat least as large as the relative displacement of two adjacent particles, so the argument in the brackets isalways nonnegative and the potential in (2) is smooth. Without precompression, however, the argument inthe brackets can switch between positive and negative signs, and the potential is not smooth at any point atwhich the argument in the brackets is 0. An important direction for future studies is to explore the effects ofsuch a non-smooth potential on exponential asymptotic analysis. If these challenges are overcome, it wouldbe useful to compare the ensuing results with the analysis of the nonlinear Hertzian woodpile systems thatwere examined in Xu et al. [89], who used fixed-point Fourier-transform methods to identify anti-resonanceconditions. In particular, it would be interesting to determine whether one can derive these conditions usingexponential asymptotic methods and whether there are connections between such an approach and the one inXu et al. A Determining the Prefactor Constants
We determine the prefactor constants for both the woodpile chain and diatomic Hertzian chains. Assumingthat the power series diverges in a factorial-over-power fashion, it must cease to be asymptotic in a narrowregion in the neighborhood of the singularities in the complex plane. Inside this region, the earlier terms inthe series are not larger asymptotically than later series terms in the limit that η →
0. Describing the solution ehavior near one of these singularities requires obtaining a local expansion of the solution near the singularpoint. Once we have obtained this local expansion, we obtain the prefactor constants in the late-order termsusing asymptotic matching to ensure that the local expansions are consistent with the power-series behavior. A.1 Woodpile chain
To determine the value of Λ + , we match the late-order expansion in the outer region with the local solutionin an inner region near the singularity at ξ = ξ by using Van Dyke’s matching principle.As ξ → ξ , we find that u ( ξ ) ∼ − α − δξ − ξ + O ( ξ − ξ ) , u ( ξ + 1) ∼ α − δ(cid:15) coth( (cid:15) ) + O ( ξ − ξ ) , (89) v ( ξ ) ∼ − α − δξ − ξ + O ( ξ − ξ ) , u ( ξ − ∼ − α − δ(cid:15) coth( (cid:15) ) + O ( ξ − ξ ) . (90)To locate the relevant inner region, we need to determine where the validity of the late-order term ansatzbreaks down. From the form of the late-order ansatz (25), we see that this occurs for η χ − = O (1) as η →
0. That is, it occurs for η ( ξ − ξ ) − = O (1); this corresponds to the inner scaling ξ − ξ = ηξ . Fromasymptotic balancing, the appropriate rescaled inner variables are u ( ξ ) = − α − δηξ + ˆ u ( ξ ) , u ( ξ +1) = ˆ u ( ξ + η − ) , u ( ξ −
1) = ˆ u ( ξ − η − ) , v ( ξ ) = − α − δηξ + ˆ v ( ξ ) η . (91)Retaining the leading-order terms as η →
0, the rescaled inner equation gives − α − δξ + d ˆ v ( ξ )d ξ = − kc (cid:15) ˆ v ( ξ ) . (92)We express ˆ v in terms of a power seriesˆ v ( ξ ) ∼ ∞ (cid:88) j =1 v n ξ j +1 as ξ → , (93)and we note that we include the leading-order singularity as part of the rescaling process (91). This yields − α − δξ + ∞ (cid:88) j =1 (2 j + 1)(2 j + 2) v j ξ j +3 = − kc (cid:15) ∞ (cid:88) j =1 v j ξ j +1 . (94)By matching orders of ξ , we obtain the recurrence relation v = 2 δc (cid:15) / ( k ( α − , (2 j + 2)(2 j + 1) c (cid:15) v j = − kv j +1 . (95)Solving the recurrence relation (95) gives v j = 1 α − δ ( − j +1 (cid:18) c (cid:15) k (cid:19) j Γ(2 j + 1) . (96)By comparing the series expression (96) with the inner limit of the late-order ansatz, we obtainΛ + = lim j →∞ v j (i √ k/c (cid:15) ) j +1 Γ(2 j + 1) = − α − δ √ kc (cid:15) . (97) A.2 Diatomic Hertzian Chain
To determine Λ − , we match the late-order term expansion in the outer region with a local expansion inan inner region near each relevant singularity. We perform this analysis in the neighborhood of ξ = ξ − .From (62), we see that the factorial-over-power ansatz breaks down when η χ − = O (1) as η →
0. Weintroduce an inner scaling ξ − ξ − = η / ¯ ξ , which gives the rescaled inner variables u ( ξ + 1) = − δη / ¯ ξ + ˆ u (¯ ξ + η − / ) , u ( ξ −
1) = ˆ u (¯ ξ − η − / ) , v ( ξ ) = − δη / ¯ ξ + ˆ v (¯ ξ ) η / . (98) (cid:15) Λ approx /δ / j
40 38 . . . Figure 9: Determining Λ approx using (101) for progressively larger values of j . The dashed line indicates the value to which Λ approx appears to be converging as we increase j . Retaining the leading-order terms as η → − δ ¯ ξ + d ˆ v (¯ ξ )d¯ ξ = 1 c (cid:15) (cid:34)(cid:18) δ ¯ ξ − ˆ v ( ξ ) (cid:19) / − (cid:18) δ ¯ ξ + ˆ v ( ξ ) (cid:19) / (cid:35) . (99)We express ˆ v in terms of the power seriesˆ v (¯ ξ ) ∼ ∞ (cid:88) j =1 a j ¯ ξ j +1 as ¯ ξ → . (100)Matching the inner expansion (100) with the outer ansatz (72) yieldsΛ − = lim j →∞ a j i j +2 j +4 δ j/ / j +1 / c j +2 (cid:15) Γ(2 j + 3 / . (101)We numerically compute values for a j by substituting the series expression (100) into (99) and solving arecurrence relation for a j in terms of the previous series coefficients. We obtain an approximation Λ approx forΛ − by evaluating the right-hand side of (101) for large values of j . When we take j to be sufficiently large,we obtain an accurate approximation of the exact value of Λ − .In Figure 9, we show the behavior of Λ approx as we increase j , and we can see that our approximation isconverging. For sufficiently large values of j , we compute thatΛ − = 38 . δ / /c (cid:15) . (102)Following the same approach as above yieldsΛ = iΛ − . (103) References [1] G. L. Alfimov, A. S. Korobeinikov, C. J. Lustri, and D. E. Pelinovsky. Standing lattice solitons in thediscrete NLS equation with saturation.
Nonlinearity , 32(9):3445–3484, 2019.[2] M. P. Allen and D. J. Tildesley.
Computer Simulation of Liquids . Clarendon Press, Oxford, UK, 1987.[3] E. ´Avalos and S. Sen. How solitary waves collide in discrete granular alignments.
Phys. Rev. E , 79:046607,2009.[4] E. ´Avalos and S. Sen. Granular chain between asymmetric boundaries and the quasiequilibrium state.
Phys. Rev. E , 89:053202, 2014.[5] E. ´Avalos, D. Sun, R. L. Doney, and S. Sen. Sustained strong fluctuations in a nonlinear chain at acousticvacuum: Beyond equilibrium.
Phys. Rev. E , 84:046610, 2011.
6] E. S. Benilov, R. Grimshaw, and E. P. Kuznetsova. The generation of radiating waves in a singularly-perturbed Korteweg–de Vries equation.
Physica D , 69(3–4):270–278, 1993.[7] M. V. Berry. Stokes’ phenomenon; smoothing a victorian discontinuity.
Pub. Math. de L’IH ´ES , 68:211–221, 1988.[8] M. V. Berry. Uniform asymptotic smoothing of Stokes’s discontinuties.
Proc. Roy. Soc. Lond. A ,422(1862):7–21, 1989.[9] M. V. Berry. Asymptotics, superasymptotics, hyperasymptotics. In H. Segur, S. Tanveer, and H. Levine,editors,
Asymptotics Beyond All Orders , pages 1–14. Plenum Publishing Corporation, Amsterdam, TheNetherlands, 1991.[10] M. V. Berry and C. J. Howls. Hyperasymptotics.
Proc. Roy. Soc. Lond. A , 430(1880):653–668, 1990.[11] N. Boechler, G. Theocharis, S. Job, P. G. Kevrekidis, M. A. Porter, and C. Daraio. Discrete breathers inone-dimensional diatomic granular crystals.
Phys. Rev. Lett. , 104:244302, 2010.[12] J. P. Boyd. A numerical calculation of a weakly non-local solitary wave: The φ breather. Nonlinearity ,3(1):177–195, 1990.[13] J. P. Boyd. Weakly non-local solitons for capillary-gravity waves: Fifth-degree Korteweg–de Vries equa-tion.
Physica D , 48:129–146, 1991.[14] J. P. Boyd.
Weakly Nonlocal Solitary Waves and Beyond-All-Orders Asymptotics: Generalized Solitonsand Hyperasymptotic Perturbation Theory , volume 442 of
Mathematics and Its Applications . KluwerPublishers, Amsterdam, The Netherlands, 1998.[15] J. P. Boyd. The devil’s invention: Asymptotic, superasymptotic and hyperasymptotic series.
Acta Appl.Math. , 56(1):1–98, 1999.[16] J. P. Boyd. Hyperasymptotics and the linear boundary layer problem: Why asymptotic series diverge.
SIAM Rev. , 47(3):553–575, 2005.[17] S. J. Chapman, J. R. King, and K. L. Adams. Exponential asymptotics and Stokes lines in nonlinearordinary differential equations.
Proc. Roy. Soc. Lond. A , 454(1978):2733–2755, 1998.[18] C. Chong, M. A. Porter, P. G. Kevrekidis, and C. Daraio. Nonlinear coherent structures in granularcrystals.
J. Phys. Cond. Matt. , 29(41):413003, 2017.[19] C. Coste, E. Falcon, and S. Fauve. Solitary waves in a chain of beads under Hertz contact.
Phys. Rev.E , 56:6104–6117, 1997.[20] C. Daraio, V. F. Nesterenko, E. B. Herbold, and S. Jin. Energy trapping and shock disintegration in acomposite granular medium.
Phys. Rev. Lett. , 96:058002, 2006.[21] C. Daraio, V. F. Nesterenko, E. B. Herbold, and S. Jin. Tunability of solitary wave properties in one-dimensional strongly nonlinear phononic crystals.
Phys. Rev. E , 73:026610, 2006.[22] G. Deng, G. Biondini, and S. Sen. Interactions of solitary waves in integrable and nonintegrable lattices.
Chaos , 30(4):043101, 2020.[23] G. Deng, G. Biondini, and S. Sen. On the generation and propagation of solitary waves in integrable andnonintegrable nonlinear lattices.
Eur. Phys. J. Plus , 135:598, 2020.[24] R. B. Dingle.
Asymptotic Expansions: Their Derivation and Interpretation . Academic Press, New York,NY, USA, 1973.[25] E. Hasco¨et and H. J. Herrmann. Shocks in non-loaded bead chains with impurities.
Eur. Phys. J. B ,14(1):183–190, 2000.[26] T. E. Faver. Nanopteron-stegoton traveling waves in spring dimer Fermi–Pasta–Ulam–Tsingou lattices.
Q. Appl. Math. , 78:363–429, 2020.[27] T. E. Faver and H. J. Hupkes. Micropterons, nanopterons and solitary wave solutions to the diatomicFermi–Pasta–Ulam–Tsingou problem. arXiv preprint arXiv:2010.02107 , 2020.
28] A. Feigel, M. Veinger, B. Sfez, A. Arsh, M. Klebanov, and V. Lyubin. Three-dimensional simple cubicwoodpile photonic crystals made from chalcogenide glasses.
Appl. Phys. Lett. , 83(22):4480–4482, 2003.[29] G. Friesecke and R. L. Pego. Solitary waves on FPU lattices: II. Linear implies nonlinear stability.
Nonlinearity , 15(4):1343–1359, 2002.[30] G. Friesecke and J. A. D. Wattis. Existence theorem for solitary waves on lattices.
Commun. Math.Phys. , 161(2):391–418, 1994.[31] N. Giardetti, A. Shapiro, S. Windle, and J. D. Wright. Metastability of solitary waves in diatomic FPUTlattices.
Math. in Eng. , 1:419–433, 2019.[32] R. Grimshaw and N. Joshi. Weakly nonlocal solitary waves in a singularly perturbed Korteweg–de Vriesequation.
SIAM J. Appl. Math. , 55(1):124–135, 1995.[33] U. Harbola, A. Rosas, A. H. Romero, and K. Lindenberg. Pulse propagation in randomly decoratedchains.
Phys. Rev. E , 82:011306, 2010.[34] E. B. Herbold, J. Kim, V. F. Nesterenko, S. Y. Wang, and C. Daraio. Pulse propagation in a linear andnonlinear diatomic periodic chain: Effects of acoustic frequency band-gap.
Acta Mech. , 205:85–103, 2009.[35] H. Hertz. ¨ueber die ber¨uhrung fester elastischer k¨orper.
J. Reine Angew. Math. , 92:156–171, 1881.[36] E. J. Hinch.
Perturbation Methods . Cambridge Texts in Applied Mathematics. Cambridge UniversityPress, Cambridge, UK, 1991.[37] E. J. Hinch and S. Saint-Jean. The fragmentation of a line of balls by an impact.
P. R. Soc. A , 455:3201–3220, 1999.[38] A. Hoffman and J. D. Wright. Nanopteron solutions of diatomic Fermi–Pasta–Ulam–Tsingou latticeswith small mass-ratio.
Physica D , 358:33–59, 2017.[39] C. Hoogeboom, Y. Man, N. Boechler, G. Theocharis, P. G. Kevrekidis, I. G. Kevrekidis, and C. Daraio.Hysteresis loops and multi-stability: From periodic orbits to chaotic dynamics (and back) in diatomicgranular crystals.
Europhys. Lett. , 101(4):44003, 2013.[40] J. K. Hunter and J. Scheurle. Existence of perturbed solitary wave solutions to a model equation forwater waves.
Physica D , 32(2):253–268, 1988.[41] G. Iooss and K. Kirchg¨assner. Travelling waves in a chain of coupled nonlinear oscillators.
Commun.Math. Phys. , 211(2):439–464, 2000.[42] K. R. Jayaprakash, Y. Starosvetsky, and A. F. Vakakis. New family of solitary waves in granular dimerchains with no precompression.
Phys. Rev. E , 83(3):036606, 2011.[43] K. R. Jayaprakash, Y. Starosvetsky, A. F. Vakakis, and O. V. Gendelman. Nonlinear resonances leadingto strong pulse attenuation in granular dimer chains.
J. Nonlinear Sci. , 23(3):363–392, 2013.[44] H. Jiang, Y. Wang, M. Zhang, Y. Hu, D. Lan, Y. Zhang, and B. Wei. Locally resonant phononic woodpile:A wide band anomalous underwater acoustic absorbing material.
Appl. Phys. Lett. , 95(10):104101, 2009.[45] S. Job, F. Melo, A. Sokolow, and S. Sen. How Hertzian solitary waves interact with boundaries in a 1Dgranular medium.
Phys. Rev. Lett. , 94:178002, 2005.[46] E. Kim, R. Chaunsali, H. Xu, J. Jaworski, J. Yang, P. G. Kevrekidis, and A. F. Vakakis. Nonlinearlow-to-high-frequency energy cascades in diatomic granular crystals.
Phys. Rev. E , 92:062201, 2015.[47] E. Kim, F. Li, C. Chong, G. Theocharis, J. Yang, and P. G. Kevrekidis. Highly nonlinear wave propagationin elastic woodpile periodic structures.
Phys. Rev. Lett. , 114(11):118002, 2015.[48] E. Kim, A. J. Martinez, S. E. Phenisee, P.G. Kevrekidis, M. A. Porter, and J. Yang. Direct measurementof superdiffusive energy transport in disordered granular chain.
Nat. Commun. , 9:640, 2018.[49] E. Kim and J. Yang. Wave propagation in single column woodpile phononic crystals: Formation oftunable band gaps.
J. Mech. Phys. Solids , 71:33–45, 2014.
50] D. J. Korteweg and G. de Vries. XLI. On the change of form of long waves advancing in a rectangularcanal, and on a new type of long stationary waves.
Lond. Edinb. Dubl. Phil. Mag. , 39(240):422–443, 1895.[51] A. N. Lazaridi and V. F. Nesterenko. Observation of a new type of solitary waves in a one-dimensionalgranular medium.
J. Appl. Mech. Tech. Phy. , 26:405–408, 1985.[52] H. Liu, J. Yao, D. Xu, and P. Wang. Characteristics of photonic band gaps in woodpile three-dimensionalterahertz photonic crystals.
Opt. Express , 15(2):695–703, 2007.[53] L. Liu, G. James, P. Kevrekidis, and A. Vainchtein. Breathers in a locally resonant granular chain withprecompression.
Physica D , 331:27–47, 2016.[54] L. Liu, G. James, P. Kevrekidis, and A. Vainchtein. Strongly nonlinear waves in locally resonant granularchains.
Nonlinearity , 29(11):3496–3527, 2016.[55] C. J. Lustri. Nanoptera and Stokes curves in the 2-periodic Fermi–Pasta–Ulam–Tsingou equation.
PhysicaD , 402:132239, 2020.[56] C. J. Lustri and M. A. Porter. Nanoptera in a period-2 Toda chain.
SIAM J. Appl. Dyn. Syst. , 17(2):1182–1212, 2018.[57] F. S. Manciu and S. Sen. Secondary solitary wave formation in systems with generalized Hertz interactions.
Phys. Rev. E , 66:016616, 2002.[58] M. Manciu, S. Sen, and A. J. Hurd. Impulse propagation in dissipative and disordered chains withpower-law repulsive potentials.
Physica D , 157(3):226–240, 2001.[59] M. Manjunath, A. P. Awasthi, and P. H. Geubelle. Wave propagation in random granular chains.
Phys.Rev. E , 85:031308, 2012.[60] M. Manjunath, A. P. Awasthi, and P. H. Geubelle. Family of plane solitary waves in dimer granularcrystals.
Phys. Rev. E , 90:032209, 2014.[61] A. J. Martinez, P. G. Kevrekidis, and M. A. Porter. Superdiffusive transport and energy localization indisordered granular crystals.
Phys. Rev. E , 93:022902, 2016.[62] A. J. Martinez, M. A. Porter, and P. G. Kevrekidis. Quasiperiodic granular chains and hofstadterbutterflies.
Philos. T. Roy. Soc. A , 376:20170139, 2018.[63] A. J. Martinez, H. Yasuda, E. Kim, P. G. Kevrekidis, M. A. Porter, and J. Yang. Scattering of waves byimpurities in precompressed granular chains.
Phys. Rev. E , 93:052224, 2016.[64] A. Molinari and C. Daraio. Stationary shocks in periodic highly nonlinear granular chains.
Phys. Rev.E , 80:056602, 2009.[65] V. F. Nesterenko. Propagation of nonlinear compression pulses in granular media.
J. Appl. Mech. Tech.Phy. , 24:733–743, 1983.[66] V. F. Nesterenko.
Dynamics of Heterogeneous Materials . Springer-Verlag, Heidelberg, Germany, 2001.[67] V. F. Nesterenko, C. Daraio, E. B. Herbold, and S. Jin. Anomalous wave reflection at the interface oftwo strongly nonlinear granular media.
Phys. Rev. Lett. , 95(15):158702, 2005.[68] Y. Okada, S. Watanabe, and H. Tanaca. Solitary wave in periodic nonlinear lattice.
J. Phys. Soc. Jpn. ,59(8):2647–2658, 1990.[69] A. B. Olde Daalhuis, S. J. Chapman, J. R. King, J. R. Ockendon, and R. H. Tew. Stokes phenomenonand matched asymptotic expansions.
SIAM J. Appl. Math. , 55(6):1469–1483, 1995.[70] L. Ponson, N. Boechler, Y. M. Lai, M. A. Porter, P. G. Kevrekidis, and C. Daraio. Nonlinear waves indisordered diatomic granular chains.
Phys. Rev. E , 82:021301, 2010.[71] M. A. Porter, C. Daraio, E. B. Herbold, I. Szelengowicz, and P. G. Kevrekidis. Highly nonlinear solitarywaves in periodic dimer granular chains.
Phys. Rev. E , 77:015601, 2008.[72] M. A. Porter, C. Daraio, I. Szelengowicz, E. B. Herbold, and P. G. Kevrekidis. Highly nonlinear solitarywaves in heterogeneous periodic granular media.
Physica D , 238(6):666–676, 2009.
73] M. A. Porter, P.G. Kevrekidis, and C. Daraio. Granular crystals: Nonlinear dynamics meets materialsengineering.
Phys. Today , 68(11):44–50, 2015.[74] R. Potekin, K. R. Jayaprakash, D. M. McFarland, K. Remick, L. A. Bergman, and A. F. Vakakis.Experimental study of strongly nonlinear resonances and anti-resonances in granular dimer chains.
Exp.Mech. , 53:861–870, 2013.[75] M. Przedborski, T. A. Harroun, and S. Sen. Granular chains with soft boundaries: Slowing the transitionto quasiequilibrium.
Phys. Rev. E , 91:042207, 2015.[76] M. Przedborski, S. Sen, and T. A. Harroun. The equilibrium phase in heterogeneous Hertzian chains.
J.Stat. Mech. — Theory E. , 2017(12):123204, 2017.[77] S. Sen, J. Hong, J. Bang, E. ´Avalos, and R. Doney. Solitary waves in the granular chain.
Phys. Rep. ,462(2):21–66, 2008.[78] S. Sen, T. R. Krishna Mohan, and J. M. M. Pfannes. The quasi-equilibrium phase in nonlinear 1d sys-tems.
Physica A , 342(1):336–343, 2004. Proceedings of the VIII Latin American Workshop on NonlinearPhenomena.[79] S. Sen, M. Manciu, and J. D. Wright. Solitonlike pulses in perturbed and driven Hertzian chains andtheir possible applications in detecting buried impurities.
Phys. Rev. E , 57:2386–2397, 1998.[80] Y. Shen, P. G. Kevrekidis, S. Sen, and A. Hoffman. Characterizing traveling-wave collisions in granularchains starting from integrable limits: The case of the Korteweg–de Vries equation and the Toda lattice.
Phys. Rev. E , 90:022905, 2014.[81] G. G. Stokes. On the discontinuity of arbitrary constants which appear in divergent developments.
Trans.Cam. Phil. Soc. , 10:106–128, 1864.[82] Y. Tabata. Stable solitary wave in diatomic Toda lattice.
J. Phys. Soc. Jpn. , 65(12):3689–3691, 1996.[83] G. Theocharis, N. Boechler, P. G. Kevrekidis, S. Job, M. A. Porter, and C. Daraio. Intrinsic energylocalization through discrete gap breathers in one-dimensional diatomic granular crystals.
Phys. Rev. E ,82:056604, 2010.[84] A. Vainchtein, Y. Starosvetsky, J. D. Wright, and R. Perline. Solitary waves in diatomic chains.
Phys.Rev. E , 93(4):042210, 2016.[85] L. Vergara. Scattering of solitary waves from interfaces in granular media.
Phys. Rev. Lett. , 95:108002,2005.[86] L. Vergara. Delayed scattering of solitary waves from interfaces in a granular container.
Phys. Rev. E ,73:066623, 2006.[87] L. Verlet. Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jonesmolecules.
Phys. Rev. , 159:98–103, 1967.[88] L. Y. Wu and L. W. Chen. Acoustic band gaps of the woodpile sonic crystal with the simple cubic lattice.
J. Phys. D Appl. Phys. , 44(4):045402, 2011.[89] H. Xu, P. G. Kevrekidis, and A. Stefanov. Traveling waves and their tails in locally resonant granularsystems.
J. Phys. A Math. Theor. , 48(19):195204, 2015., 48(19):195204, 2015.