# Nanoptera in weakly nonlinear woodpile chains and diatomic granular chains

G. Deng1, C. J. Lustri1111Corresponding Author. Electronic address:
<EMAIL_ADDRESS>and Mason A. Porter2

(1Department of Mathematics and Statistics, 12 Wally’s Walk, Macquarie
University, New South Wales 2109, Australia
2Department of Mathematics, 520 Portola Plaza, University of California, Los
Angeles, California 90024, USA )

###### Abstract

We study “nanoptera”, which are non-localized traveling waves with
exponentially small but non-decaying oscillations, in two singularly-perturbed
Hertzian chains with precompression. These two systems are woodpile chains
(which we model as systems of Hertzian particles and springs) and diatomic
Hertzian chains with alternating masses. We demonstrate that nanoptera arise
from Stokes phenomena and appear as special curves, which are called “Stokes
curves”, are crossed in the complex plane. We use techniques from exponential
asymptotics to obtain approximations of the oscillation amplitudes. Our
analysis demonstrates that traveling waves in a singularly perturbed woodpile
chain have a single Stokes curve, across which oscillations appear. Comparing
these asymptotic approximations with numerical simulations reveals that they
accurately describe the non-decaying oscillatory behavior in a woodpile chain.
We perform a similar analysis of a diatomic Hertzian chain and show that its
nanopteron solution has two distinct exponentially small oscillatory
contributions. We demonstrate that there exists a set of mass ratios for which
these two contributions cancel to produce localized solitary waves. This
result builds on prior experimental and numerical observations that there
exist mass ratios that support localized solitary waves in diatomic Hertzian
chains without precompression. Comparing asymptotic and numerical results in a
diatomic Hertzian chain with precompression reveals that our exponential
asymptotic approach accurately predicts the oscillation amplitude for a wide
range of system parameters but that it fails to identify several values of the
mass ratio that correspond to localized solitary-wave solutions.

##  1 Introduction

###  1.1 Particle Chains

The behavior of a particle chain under compression depends on the character of
the interactions between particles. A key example is the Hertzian interaction
[37], which describes the relationship between the repelling force from the
contact area of adjacent frictionless spherical particles and the compression
between them. Hertzian interactions were first used to describe static chain
configurations. Nesterenko [66] demonstrated that it is also reasonable to
model dynamical processes such as wave propagation in particle chains using
Hertzian interactions between adjacent particles if certain constraints on the
particle contact areas, the stresses at the contact points, and the
characteristic time scales of the dynamics are satisfied. These constraints
were validated in experimental studies [53, 19]. Chains of particles with
Hertzian interactions, which are typically called “Hertzian chains”, have been
the subject of numerous investigations in the last few decades [71, 83, 18]
both because they possess rich dynamics and because of their potential uses
for practical engineering applications [79].

A general equation that governs the motion of particles in a chain is the
following system of differential–difference equations:

$m(n)\ddot{x}(n,t)=\phi^{\prime}(x(n+1,t)-x(n,t))-\phi^{\prime}(x(n,t)-x(n-1,t))\,,$
(1)

where $n\in\mathbb{Z}$, the quantity $m(n)$ is the mass of the $n$-th
particle, $x(n,t)$ is the position of the $n$-th particle at time $t$, a dot
denotes differentiation with respect to time, a prime denotes differentiation
with respect to space, and the interaction potential between adjacent
particles is

$\displaystyle\phi(r)=\begin{cases}c(\delta_{0}-r)^{\alpha+1}\,,&r\leq\delta_{0}\cr
0\,,&r>\delta_{0}\end{cases}\,,~{}\quad c=\mathrm{constant}\,,$ (2)

where $\alpha>1$ and $\delta_{0}$ is the equilibrium overlap of adjacent
particles that arise from the precompression that is induced by an external
force. In Fig. 1, we illustrate a particle chain with precompression. The
power-law interaction potential (2) is $0$ when particles are not in contact,
and it cannot take negative values. For algebraic convenience, we set
$c=1/(\alpha+1)$ in all of our examples. The exponent $\alpha$ in the power-
law interaction potential (2) depends on the contact geometry between adjacent
particles [88]. As was illustrated in [28, 90], $\alpha$ is tunable
experimentally. The choice of $\alpha=3/2$ gives the widely-studied Hertzian
interaction potential, but it is also relevant to consider other values of the
exponent $\alpha$ [18]. For example, Nesterenko [69] investigated transverse
vibrations of unstressed linear elastic fibers, for which $\alpha=3$. Sen and
Manciu [85] and Avalos and Sen [3] investigated chains with power-law
interaction potentials for several different values of $\alpha$. The
transition from interactions with $\alpha=2$ to Hertzian interactions was
investigated in [32].

The existence of solitary waves was first reported for Hertzian chains with
precompression in [66]. It was subsequently shown that particle chains with
other interaction exponents $\alpha>1$ also support the propagation of
solitary waves [68, 67]. Friesecke and Wattis [30] showed that solitary-wave
solutions can occur in nonlinear lattices with a broad class of interaction
potentials that are known as “superquadratic potentials”222A interaction
potential $\phi(r)$ is “superquadratic” when $\phi(r)/r^{2}$ strictly
increases with $|r|$ for all $r\in\Lambda$ for either $\Lambda=(-\infty,0)$ or
$\Lambda=(0,\infty)$., which includes power-law interactions with $\alpha>1$.
See [71, 83, 18] for reviews of the properties of solitary waves in granular
chains. Particle chains with $\alpha<1$ do not support solitary waves; see
[70, 99, 35] for further discussion of such chains.

Since the first prediction of solitary waves in Hertzian chains [66], there
have been many theoretical, numerical, and experimental studies of the
properties of solitary waves in Hertzian chains [71, 83, 79, 18, 72]. These
studies have concerned a wide variety of topics, including the generation [23,
53, 19, 21, 39], propagation [23, 53, 19, 21, 39], interaction [58, 47, 3, 5,
4, 22], and long-time dynamics [5, 4, 84, 81, 82] of solitary waves in
monoatomic Hertzian chains, in which every particle is the same.

Figure 1: Schematic illustration of a particle chain with precompression. The
chain consists of $N$ identical, aligned spheres with radius $R$. The
interaction potential between adjacent spheres is governed by (2). Applying a
compressive force $F$ at both ends of the chain compresses each sphere at the
contact point by a uniform distance $\delta_{0}$.

The system (1)–(2) that describes the dynamics of a particle chain is both
nonlinear and nonintegrable, so typically one cannot study its solutions using
exact, fully analytical methods. However, under the assumption of small
deformations, such that $|x(n-1,t)-x(n,t)|/\delta_{0}\ll 1$, one can expand
the right-hand side of (1) as a Taylor series in the normalized deformation
size [66]. Neglecting terms of order $O((|x(n-1,t)-x(n,t)|/\delta_{0})^{4})$
produces the frequently studied Fermi–Pasta–Ulam–Tsingou (FPUT) chain. In the
long-wavelength limit, in which the characteristic spatial size of a solution
is much larger than the lattice spacing, one can approximate the lattice
system (1)–(2) by a continuous system. In the limit of both small deformations
and long wavelengths in comparison to particle size, the monoatomic Hertzian
system with precompression reduces to the Kortweg–de Vries (KdV) equation
[71]. In this long-wavelength limit, one can approximate the solitary-wave
behavior in a Hertzian chain as a soliton solution of the KdV equation [52].
Building on this idea, the interaction of solitary waves in a Hertzian chain
in the long-wavelength limit was described in [87] using a two-soliton
solution of the KdV equation.

(a) A physical woodpile system (b) A woodpile-chain model (c) A diatomic chain
of particles

Figure 2: The chains of particles that we examine in the present study. The
schematic illustration in (a) shows the physical configuration of orthogonal
rods that is known as a “woodpile chain”. The schematic in (b) shows an
idealized mathematical model of the physical configuration in (a). This model
consists of heavy spherical particles in physical contact via interactions
governed by (2) and light spherical particles (which are sometimes called
“resonators”) that are attached to each heavy particle by a spring. The
schematic illustration in (c) shows a diatomic chain of alternating heavy and
light spheres in physical contact via Hertzian interactions. We analyze the
models in (b) and (c).

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 homogeneous
Hertzian system, such as by incorporating different types of particle
heterogeneities. Such modifications can produce behavior that is very
different from that of homogeneous Hertzian chains. For example, researchers
have explored the behavior of Hertzian chains with impurities [86, 25, 64],
disordered particle arrangements [50, 62, 59, 34, 60, 66], quasiperiodicity
[63], and particles that are composed of segments with different materials
(i.e., so-called “compound chains”) [94, 95, 73, 20]. It has been demonstrated
that impurities can cause solitary waves to scatter [86, 25, 64]. Transport
and localization of energy in Hertzian systems with disorder and
quasiperiodicity were studied in [50, 62, 63]. In Hertzian chains with
randomly arranged masses, solitary waves can delocalize, such that the energy
that is carried by a solitary wave spreads among many particles in a chain
[59, 34, 60]. The scattering of solitary waves at an interface of a compound
Hertzian system was studied in [94, 73]. Other phenomena that have been
studied in numerical and experimental investigations of heterogeneous Hertzian
chains include energy trapping [95, 20], shock disintegration [20], and the
generation of secondary solitary waves [94].

It is also interesting to consider other types of heterogeneous chains,
including (1) particle–rod systems that model so-called “woodpile chains” [18,
49, 98, 55, 54, 46, 97, 51] and (2) diatomic chains [18, 92, 11, 65, 76, 77,
78, 41, 36, 44, 45, 48, 80]. It is feasible to study both woodpile chains and
diatomic chains in laboratory experiments.

A woodpile chain consists of orthogonally stacked slender rigid cylinders with
mass $m_{1}$ [see Fig. 2(a)]. The interaction force along the direction of the
stack is determined by (2), and one can model the elastic deformation along
the direction that is perpendicular to the stack direction using internal
“resonators”, where the mass and coupling constant of these resonators come
from the material and shape of the cylinders. When considering a stack of
identical cylinders, each resonator has the same mass $m_{2}$ and elastic
constant $k$. We model such a woodpile chain as a monoatomic granular chain in
which each particle is connected to an external particle by a rod, which we
treat as a linear spring. Each particle in the chain has mass $m_{1}$, each
external particle has mass $m_{2}$, and each linear spring has spring constant
$k$. We illustrate this woodpile chain model in Fig. 2(b). For the rest of the
present study, we use the term “woodpile chain” to refer to this idealized
system. Much of the existing work on woodpile chains has concentrated on their
linear elastic responses [46, 97, 51]. Some recent work studied the behavior
of traveling waves in woodpile chains with no precompression [49, 98] and
hence in a strongly nonlinear regime. Other studies have investigated the
existence of discrete breathers in woodpile chains both without precompression
[55] and with precompression [54].

A diatomic chain consists of particles with some characteristic, such as
particle mass, that alternates in adjacent particles. We consider chains [see
Fig. 2(c)] in which even particles have mass $m_{1}$ and odd particles have
mass $m_{2}$. Diatomic Hertzian chains have rich dynamics [18]. They can have
both discrete breather solutions [92, 11] and stationary shock waves [65].
Previous studies have investigated the propagation and scattering of nonlinear
waves in both ordered and disordered diatomic Hertzian chains [76, 77, 78, 36]
and chaotic dynamics in diatomic Hertzian chains that are damped and driven
[41].

###  1.2 Nanoptera

Typically, traveling-wave solutions in woodpile and diatomic Hertzian systems
are not truly localized [44, 49, 98]; instead, they also include a small non-
localized oscillatory tail 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 one or both
sides of the central wave333Several studies have used 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.. These oscillations
typically have an exponentially small amplitude with respect to a singular
perturbation parameter. (In the present study, this parameter is the mass
ratio between one type of particle and the other type of particle.) We show
examples of nanoptera in Fig. 3. A genuine solitary wave is exponentially
localized, but a nanopteron is localized only up to algebraic orders in the
small parameter.

(a) Solitary wave (b) One-sided nanopteron (c) Two-sided nanopteron

Figure 3: Comparison of the profiles of (a) a standard solitary wave, (b) a
one-sided nanopteron, and (c) a two-sided nanopteron. The waves propagate at
speeds of $c_{0}$, $c_{1}$, and $c_{2}$ respectively. The solitary wave is
localized spatially, whereas the nanoptera have non-decaying oscillatory tails
on (b) one side or (c) both sides of the wave front. The waves in (a) and (c)
propagate without decaying, but the wave in (b) cannot propagate indefinitely.
In (b), the nanoptera decay very slowly because the one-sided radiation draws
energy from the wave front. Additionally, because this decay occurs over a
long time scale, the one-sided nanopteron in (b) is said to be “metastable”.

Jayaprakash _et al._ [44] showed that for uncompressed diatomic Hertzian
chains non-decaying oscillations are absent at certain values of the system
parameters. These sets of parameters arise from satisfying an “anti-resonance
condition”. When such anti-resonance condition is satisfied, nanoptera in
uncompressed diatomic Hertzian chains become solitary waves that propagate
without attenuation. Other studies that have examined anti-resonance
conditions in diatomic Hertzian chains include [48, 80]. Hertzian chains can
also have nanopteron solutions that travel with particularly strong dispersion
and attenuation. Such behavior occurs if a chain satisfies a “resonance
condition” [45, 48, 80], which one obtains in a similar fashion to an anti-
resonance condition. The study of nanopteron solutions in one-dimensional
Hertzian chains [44] was extended to two-dimensional Hertzian systems by
Manjunath _et al._ [61], who obtained anti-resonance conditions for these
systems and examined the properties of solitary-wave solutions when these
anti-resonance conditions are satisfied. Anti-resonance conditions were also
identified for uncompressed woodpile chains by Xu et al. [98].

Nanoptera have been studied in many other particle chains, such as diatomic
Toda chains [93, 74, 91, 57], diatomic FPUT chains [93, 26, 40, 56], and
chains with an on-site nonlinear potential in which each particle is linearly
coupled to its neighbors [43]. The existence of nanopteron solutions in a
diatomic FPUT chain was proven rigorously in [26, 40], and nanoptera in
periodic Toda chains were studied using a multiple-scale approach in [93].
Recently, [27] studied diatomic FPUT chains using numerical continuation. They
explored the relation between nanopteron solutions for chains in which one
mass is much larger than the other (i.e., the so-called “small-mass-ratio”
regime) and traveling-wave solutions for chains with a mass ratio close to
$1$.

Anti-resonance conditions, at which the small oscillations disappear entirely
and a nanopteron becomes a genuine solitary wave, have been identified in
diatomic Toda chains [93, 91, 57], diatomic FPUT chains [57], woodpile chains
without precompression [98], and a discrete nonlinear Schrödinger equation
[1].

An exponential asymptotic method was applied recently by Lustri _et al._ [57,
56] to study exponentially small oscillations in diatomic Toda and FPUT
chains. They showed that the tailing oscillations in these diatomic systems
are 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. They demonstrated that traveling-wave solutions in these
diatomic chains possess two Stokes curves, which generate two distinct
oscillations with the same amplitude but different phases. For particular
values of the mass ratio, the oscillations are precisely out of phase. The
oscillations thus cancel, and the solution is a localized solitary wave. In
[57, 56], Lustri _et al._ also derived asymptotic anti-resonance conditions
for mass ratios for these systems. In the present study, we use exponential
asymptotic methods to examine nanoptera in singularly perturbed woodpile and
diatomic Hertzian chains under precompression.

As we illustrated in Fig. 3, nanoptera can have oscillations on either one or
both sides of a central traveling wave. The existence of symmetric non-
decaying two-sided nanopteron solutions has been proven for a fifth-order KdV
equation [42]. Subsequently, symmetric nanopteron solutions of this fifth-
order KdV equation were constructed numerically using perturbation series in
Boyd [13]; Boyd then showed that his numerical approach cannot be used to
compute one-sided nanopteron solutions in that system. It was proven in [6]
that one-sided nanopteron solutions for the fifth-order KdV equation cannot
propagate indefinitely without changing form. The one-sided oscillation draws
energy from the central wave, which causes the amplitude of the central core
to decay. The eventual decay of nanopteron solutions of the fifth-order KdV
equation is a result of energy conservation, and one can make a similar
argument against the existence of one-sided nanoptera in any energy-conserving
system. The instability of one-sided nanopteron solutions in the fifth-order
KdV equation was discussed further in [33].

Giardetti _et al._ [31] simulated one-sided nanoptera in a diatomic FPUT model
with small mass ratios for a long enough time to obtain a solution that
appeared to be a stable traveling wave. In their simulations, both the central
solitary wave and the oscillatory tail traveled without any apparent decay.
However, the energy of the system that they studied is conserved.
Consequently, it is impossible for both the central core and the one-sided
oscillation to have a constant amplitude indefinitely because the one-sided
oscillation constantly draws energy from the central core [31, 44, 93].
Therefore, it was conjectured in [31] that the decay of one-sided nanoptera in
diatomic FPUT chains occurs on a time scale that is exponentially large in the
limit that the mass ratio parameter tends to $0$. Solutions that decay on
exponentially slow time scales are sometimes known as “metastable” solutions
(they are also sometimes called “quasistable” solutions), which describe
asymptotic solutions that appear to be stable over a duration that is large in
comparison to any inverse power of a small parameter before eventually
decaying. In experimental settings, nanoptera are typically generated by a
pulse moving through an undisturbed chain. This produces one-sided nanoptera,
with no oscillations in the undisturbed region ahead of the pulse. To the best
of our knowledge, there do not currently exist experimental realizations of
two-sided nanoptera. The nanoptera that we consider are all one-sided, and
therefore metastable; however, as we will note at the end of Section 3.3, it
is straightforward to extend our analysis to the stable case of two-sided
nanoptera, which is of theoretical interest.

In the present study, we compute the behavior of nanopteron solutions in
precompressed woodpile and diatomic Hertzian chains in the asymptotic limit
$m_{2}/m_{1}\to 0$. In this limit, both the woodpile and the diatomic Hertzian
chains are singularly perturbed around a monoatomic Hertzian chain.
Jayaprakash _et al._ [44] investigated this limit for uncompressed diatomic
Hertzian chains, but existing studies of woodpile chains have not focused on
the small-mass-ratio regime.

In the present paper, we apply an exponential asymptotic method that was
developed in [17, 75] and was used in [57, 56] to study traveling waves in
diatomic Toda and FPUT chains with small mass ratios. In a typical asymptotic
power-series analysis, one expands a solution as an algebraic series in some
small parameter. One then computes the terms in the series using a recursion
relation that one obtains by asymptotic matching. This process can never
capture 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 describing asymptotic
behavior on this exponentially small scale. Additionally, it only requires the
direct computation of the leading-order series behavior. We describe this
exponential asymptotic approach in detail in Section 2.

###  1.3 Outline of the Paper

We first investigate traveling-wave solutions in a woodpile chain with
precompression in the limit that the wavelength of the central traveling wave
is large in comparison to the particle size; this is a weakly nonlinear
regime. 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 arises because of the Stokes phenomenon. Using an exponential
asymptotic approach, we obtain simple asymptotic expressions for the
exponentially small, constant-amplitude waves in the wake of the leading-order
solitary wave. We demonstrate that the solution has a single oscillatory
contribution in the wake of the leading wave and thus that it can never be
canceled out. Consequently, the woodpile chains that we examine do not have an
anti-resonance condition. This contrasts with the results of Xu _et al._ [98],
who identified anti-resonance conditions for a woodpile system with no
precompression.

We then study traveling waves in a diatomic Hertzian chain using the same
exponential asymptotic approach. The key difference between diatomic Hertzian
chains and woodpile chains is that traveling waves in diatomic Hertzian chains
have two Stokes curves; this yields solutions with two oscillatory
contributions. For certain mass ratios, the two oscillatory contributions
cancel precisely and the traveling-wave solution becomes a genuine localized
solitary wave.

For both woodpile chains and diatomic Hertzian chains, we use the long-
wavelength approximation from [71] to describe the leading-order solitary wave
behavior (as in [56]). This differs from a diatomic Toda system [57], for
which there exists an analytical expression for the leading-order wave
behavior. Therefore, our analysis of woodpile and diatomic Hertzian chains has
two distinct small parameters, which correspond to the small mass ratio and
the large traveling-wave length scale.

Our paper proceeds as follows. In Section 2, we discuss the exponential
asymptotic method that we employ for our analysis. In Sections 3.1–3.3, we use
this exponential asymptotic method to obtain an asymptotic approximation of
nanopteron solutions in a singularly perturbed woodpile chain. In Section 3.4,
we compare the results of this approximation with numerical computations. We
obtain strong agreement between the asymptotic and numerical results. In
Sections 4.1–4.3, we perform a similar exponential asymptotic analysis on a
singularly perturbed diatomic Hertzian chain. In Section 4.4, we compare our
asymptotic and numerical results for the diatomic Hertzian chain. From this
comparison, we see that the employed exponential asymptotic method is useful
for approximating the solution behavior for a wide range of mass ratios.
However, it fails to detect some important features that arise from the
interference between oscillatory contributions. In Section 5, we conclude and
further discuss our results.

##  2 Exponential Asymptotics

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 $\eta$ to denote the associated small
parameter, which is related to mass ratio and satisfies
$\eta^{2}=m_{2}/m_{1}$. It is impossible to determine the form of these
oscillations using only classical asymptotic power series because the
oscillation amplitude is exponentially small in the asymptotic limit $\eta\to
0$. In particular, the amplitude is smaller than any algebraic power of
$\eta$. 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^{\prime}(x),g^{\prime\prime}(x),\ldots;\eta)=0\,,$ (3)

where $\eta$ is a small parameter. To find an asymptotic solution of (3), we
typically expand a solution about some leading-order approximation as a power
series in the perturbation parameter $\eta$. Once we have obtained this
leading-order solution, we analytically continue it into the complex plane.
Typically, this solution is singular at some set of points, which are the end
points of Stokes curves [89]. Stokes curves play an important role in the
present study: as a Stokes curve is crossed, the amplitude of the
exponentially small contribution undergoes a smooth but rapid change in value.
In the present work, we require the exponentially small contributions to
vanish on one side of the Stokes curve, so they “switch on” as the Stokes
curves are crossed.

In general, an asymptotic power series of the solution of a singularly
perturbed system is divergent [24]. Nevertheless, one can truncate a divergent
asymptotic series to approximate the solution [16]. If one chooses the
truncation point to minimize the difference between the approximation and the
exact solution, the approximation error is typically exponentially small in
the asymptotic limit [15]. This is called “optimal truncation”; we use
$N_{\mathrm{opt}}$ to denote an optimal truncation point. We express the
solution of a singularly perturbed system as the sum of an optimally truncated
power series and an exponentially small contribution. Substituting this sum
into (3) produces an equation that describes the behavior of the exponentially
small contribution.

Early work that applied this idea includes [7, 8], which determined the Stokes
switching behavior in several important special functions. The analysis in [8]
demonstrated that the switching behavior depends predictably on the manner in
which the asymptotic power series diverges. Subsequent work in [10]
established techniques known as “hyperasymptotics” to further reduce the
exponentially small error that is generated by truncating asymptotic power
series. See [9] for a summary and discussion of the results in [8, 10].

In the present study, we apply an exponential asymptotic method that was
developed in [17, 75]. We express the solution $g$ of the governing equation
as an asymptotic power series

$g\sim\sum_{j=0}^{\infty}\eta^{rj}g_{j}\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (4)

where $r$ is the number of times that one needs to differentiate $g_{j-1}$ to
obtain $g_{j}$.

We first substitute the series (4) into the governing equations (3), and we
then asymptotically match terms in 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. If the series terms have singular points, this
repeated differentiation ensures that the series terms diverge in a
predictable fashion that is called a “factorial-over-power divergence” [24].
As $j$ becomes large, behavior of this form dominates the series terms.

To capture the factorial-over-power divergence, we write an ansatz for the
behavior of $g_{j}$ in the limit $j\rightarrow\infty$. This yields the so-
called “late-order terms” of an asymptotic series. As $j\rightarrow\infty$,
factorial-over-power divergent behavior dominates the series terms. Motivated
by this pattern, Chapman et. al [17] proposed applying a late-order ansatz to
approximate the form of the late-order terms, even if computing earlier terms
in a series is challenging or intractable. This ansatz is a sum of terms of
the form

$g_{j}\sim\frac{G\Gamma(rj+\gamma)}{\chi^{rj+\gamma}}\quad\mathrm{as}\quad
j\to\infty\,,$ (5)

where the parameter $\gamma$ is constant and $G$ and $\chi$ are functions of
any variables but are independent of $j$. The function $\chi$ is known as the
“singulant” and it is equal to $0$ at one or more singularities of the
leading-order solution $g_{0}$. This ensures that the late-order ansatz for
$g_{j}$ is singular at the same location(s) as the leading order, with a
singularity strength that grows as $j$ increases. Substituting (5) into the
differential equation (3) and matching orders of $j$ allows us to determine
the functional forms of $\chi$ and $G$. We determine $\gamma$ by requiring
that the late-order behavior is consistent with the local behavior of the
leading order in the neighborhood of singular points.

We optimally truncate the power series (4) using the additional information
from the late-order ansatz (5). Optimal truncation points typically occur
after an asymptotically large number of terms, so we apply the late-order
ansatz (5) to obtain a simplified expression using the heuristic that was
described in [15]. This heuristic requires truncating the series after the
value of $N_{\mathrm{opt}}$ for which the term
$\eta^{N_{\mathrm{opt}}}g_{N_{\mathrm{opt}}}$ has the smallest magnitude.

We use the truncated series to approximate the exact solution with
exponentially small error. We write

$g=\sum_{j=0}^{N_{\mathrm{opt}}-1}\eta^{rj}g_{j}+g_{\exp}\,,$ (6)

where $g_{\exp}$ denotes the exponentially small error term and
$N_{\mathrm{opt}}$ is the optimal truncation point that we calculate using the
late-order ansatz (5).

Substituting the truncated series expression (6) into the differential
equation (3) produces an equation for the exponentially small remainder term
[75]. Away from Stokes curves, one can determine this remainder using a
straightforward application of the WKB method [38]. In the neighborhood of the
Stokes curves, we apply an exponential ansatz for the exponentially small term
$g_{\exp}$ and write

$g_{\exp}\sim\mathcal{S}G\mathrm{e}^{-\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (7)

where $\mathcal{S}$ is a function that is known as the “Stokes multiplier”.
The Stokes multiplier is approximately constant except in the neighborhood of
Stokes curves, where it varies rapidly in a neighborhood of width
$\mathcal{O}(\sqrt{\eta})$ as $\eta\to 0$. This behavior is known as “Stokes
switching”. The locations of Stokes curves are determined completely by the
form of the singulant $\chi$. As was shown in [10], Stokes curves occur only
in locations where $\chi$ is real and positive. We use exponential asymptotics
to directly calculate the exponentially small contributions to the asymptotic
behavior of solutions that appear as the Stokes curves are crossed. One cannot
express these contributions using classical asymptotic power series.

It is often the case that this exponential asymptotic approach requires the
explicit calculation of only the leading-order solution in (4) to determine
the exponentially small contributions. This is convenient because it can be
complicated or even intractable to compute series terms beyond a leading-order
expression in nonlinear problems. See [14, 15] for more details on exponential
asymptotics and their applications to nonlocal solitary waves, [8, 9] for
examples of other studies of exponential asymptotics, and [17, 75] for more
details about the particular methodology that we apply in the present paper.

##  3 A Singularly Perturbed Woodpile Chain

We consider a precompressed woodpile chain, such as the one in Fig. 2(a). We
model it by a particle chain; each particle has mass $m_{1}$ and is connected
to an outside particle of mass $m_{2}$ by a linear spring with spring constant
$k$. The governing equations of this idealized model, which we illustrated in
Fig. 2(b), are

$\displaystyle m_{1}\ddot{u}(n,t)$
$\displaystyle=[\delta_{0}+u(n-1,t)-u(n,t)]_{+}^{\alpha}-[\delta_{0}+u(n,t)-u(n+1,t)]_{+}^{\alpha}-k[u(n,t)-v(n,t)]\,,$
(8) $\displaystyle m_{2}\ddot{v}(n,t)$ $\displaystyle=k[u(n,t)-v(n,t)]\,,$ (9)

where $\delta_{0}$ is the precompression and $u(n,t)$ and $v(n,t)$,
respectively, denote the displacements of the $n$-th particle of mass $m_{1}$
and $m_{2}$ at time $t$. The choice $\alpha=3/2$ corresponds to a classical
Hertzian chain. Our analysis in this section is valid for any choice of
$\alpha$ for which one can approximate the leading-order solution by a soliton
solution of the KdV equation. In practice, this corresponds to $\alpha>1$
[30]. The subscript $+$ indicates that we evaluate the bracketed term only if
its argument is positive; it is equal to $0$ when its argument is negative.
That is, we only have a Hertzian interaction between particles that are in
physical contact with each other.

We scale the system (8)–(9) using $u=m_{1}^{1/(\alpha-1)}\hat{u}$ and
$v=m_{1}^{1/(\alpha-1)}\hat{v}$, and we rewrite the governing equations in
terms of $\hat{\delta}=\delta_{0}/m_{1}^{1/(\alpha-1)}$ and $\hat{k}=k/m_{1}$.
Setting $\eta^{2}=m_{2}/m_{1}$ gives the following scaled governing equations:

$\displaystyle\ddot{\hat{u}}(n,t)$
$\displaystyle=[\hat{\delta}+\hat{u}(n-1,t)-\hat{u}(n,t)]_{+}^{\alpha}-[\hat{\delta}+\hat{u}(n,t)-\hat{u}(n+1,t)]_{+}^{\alpha}-\hat{k}[\hat{u}(n,t)-\hat{v}(n,t)]\,,$
(10) $\displaystyle\eta^{2}\ddot{\hat{v}}(n,t)$
$\displaystyle=\hat{k}[\hat{u}(n,t)-\hat{v}(n,t)]\,.$ (11)

In our study, we assume that the heavy particles always maintain physical
contact with neighboring particles. Therefore, the argument inside a bracket
is always non-negative, and we omit the subscript $+$ in our subsequent
analysis. One can validate this assumption by checking the solution behavior
directly. For convenience, we perform our analysis on the scaled system
(10)–(11); in our subsequent notation, we omit the hats that indicate this
scaling.

We are interested in the asymptotic behavior of traveling-wave solutions of
(10)–(11) when $0<\eta\ll 1$. These solutions consist of a localized wave
core, which we generate using an asymptotic power series in $\eta$, and
exponentially small but non-decaying oscillations that we compute using
exponential asymptotic techniques.

###  3.1 Leading-Order Solution

We expand $u(x,t)$ and $v(x,t)$ as an asymptotic power series in $\eta^{2}$ in
the limit $\eta\rightarrow 0$ and write

$\displaystyle u(n,t)\sim\sum_{j=0}^{\infty}\eta^{2j}u_{j}(n,t)\,,\quad
v(n,t)\sim\sum_{j=0}^{\infty}\eta^{2j}v_{j}(n,t)\,.$ (12)

To 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 $\eta\to 0$ to obtain

$\displaystyle u_{0}(n,t)=v_{0}(n,t)\,.$ (13)

Inserting the relation (13) into (10) gives

$\displaystyle\ddot{u}_{0}(n,t)=\left[\delta+u_{0}({n-1},t)-u_{0}({n},t)\right]^{\alpha}-\left[\delta+u_{0}({n},t)-u_{0}({n+1},t)\right]^{\alpha}\,.$
(14)

One can obtain an approximation to the leading-order behavior by following the
analysis in [71]. Assume that the deformation in the chain is sufficiently
small such that $|u_{0}(n-1,t)-u_{0}(n,t)|\ll\delta$. We expand (14) as

$\displaystyle\ddot{u}_{0}({n},t)$
$\displaystyle=\alpha\delta^{\alpha-1}[u_{0}({n-1},t)-2u_{0}({n},t)+u_{0}({n+1},t)]$
$\displaystyle\quad+\frac{\alpha(\alpha-1)}{2}\delta^{\alpha-2}[(u_{0}({n-1},t)-u_{0}({n},t))^{2}-(u_{0}({n},t)-u_{0}({n+1},t))^{2}]$
$\displaystyle\quad+\mathcal{O}((|u_{0}({n-1},t)-u_{0}({n},t)|/\delta)^{3})\,.$
(15)

This expansion has only a single nonlinear term in the retained orders, so we
refer to the physical regime in which this expansion is valid as a “weakly
nonlinear regime”.

In the long-wavelength limit, in which the characteristic size $L$ of a wave
is large in comparison to the particle radius $R$, we write

$u_{0}(n,t)=u(x,t)\,,\quad x=2Rn\,,\quad u_{0}(n\pm 1)=\mathrm{e}^{\pm
2R\frac{\partial}{\partial x}}{u_{0}(n)}\,.$ (16)

Substituting (16) into (15) and neglecting terms of order
$O(u\delta^{\alpha-1}(R/L)^{5}[(R/L)+(u/\delta)])$ yields

$u_{tt}=c_{0}^{2}u_{xx}+2c_{0}\gamma u_{xxxx}-\sigma u_{x}u_{xx}\,,$ (17)

where

$c_{0}=2R\sqrt{\alpha\delta^{\alpha-1}}\,,\quad\gamma=\frac{c_{0}R^{2}}{6}\,,\quad\sigma=2(\alpha-1)\frac{c_{0}^{2}R}{\delta}\,.$
(18)

Introducing the transformations $\tau=c_{0}^{3}t$ and $\xi=x-c_{0}t$ and
neglecting the higher-order term $u_{\tau\tau}$ gives the KdV equation

$c_{0}^{3}w_{\tau}+\gamma w_{\xi\xi\xi}+\frac{\sigma}{2c_{0}}ww_{\xi}=0\,,$
(19)

where $w=-u_{\xi}$.

We use solutions of (19) to approximate the leading-order behavior of $u(x,t)$
and then obtain the leading-order behavior of $v(x,t)$ using (13). We are
particularly interested in solutions that are perturbations of a leading-order
solitary wave, and we thus select $w(\xi,\tau)$ in (19) to be the KdV
solitary-wave solution

$w(\xi,\tau)=A\,\mathrm{sech}^{2}\bigg{[}\sqrt{\frac{\sigma
A}{24c_{0}\gamma}}\bigg{(}\xi-\frac{\sigma}{6c_{0}^{4}}A\tau\bigg{)}\bigg{]}\,.$
(20)

Integrating (20) and transforming $\xi$, $\tau$ into $n$, $t$ yields the
leading-order behavior

$\displaystyle
u_{0}(n,t)=-\frac{\delta\epsilon}{\alpha-1}\tanh\left(\epsilon(n-c_{\epsilon}t)\right)+\mathcal{O}(\epsilon^{3/2})\,,$
(21)

where

$\epsilon=\frac{2R}{L}\,,\quad
c_{\epsilon}=\delta^{(\alpha-1)/2}\sqrt{\alpha}+\frac{\sqrt{\alpha}}{6}\epsilon^{2}\delta^{(\alpha-1)/2}\,,\quad
L=\sqrt{\frac{24c_{0}\gamma}{\sigma A}}.$ (22)

The long-wavelength limit corresponds to $\epsilon\rightarrow 0$, and we thus
refer to $\epsilon$ as the “long-wavelength parameter”.

Motivated by (21), we use the co-moving frame $\xi=n-c_{\epsilon}t$. In terms
of $\xi$, the leading-order behavior is

$\displaystyle u_{0}(\xi)$
$\displaystyle=-\frac{\delta\epsilon}{\alpha-1}\tanh\left(\epsilon\xi\right)+\mathcal{O}(\epsilon^{3/2})\quad\mathrm{as}\quad\epsilon\rightarrow
0\,,$ (23) $\displaystyle v_{0}(\xi)$ $\displaystyle=u_{0}(\xi)\,.$ (24)

The leading-order behavior has singularities at
$\xi_{N\pm}=\pm\mathrm{i}(2N+1)\pi/2\epsilon$ for $N\in\mathbb{Z}$. The
singularities that are closest to the real axis occur at
$\xi_{0\pm}=\pm\mathrm{i}\pi/2\epsilon$ and are associated with singulants
with the smallest values of $|\chi|$ when evaluated at real $\xi$.
Consequently, these terms dominate the late-order behavior as
$n\rightarrow\infty$. Near these singularities, we calculate that

$\displaystyle
v_{0}(\xi)\sim-\frac{\delta(\alpha-1)^{-1}}{\xi-\xi_{0\pm}}\quad\mathrm{as}\quad\xi\rightarrow\xi_{0\pm}$
(25)

and that $u_{0}$ has the same behavior as $v_{0}$ in the neighborhood of the
singularity.

###  3.2 Late-Order Terms

Writing the governing equations (10)–(11) in terms of $\xi$ and matching at
each order of $\eta$ gives

$\displaystyle c_{\epsilon}^{2}u_{j}^{\prime\prime}(\xi)$
$\displaystyle=\alpha[u_{j}(\xi-1)-u_{j}(\xi))(\delta+u_{0}(\xi-1)-u_{0}(\xi)]^{\alpha-1}$
$\displaystyle\quad-\alpha[u_{j}(\xi)-u_{j}(\xi+1)][\delta+u_{0}(\xi)-u_{0}(\xi+1)]^{\alpha-1}-k[u_{j}(\xi)-v_{j}(\xi)]+\ldots\,,$
(26) $\displaystyle c_{\epsilon}^{2}v_{j-1}^{\prime\prime}(\xi)$
$\displaystyle=k[u_{j}(\xi)-v_{j}(\xi)]\,,$ (27)

where we omit the terms that are products that include $u_{j-k}$ with $k>1$.
These terms are subdominant with respect to the terms that we retain in the
limit $j\to\infty$. We retain all terms that include $u_{j}$, $v_{j}$, and
derivatives of $u_{j}$ and $v_{j-1}$. Given the general form of the factorial-
over-power ansatz (5), we conclude that the 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
(28).

In principle, one can apply (3.2) and (27) recursively to obtain terms in the
series (12) up to arbitrarily large values of $j$. This process is challenging
technically because it requires solving both a differential–difference
equation (3.2) and an algebraic equation (27) at each order. Fortunately, this
is not necessary for our analysis. Additionally, obtaining terms up to
arbitrary order does not reveal the presence of oscillations in the far field,
where $\xi\to-\infty$, because the oscillations are exponentially small in the
singularly perturbed limit $\eta\to 0$. Instead, we obtain the asymptotic form
of the late-order terms as part of the exponential asymptotic process that we
use to calculate the behavior of these oscillations.

The late-order ansatz consists of a sum of terms of the forms

$u_{j}\sim\frac{U(\xi)\Gamma(2j+\beta_{1})}{\chi(\xi)^{2j+\beta_{1}}}\,,\quad
v_{j}\sim\frac{V(\xi)\Gamma(2j+\beta_{2})}{\chi(\xi)^{2j+\beta_{2}}}\quad\mathrm{as}\quad
j\rightarrow\infty\,.$ (28)

To ensure that late-order terms have singularities at the same locations as
the leading-order solution, we set $\chi=0$ at a particular choice of
$\xi=\xi_{0\pm}$. The full late-order expression is then the sum of the late-
order contributions from each of the singularities.

For sufficiently large $j$, the terms of the asymptotic series (28) diverge in
a factorial-over-power fashion, which confirms that $u_{j}\gg u_{j-k}$ and
$v_{j}\gg v_{j-k}$ as $j\to\infty$ for $k>0$. Additionally, inserting the
late-order ansatz (28) into (3.2) shows that only $\beta_{1}+2=\beta_{2}$
produces a nontrivial asymptotic balance, implying that
$u_{j}=\mathcal{O}(v_{j-1})$ as $j\rightarrow\infty$ and hence that $v_{j}\gg
u_{j}$ as $j\to\infty$.

Inserting the late-order ansatz (28) into (27) gives

$\displaystyle\frac{c_{\epsilon}^{2}(\chi^{\prime}(\xi))^{2}V(\xi)\Gamma(2j+\beta_{2})}{\chi(\xi)^{2j+\beta_{2}}}$
$\displaystyle-\frac{2c_{\epsilon}^{2}\chi^{\prime}(\xi)V^{\prime}(\xi)\Gamma(2j+\beta_{2}-1)}{\chi(\xi)^{2j+\beta_{2}-1}}$
$\displaystyle-\frac{c_{\epsilon}^{2}\chi^{\prime\prime}(\xi)V(\xi)\Gamma(2j+\beta_{2}-1)}{\chi(\xi)^{2j+\beta_{2}-1}}+\cdots=-\frac{kV(\xi)\Gamma(2j+\beta_{2})}{\chi(\xi)^{2j+\beta_{2}}}+\cdots\,,$
(29)

where the omitted terms are no larger than $\mathcal{O}(v_{j-1})$ in the
$j\rightarrow\infty$ limit.

Matching terms at $\mathcal{O}(v_{j})$ in the $j\rightarrow\infty$ limit gives
the singulant equation $c_{\epsilon}^{2}(\chi^{\prime}(\xi))^{2}=-k$. This
implies that $\chi^{\prime}(\xi)=\pm\mathrm{i}\sqrt{k}/c_{\epsilon}$, which we
integrate to obtain

$\chi(\xi)=\pm\frac{\mathrm{i}\sqrt{k}(\xi-\xi_{0\pm})}{c_{\epsilon}}\,.$ (30)

Stokes switching can occur only if $\mathrm{Re}(\chi)>0$, which corresponds to
the positive sign choice for $\xi_{0+}$ and the negative sign choice for
$\xi_{0-}$. Consequently, we retain these solutions and ignore the solutions
that are associated with the remaining sign choices, as those can never appear
in the asymptotic solution.

Matching terms at $\mathcal{O}(v^{\prime}_{j-1})$ gives the prefactor equation
$2V^{\prime}(\xi)\chi^{\prime}(\xi)=0$, so $V$ is a constant, with a value
that depends on the choice of singularity. For clarity, we subsequently use
$\Lambda_{\pm}$ to denote the constant prefactor that is associated with a
singularity at $\xi=\xi_{0\pm}$. For the singular late-order behavior to be
consistent with the local behavior of the leading-order solution (25), we
calculate that $\beta_{2}=1$. In Appendix A.1, we perform a local expansion of
the solutions $u(\xi)$ and $v(\xi)$ in a neighborhood of size
$\mathcal{O}(\epsilon)$ near the singularity, and we use asymptotic matching
to obtain

$\Lambda_{+}=-\frac{\mathrm{i}\delta\sqrt{k}}{(\alpha-1)c_{\epsilon}}\,.$ (31)

We thereby fully determine the asymptotic behavior of $v_{j}$ in the
$j\rightarrow\infty$ limit.

In Section 3.3, we give a detailed calculation of the exponentially small
oscillation that is associated with the singularity at $\xi=\xi_{0+}$. At the
conclusion of our analysis in Section 3.3, we state the corresponding
contribution that arises from the singularity at $\xi=\xi_{0-}$; this
contribution is the complex conjugate of the contribution from $\xi=\xi_{0+}$.

###  3.3 Stokes Switching

We truncate the asymptotic series after $N$ terms to obtain

$u(\xi)=\sum_{j=0}^{N-1}\eta^{2j}u_{j}(\xi)+S_{N}(\xi)\,,\quad
v(\xi)=\sum_{j=0}^{N-1}\eta^{2j}v_{j}(\xi)+R_{N}(\xi)\,,$ (32)

where $S_{N}$ and $R_{N}$ are the remainder terms that we obtain by truncating
the series. These remainder terms are exponentially small if we optimally
truncate the series; we again denote the optimal truncation point as
$N=N_{\mathrm{opt}}$.

To optimally truncate the power series (12), we follow the heuristic approach
of [15]. This approach requires truncating the series at the smallest term of
the series. See [15] for a detailed discussion of the validity of this
approach. We locate the smallest term by taking the derivative of the term
$\eta^{2N}v_{N}(\xi)$ with respect to $N$ to obtain

$\frac{\partial}{\partial N}\left|\eta^{2N}v_{N}\right|\sim
2\eta^{2N}\frac{|V|\Gamma(2N+1)}{|\chi|^{2N+1}}\left(\log\eta-\log|\chi(\xi)|+\log(2N+1)\right)\,.$
(33)

The smallest term occurs at the point at which the derivative term in (33) is
equal to $0$. This requires that

$\log(2N+1)\sim\log|\chi|-\log\eta\,.$ (34)

We thus obtain

$N_{\mathrm{opt}}=|\chi|/2\eta+\omega\,,$ (35)

where we choose $\omega\in[0,1)$ such that $N_{\mathrm{opt}}$ is an integer.

Inserting (32) into the governing equations (10)–(11) yields

$\displaystyle c_{\epsilon}^{2}S^{\prime\prime}(\xi)$ $\displaystyle\sim-
kR_{N}(\xi)\,,$ (36)
$\displaystyle\eta^{2}c_{\epsilon}^{2}R^{\prime\prime}(\xi)+\eta^{2N}c_{\epsilon}^{2}v_{N-1}^{\prime\prime}(\xi)$
$\displaystyle\sim-kR_{N}(\xi)\quad\mathrm{as}\quad\eta\rightarrow 0\,,$ (37)

where the omitted terms are smaller than those that we retain in the limit
$\eta\rightarrow 0$. Equation (37) decouples from (36), so we can study it
independently. Applying the late-order ansatz and rearranging gives

$\eta^{2}c_{\epsilon}^{2}R_{N}^{\prime\prime}+kR_{N}\sim-\frac{\Lambda\eta^{2N}(\chi^{\prime})^{2}\Gamma(2N+1)}{\chi^{2N+1}}\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (38)

The right-hand side of (38) is exponentially small except in a neighborhood
around the Stokes curve. Away from the Stokes curve, we use the WKB method to
obtain

$\displaystyle R_{N}\sim
C\mathrm{e}^{-\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow 0\,,$ (39)

where $C$ is a constant that we need to determine.

To capture the variation in the neighborhood of a Stokes curve, we adapt the
form of (39) to include a Stokes-switching parameter $A(\xi)$, which is
constant except in the neighborhood of the Stokes curve. In such a
neighborhood, $R_{N}$ takes the form

$R_{N}(\xi)\sim
A(\xi)\mathrm{e}^{-\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow 0\,.$ (40)

Inserting (40) into (38) and rearranging yields

$\frac{\mathrm{d}A}{\mathrm{d}\xi}\sim\frac{\Lambda\chi^{\prime}\eta^{2N-1}\Gamma(2N+1)}{2\chi^{2N+1}}\,\mathrm{e}^{\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (41)

Writing $N_{\mathrm{opt}}$ in terms of $\chi$, expanding the gamma function
using Stirling’s formula, and transforming to make $\chi$ the independent
variable gives

$\frac{\mathrm{d}A}{\mathrm{d}\chi}\sim\frac{\Lambda\sqrt{\pi}\eta^{|\chi|/\eta+2\omega-1}(|\chi|/\eta)^{|\chi|/\eta+2\omega-1/2}\mathrm{e}^{-|\chi|/\eta}}{\sqrt{2}\chi^{|\chi|/\eta+2\omega+1}}\,\mathrm{e}^{\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (42)

We transform (42) into polar coordinates using
$\chi=\rho\mathrm{e}^{\mathrm{i}\theta}$ and consider variations in the
angular direction. After some simplification, we obtain

$\frac{\mathrm{d}A}{\mathrm{d}\theta}\sim{\mathrm{i}\Lambda}\sqrt{\frac{\pi\rho}{2\eta^{3}}}\exp\left(\frac{\rho}{\eta}(\mathrm{e}^{\mathrm{i}\theta}-1)-\frac{\mathrm{i}\theta\rho}{\eta}-2\mathrm{i}\omega\theta\right)\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (43)

The right-hand side of (43) is exponentially small in $\eta$, except in the
neighborhood of $\theta=0$. Defining an inner region
$\theta=\eta^{1/2}\bar{\theta}$, we find that

$\displaystyle\frac{\,\mathrm{d}A}{\,\mathrm{d}\bar{\theta}}\sim\frac{\mathrm{i}\Lambda}{\eta}\sqrt{\frac{\pi\rho}{2}}\mathrm{e}^{-\rho\bar{\theta}^{2}/2}\,.$
(44)

By integrating (44), we see that the behavior of $A$ as the Stokes curve is
crossed is

$\displaystyle
A\sim\frac{\mathrm{i}\Lambda}{\eta}\sqrt{\frac{\pi}{2}}\int_{-\infty}^{\sqrt{\rho}\bar{\theta}}\mathrm{e}^{-s^{2}/2}\,\mathrm{d}s\quad\mathrm{as}\quad\eta\to
0\,.$ (45)

We evaluate the integral in (45) and find that the difference between the
values of $A$ on the two sides of the Stokes curve is

$[A]_{-}^{+}\sim\frac{\mathrm{i}\pi\Lambda}{\eta}\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (46)

where $[A]_{-}^{+}$ denotes the change in $A$ as the Stokes curve is crossed
from $\theta<0$ to $\theta>0$. Recalling (31) and (40), we find that the
exponentially small contribution from $\xi=\xi_{0+}$ is

$[R_{N}]_{-}^{+}\sim\frac{\delta\pi\sqrt{k}}{(\alpha-1)c_{\epsilon}\eta}\mathrm{e}^{-\mathrm{i}\sqrt{k}(\xi-\mathrm{i}\pi/2\epsilon)/c_{\epsilon}\eta}\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (47)

The exponentially small contribution from $\xi=\xi_{0-}$ is given by the
complex conjugate of (47). Therefore, the total exponentially small
contribution is

$[R_{N}]_{-}^{+}\sim\frac{\delta\pi\sqrt{k}}{(\alpha-1)c_{\epsilon}\eta}\mathrm{e}^{-\mathrm{i}\sqrt{k}(\xi-\mathrm{i}\pi/2\epsilon)/c_{\epsilon}\eta}+\mathrm{c.c.}\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (48)

where c.c. denotes the complex conjugate. We express (48) in terms of
trigonometric functions and thus write

$[R_{N}]_{-}^{+}\sim\frac{2\delta\pi\sqrt{k}}{(\alpha-1)\eta
c_{\epsilon}}\exp\left(-\frac{\pi\sqrt{k}}{2c_{\epsilon}\epsilon\eta}\right)\cos\left(\frac{\xi\sqrt{k}}{c_{\epsilon}\eta}\right)\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (49)

We are considering one-sided nanoptera for which $R_{N}=0$ in the undisturbed
upstream region. This corresponds to requiring that $R_{N}=0$ for $\theta<0$
and thus that $R_{N}\sim[R_{N}]_{-}^{+}$ for $\theta>0$. If we were studying
symmetric two-sided nanoptera, we would instead set
$R_{N}\sim-\tfrac{1}{2}[R_{N}]_{-}^{+}$ for $\theta<0$ and
$R_{N}\sim\tfrac{1}{2}[R_{N}]_{-}^{+}$ for $\theta>0$. This is consistent with
the jump in (49), and it yields an oscillation amplitude that is half of that
of the corresponding one-sided nanopteron.

###  3.4 Comparison of our Asymptotic and Computational Results

We compare the amplitude of the oscillations that we predict using our
asymptotic analysis to that from numerical simulations. We employ a symplectic
integrator in the form of a velocity Verlet algorithm [96, 2], which is
convenient for studying chains of particles. This approach 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

$\boldsymbol{x}_{n+1}=\boldsymbol{x}_{n}+\boldsymbol{v}_{n}\Delta
t+\frac{1}{2}\boldsymbol{a}_{n}(\Delta
t)^{2}\,,\quad\boldsymbol{v}_{n+1}=\boldsymbol{v}_{n}+\frac{1}{2}(\boldsymbol{a}_{n}+\boldsymbol{a}_{n+1})\Delta
t\,,$ (50)

where $\Delta t$ is the size of the time step and $\boldsymbol{x}_{n}$,
$\boldsymbol{v}_{n}$, and $\boldsymbol{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 vector
$\boldsymbol{a}_{n}$ as an algebraic function of the displacement vector
$\boldsymbol{x}_{n}$ using (10)–(11).

We truncate the domain and impose periodic boundary conditions. Instead of
directly computing the absolute displacements and velocities of the particles,
we perform our computations using the relative displacements and velocities,
which we compute by calculating the differences in the positions and
velocities of adjacent particles. This is convenient for our computations
because the far-field behavior of the leading-order traveling wave approaches
$0$ in both directions in these coordinates. The relative displacements of the
heavy particles and light particles are $r_{1}$ and $r_{2}$, respectively,
where

$r_{1}(n,t)=u(n+1,t)-u(n,t)\,,\quad r_{2}(n,t)=v(n+1,t)-v(n,t)\,.$ (51)

In relative coordinates, the governing equations are

$\displaystyle\ddot{r_{1}}(n,t)$ $\displaystyle=2[\delta-
r_{1}(n,t)]_{+}^{\alpha}-[\delta-r_{1}(n-1,t)]_{+}^{\alpha}-[\delta-
r_{1}(n+1,t)]_{+}^{\alpha}-k[r_{1}(n,t)-r_{2}(n,t)]\,,$ (52)
$\displaystyle\eta^{2}\ddot{r_{2}}(n,t)$
$\displaystyle=k[r_{1}(n,t)-r_{2}(n,t)]\,.$ (53)

The truncated domain has $M=2^{10}$ particles, whose indices are
$n\in\\{-M/2+1,\ldots,M/2\\}$. The initial condition is the leading-order
traveling-wave solution (23)–(24). To avoid interactions between the far-field
oscillations and the leading-order traveling wave in the periodic domain, we
apply a window function to the solution, as in [31, 56], at each iteration. We
obtain the window by multiplying $r(n,t)$ and $\dot{r}(n,t)$ by a function
$W(n-n_{\max}+M/8)$, where

$W(k)=\left\\{\begin{array}[]{ll}1\,,&\quad|K|\leq\tfrac{5M}{16}\\\
1-\tfrac{8}{M}\left(|K|-\tfrac{5M}{16}\right)\,,&\quad\tfrac{5M}{16}<|K|\leq\tfrac{7M}{16}\\\
0\,,&\quad\tfrac{7M}{16}<|K|\leq\tfrac{M}{2}\end{array}\right.$ (54)

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 trailing oscillations in simulations of FPUT systems, as
any disturbances that are associated with the windowing must travel more
slowly than the leading-order solitary wave and hence than the window itself
(see [29]). We do not perform a comparable analysis for a woodpile chain, but
we do not observe any discernible differences in the amplitudes of the
trailing oscillations as a result of our windowing.

Figure 4: Numerically-calculated relative displacement of a single light
particle at some index $n$ in a woodpile chain with an interaction exponent of
$\alpha=1.5$. The mass-ratio parameter is $\eta=0.4$, the precompression of
the chain gives $\delta=5$, and the long-wavelength parameter is
$\epsilon=0.3$.

In Fig. 4, we show a numerically-calculated profile of one particle in a
woodpile chain. In this simulation, we use an interaction exponent of
$\alpha=1.5$, a mass-ratio parameter of $\eta=0.4$, a precompression of
$\delta=5$, and a long-wavelength parameter of $\epsilon=0.3$. The curve in
the figure shows the relative displacement of a single light particle over a
range of times $t$. The particle is 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 small amplitude. In Fig. 5, we measure this oscillation
amplitude from numerical simulations and compare it to our asymptotic
approximations for a range of parameter choices.

(a)  (b)  (c)  (d)

(e)  (f)  (g)  (h)

(i)  (j)  (k)  (l)

Figure 5: Comparison of our asymptotic approximations and numerical
simulations of the amplitude of exponentially small oscillations in the far
field (i.e., as $\xi\to-\infty$) for woodpile chains with interaction
exponents of (a) $\alpha=1.5$, (b) $\alpha=2$, and (c) $\alpha=2.5$. For each
$\alpha$, we show simulations for precompressions of $\delta=5$, $\delta=7.5$,
$\delta=10$, and $\delta=15$. For each choice of $\delta$, we set the long-
wavelength parameter to be $\epsilon=0.2$, $\epsilon=0.25$, and
$\epsilon=0.3$. The accuracy of our asymptotic approximation improves for
progressively smaller values of the mass-ratio parameter $\eta$, as expected
for an $\eta\rightarrow 0$ asymptotic approximation. Additionally, the
approximation improves as we increase $\delta$. This is a consequence of using
a leading-order approximation for the weakly nonlinear regime; decreasing
$\delta$ approaches a configuration that needs to be treated as strongly
nonlinear. Our approximation also improves as we increase $\alpha$, although
the reason for this improvement is not apparent.

In our numerical simulations, we calculate the amplitude of the far-field
waves for interaction exponents of $\alpha=1.5$, $\alpha=2$ and $\alpha=2.5$.
For each $\alpha$, we show simulations with precompressions of $\delta=5$,
$\delta=7.5$, $\delta=10$, and $\delta=15$. For each choice of $\delta$, we
set the long-wavelength parameters to be $\epsilon=0.2$, $\epsilon=0.25$, and
$\epsilon=0.3$. We choose the values of $\eta$ so that the far-field waves are
large enough to be detected numerically, but small enough that our
small-$\eta$ asymptotic approximation accurately describes the behavior of the
solution. We run each of our simulations for enough time so that any time
variation in the far-field oscillation amplitude is not visible in the
solution. We show the results of these simulations, which we compare to the
asymptotic results from (49), in Fig. 5.

It is apparent from Fig. 5 that our asymptotic approximation captures the far-
field wave amplitude effectively for a wide range of the parameters and that
the approximation error decreases as $\eta\rightarrow 0$, as we expect from
the asymptotic nature of the approximation. For a fixed $\epsilon$, our
asymptotic approximation becomes progressively more accurate for progressively
larger $\delta$. However, for progressively smaller $\delta$, the system
approaches the strongly nonlinear regime and the long-wavelength approximation
becomes less accurate, which in turn causes the late-order approximation to
become less accurate. Our asymptotic approximation also appears to become more
accurate for progressively larger $\alpha$, 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 the oscillations that appear in woodpile chains
as a consequence of Stokes phenomena and that the employed exponential
asymptotic analysis is a useful technique for studying nanoptera in these
chains.

##  4 A Singularly Perturbed Diatomic Hertzian Chain

We now consider a singularly perturbed diatomic Hertzian chain. This system is
related to the woodpile chains that we analyzed in Section 3. Once we derive a
leading-order solution for a traveling wave in a diatomic Hertzian chain, we
can perform a similar exponential asymptotic analysis. Obtaining the singulant
equation for diatomic Hertzian chains requires solving a complicated integral
expression; this is significantly more complicated than the corresponding step
in our analysis of woodpile chains.

We consider a diatomic Hertzian chain, which is governed by (55)–(56). We
denote the masses in the chain by $m(j)=m_{1}$ for even $j$ and $m(j)=m_{2}$
for odd $j$. We consider a small ratio between the masses of the two types of
particles. This choice amounts to a precompressed Hertzian analog of the
diatomic Toda and FPUT chains that were studied in [57, 56]. Those previous
studies used far-field asymptotic behavior to identify an orthogonality
condition, which yields values of the mass-ratio parameter $\eta$ that cause
the far-field waves to cancel and produce a genuine solitary wave.

We derive an asymptotic approximation in a diatomic Hertzian chain to identify
an orthogonality condition that finds some (but not all) values of $\eta$ that
produce solitary waves. Specifically, our asymptotic approach systematically
fails to identify every second value of $\eta$ that causes the far-field
oscillations to cancel.

Much of our analysis in this section is similar to that in Section 3 and in
[57, 56]. Therefore, we present only an outline of our calculations.

The governing equations for a diatomic chain of particles is

$\displaystyle m_{1}\ddot{u}(2n,t)$
$\displaystyle=[\delta_{0}+v(2n-1,t)-u(2n,t)]_{+}^{\alpha}-[\delta_{0}+u(2n,t)-v(2n+1,t)]_{+}^{\alpha}\,,$
(55) $\displaystyle m_{2}\ddot{v}(2n+1,t)$
$\displaystyle=[\delta_{0}+u(2n,t)-v(2n+1,t)]_{+}^{\alpha}-[\delta_{0}+v(2n+1,t)-u(2n+2,t)]_{+}^{\alpha}\,,$
(56)

where $u$ and $v$, respectively, represent the displacement of even and odd
particles, the $+$ subscript has the same meaning as in (10)–(11), and we set
$\alpha=3/2$ so that we consider Hertzian interactions.444 For diatomic
Hertzian chains, the form of the singulant depends on the value of $\alpha$.
By contrast, for woodpile chains, the singulant is independent of $\alpha$.
Therefore, in Section 3, we permitted $\alpha$ to take any value that produces
leading-order solitary waves, as it did not change our subsequent analysis.
For diatomic chains, we must pick a specific value of $\alpha$. Therefore, we
restrict our attention to Hertzian interactions, which are the most common
case realized in existing laboratory experiments. We apply the same scalings
as we did for woodpile chains, so we let $u=m_{1}^{2}\hat{u}$ and
$v=m_{1}^{2}\hat{v}$. We rewrite (55)–(56) in terms of our scaled variables
and a scaled precompression parameter $\hat{\delta}=\delta_{0}/m_{1}^{2}$.
Setting $\eta^{2}=m_{2}/m_{1}$ yields the scaled system

$\displaystyle\ddot{\hat{u}}(2n,t)$
$\displaystyle=[\hat{\delta}+\hat{v}(2n-1,t)-\hat{u}(2n,t)]_{+}^{3/2}-[\hat{\delta}+\hat{u}(2n,t)-\hat{v}(2n+1,t)]_{+}^{3/2}\,,$
(57) $\displaystyle\eta^{2}\ddot{\hat{v}}(2n+1,t)$
$\displaystyle=[\hat{\delta}+\hat{u}(2n,t)-\hat{v}(2n+1,t)]_{+}^{3/2}-[\hat{\delta}+\hat{v}(2n+1,t)-\hat{u}(2n+2,t)]_{+}^{3/2}\,.$
(58)

As in the woodpile chains, we assume that the precompression is sufficiently
strong that the particles remain in contact and thus that the quantities in
the square brackets in (57)–(58) are never negative. Therefore, we omit the
$+$ subscript in our subsequent notation. In our analysis of the scaled system
(57)–(58), we also omit the hats from the variables and the parameters for
notational convenience.

As in the woodpile chains of Section 3, we study traveling-wave solutions of
(57)–(58) when the mass-ratio parameter $\eta$ is small. Specifically, we
consider $0<\eta\ll 1$. As before, the solution consists of a localized wave
core and exponentially small, non-decaying oscillations. We again use
exponential asymptotic methods to study these oscillations.

###  4.1 Leading-Order Solution

The first step to construct a nanopteron solution of (57)–(58) 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 $\eta^{2}$ in the
limit $\eta\rightarrow 0$. This yields

$\displaystyle u(n,t)\sim\sum_{j=0}^{\infty}\eta^{2j}u_{j}(n,t)\,,\quad
v(n,t)\sim\sum_{j=0}^{\infty}\eta^{2j}v_{j}(n,t)\,.$ (59)

Inserting the series expansions (59) into (58) and matching at leading order
in the $\eta\rightarrow 0$ limit gives

$\displaystyle v_{0}(2n+1,t)=\tfrac{1}{2}[u_{0}(2n,t)+u_{0}(2n+2,t)]\,.$ (60)

Inserting (60) into (57) yields

$\displaystyle\ddot{u}_{0}(2n,t)=\bigg{(}\delta+\frac{u_{0}({2n-2},t)-u_{0}({2n},t)}{2}\bigg{)}^{3/2}-\bigg{(}\delta+\frac{u_{0}({2n},t)-u_{0}({2n+2},t)}{2}\bigg{)}^{3/2}\,.$
(61)

Using a similar analysis as in Section 3.1, we obtain an approximation to the
leading-order solution by assuming that the displacement of a particle is
small in comparison to the precompression parameter $\delta$ and that the
characteristic length scale of the leading-order wave is large in comparison
to the radius of the spherical particles. We again define a long-wavelength
parameter $\epsilon$. The leading-order solution is approximately

$\displaystyle
u_{0}(\xi)=-4\delta\epsilon\tanh\left(\epsilon\xi/2\right)+\mathcal{O}(\epsilon^{3/2})\,,\quad
v_{0}(\xi)=\tfrac{1}{2}(u_{0}(\xi+1)+u_{0}(\xi-1))\,,$ (62)

where the variable $\xi$ defines the co-moving frame

$\displaystyle\xi=n-c_{\epsilon}t\,,\quad
c_{\epsilon}=\delta^{1/4}\sqrt{3}+\frac{\delta^{1/4}\epsilon^{2}}{2\sqrt{3}}\,.$
(63)

We now extend the leading-order solution (62) into the complex plane. We see
that $v_{0}(\xi)$ is singular at $\xi=\xi_{N\pm}$ in the complex plane, where

$\displaystyle\xi_{N\pm}=\frac{{(2N-1)}\pi\mathrm{i}}{\epsilon}\pm 1\,,\quad
N\in\mathbb{Z}\,.$ (64)

The behavior near the singularities is

$\displaystyle
v_{0}(\xi)\sim-4\delta(\xi-\xi_{N\pm})^{-1}\quad\mathrm{as}\quad\xi\to\xi_{N\pm}\,.$
(65)

As with woodpile chains, the dominant contribution to the late-order behavior
arises from the singularities that are closest to the real axis; this occurs
when $N=0$ and $N=1$ in (64). We illustrate the positions of these
singularities in Fig. 6. In our subsequent analysis, we show the details for
calculating the contribution from the singularity at $\xi=\xi_{1-}$. We then
state the analogous results for the contributions from the singularities at
$\xi=\xi_{1+}$, $\xi_{0-}$, and $\xi_{0+}$.

Figure 6: Singularities of $v_{0}(\xi)$ in (64) that contribute to the
asymptotic form of the far-field oscillations of (90). We give a detailed
discussion of the analysis of the oscillations that arise from the singularity
at $\xi=\xi_{1-}$. We then state the contributions from the other three
singularities.

###  4.2 Terms in the Late-Order Series

Inserting the series (59) into the system (57)–(58) and matching at each order
of $\eta$ gives a recurrence relation for $j\geq 2$. The recurrence relation
is

$\displaystyle c_{\epsilon}^{2}u_{j}^{\prime\prime}(\xi)=$
$\displaystyle\tfrac{3}{2}(v_{j}(\xi-1)-u_{j}(\xi))(\delta+v_{0}(\xi-1)-u_{0}(\xi))^{1/2}$
$\displaystyle-\tfrac{3}{2}(u_{j}(\xi)-v_{j}(\xi+1))(\delta+u_{0}(\xi)-v_{0}(\xi+1))^{1/2}$
$\displaystyle+\tfrac{3}{4}(v_{j-1}(\xi-1)-u_{j-1}(\xi))(v_{1}(\xi-1)-u_{1}(\xi))[\delta+(v_{0}(\xi-1)-u_{0}(\xi))]^{-1/2}$
$\displaystyle-\tfrac{3}{4}(u_{j-1}(\xi)-v_{j-1}(\xi+1))(u_{1}(\xi)-v_{1}(\xi+1))[\delta+(u_{0}(\xi)-v_{0}(\xi+1))]^{-1/2}+\cdots\,,$
(66) $\displaystyle c_{\epsilon}^{2}v_{j-1}^{\prime\prime}(\xi)=$
$\displaystyle\tfrac{3}{2}(u_{j}(\xi-1)-v_{j}(\xi))(\delta+u_{0}(\xi-1)-v_{0}(\xi))^{1/2}$
$\displaystyle-\tfrac{3}{2}(v_{j}(\xi)-u_{j}(\xi+1))(\delta+v_{0}(\xi)-u_{0}(\xi+1))^{1/2}$
$\displaystyle+\tfrac{3}{4}(u_{j-1}(\xi-1)-v_{j-1}(\xi))(u_{1}(\xi-1)-v_{1}(\xi))[\delta+(u_{0}(\xi-1)-v_{0}(\xi))]^{-1/2}$
$\displaystyle-\tfrac{3}{4}(v_{j-1}(\xi)-u_{j-1}(\xi+1))(v_{1}(\xi)-u_{1}(\xi+1))[\delta+(v_{0}(\xi)-u_{0}(\xi+1))]^{-1/2}+\cdots\,,$
(67)

where we omit the terms that are products that include both $u_{j-k}$ and
$v_{j-k}$ with $k>1$. The terms that we retain are those that contain
$u_{j-1}$, $v_{j-1}$, $u_{j}$, and $v_{j}$. All of the omitted terms are
subdominant in comparison to the retained terms as $j\to\infty$.

We again use a factorial-over-power ansatz to approximate the late-order
terms, so we write

$u_{j}\sim\frac{U(\xi)\Gamma(2j+\beta_{1})}{\chi(\xi)^{2j+\beta_{1}}}\,,\quad
v_{j}\sim\frac{V(\xi)\Gamma(2j+\beta_{2})}{\chi(\xi)^{2j+\beta_{2}}}\quad\mathrm{as}\quad
j\rightarrow\infty\,,$ (68)

where $\beta_{1}$ and $\beta_{2}$ are constants. Inserting the late-order
ansatz (68) into (4.2) shows that only $\beta_{1}=\beta_{2}-2$ produces a
nontrivial asymptotic balance.

####  4.2.1 Calculating $\chi$

Inserting the ansatz (68) into the recurrence relation (4.2)–(67) and matching
at $\mathcal{O}(v_{j})$ as $j\to\infty$ gives

$c_{\epsilon}^{2}(\chi^{\prime})^{2}=-3\bigg{(}\delta+\frac{u_{0}(\xi-1)-u_{0}(\xi+1)}{2}\bigg{)}^{1/2}\,.$
(69)

Integrating (69) and recalling that $\chi=0$ at the singularity location
$\xi=\xi_{1-}$, we obtain

$\chi=\pm\frac{\sqrt{3}\mathrm{i}}{c_{\epsilon}}\int_{\mathcal{C}}\bigg{(}\delta+\frac{(u_{0}(s-1)-u_{0}(s+1))}{2}\bigg{)}^{1/4}\,\mathrm{d}s\,,$
(70)

where $\mathcal{C}$ is a contour from $\xi_{1-}$ to $\xi$. Although one can
select any contour, it is convenient to divide the contour into a vertical
component $\mathcal{C}_{1}$ and a horizontal component $\mathcal{C}_{2}$. We
show our contour in Fig. 7.

The contribution to the integral in (70) from $\mathcal{C}_{1}$ has both real
and imaginary components, but the contribution to the integral from
$\mathcal{C}_{2}$ is purely imaginary. The real and imaginary parts of $\chi$
are

$\displaystyle\mathrm{Re}(\chi)=\pm\mathrm{Re}\bigg{\\{}\frac{\sqrt{3}\mathrm{i}}{c_{\epsilon}}\int_{\mathcal{C}_{1}}\bigg{(}\delta+$
$\displaystyle\frac{u_{0}(s-1)-u_{0}(s+1)}{2}\bigg{)}^{1/4}\,\mathrm{d}s\bigg{\\}}\,,$
(71)
$\displaystyle\mathrm{Im}(\chi)=\pm\mathrm{Im}\bigg{\\{}\frac{\sqrt{3}\mathrm{i}}{c_{\epsilon}}\int_{\mathcal{C}_{1}}\bigg{(}\delta+$
$\displaystyle\frac{u_{0}(s-1)-u_{0}(s+1)}{2}\bigg{)}^{1/4}\,\mathrm{d}s\bigg{\\}}$
$\displaystyle\pm\frac{\sqrt{3}\mathrm{i}}{c_{\epsilon}}\int_{\mathcal{C}_{2}}\bigg{(}\delta+\frac{u_{0}(s-1)-u_{0}(s+1)}{2}\bigg{)}^{1/4}\,\mathrm{d}s\,,$
(72)

where the sign choices are either all positive or all negative. We see that
$\mathrm{Re}(\chi)$ is constant for real-valued $\xi$, so the far-field
oscillations have a constant amplitude. Because
$|u_{0}(\xi-1)-u_{0}(\xi+1)|\rightarrow 0$ as $|\xi|\rightarrow\infty$, we
also see that $\chi^{\prime}(\xi)$ tends to a constant value and that
$\mathrm{Im}(\chi)$ depends linearly on $\xi$ in this limit. Consequently, the
associated far-field oscillations must tend to a constant wavelength.

Figure 7: The contour of integration for (70) that connects the singularity at
$s=\xi_{1-}$ with the point $s=\xi$. We indicate the location of the
singularity using a cross, and we show the contour using a thick black curve.
We divide the contour into a vertical component $\mathcal{C}_{1}$ and a
horizontal component $\mathcal{C}_{2}$. The contribution to the contour
integral from $\mathcal{C}_{1}$ has both real and imaginary components,
whereas the contribution to the integral from $\mathcal{C}_{2}$ is purely
imaginary. This implies that $\mathrm{Re}(\chi)$ is constant for real-valued
$\xi$.

Stokes switching can occur only if $\mathrm{Re}(\chi)>0$, which corresponds to
the positive signs for $\xi_{1-}$ in both (71) and (72). Therefore, we
restrict our analysis to this choice of signs. For our analysis, it is helpful
to have an asymptotic expression for the local behavior of $\chi$ near the
singular point $\xi_{1-}$. From direct calculation using local expressions for
$u_{0}(\xi+1)$ and $u_{0}(\xi-1)$, which we give explicitly in (104), we
obtain

$\displaystyle
c_{\epsilon}^{2}(\chi^{\prime})^{2}\sim-6\sqrt{\delta}(\xi-\xi_{1-})^{-1/2}\quad\mathrm{as}\quad\xi\to\xi_{1-}\,.$
(73)

We rearrange this local expression to give

$\displaystyle\chi^{\prime}\sim\pm\frac{\sqrt{6}\delta^{1/4}\mathrm{i}}{c_{\epsilon}}(\xi-\xi_{1-})^{-1/4}\quad\mathrm{as}\quad\xi\to\xi_{1-}\,.$
(74)

We choose the positive sign to be consistent with the full expression for
$\chi$, and we integrate to obtain

$\displaystyle\chi\sim\frac{4\sqrt{6}\delta^{1/4}\mathrm{i}}{3c_{\epsilon}}(\xi-\xi_{1-})^{3/4}\quad\mathrm{as}\quad\xi\to\xi_{1-}\,.$
(75)

####  4.2.2 Calculating $V$ and $\beta_{2}$

Matching (4.2)–(67) at order $\mathcal{O}(v_{j-1}^{\prime})$ as
$j\rightarrow\infty$ gives

$\displaystyle V=\frac{\Lambda_{1-}}{\sqrt{\chi^{\prime}(\xi)}}\,,$ (76)

where $\Lambda_{1-}$ is a constant that we can determine by taking an inner
expansion of the solution in the neighborhood of $\xi=\xi_{1-}$ and matching
the outer limit of this expansion with the inner limit of the late-order
ansatz (68). We need to determine the value of $\beta_{2}$ in (68) to perform
this inner analysis.

To determine $\beta_{2}$, we combine (74) and (76) to obtain a local
expression for $V$. This expression is

$\displaystyle
V\sim\frac{\Lambda_{1-}c_{\epsilon}^{1/2}}{6^{1/4}\delta^{1/8}\sqrt{\mathrm{i}}}(\xi-\xi_{1-})^{1/8}\quad\mathrm{as}\quad\xi\to\xi_{1-}\,.$
(77)

Substituting (75) and (77) into (68) yields a local expression for the late-
order ansatz:

$\displaystyle
v_{j}(\xi)\sim\frac{\Lambda_{1-}c_{\epsilon}^{1/2}(\xi-\xi_{1-})^{1/8}\Gamma(2j+\beta_{2})}{6^{1/4}\delta^{1/8}\sqrt{\mathrm{i}}(4\sqrt{6}\delta^{1/4}\mathrm{i}(\xi-\xi_{1-})^{3/4}/3c_{\epsilon})^{2j+\beta_{2}}}\,.$
(78)

The leading-order solution $v_{0}$ (62) has a singularity of order $1$ at
$\xi=\xi_{1-}$. For (78) to be consistent with repeated differentiation of the
leading-order behavior, we require that $3\beta_{2}/4-1/8=1$, yielding
$\beta_{2}=3/2$.

Finally, we need to determine $\Lambda_{1-}$. The late-order expansion (68)
breaks down in the neighborhood of $\xi=\xi_{1-}$. In this neighborhood, we
determine a local expansion of $v_{j}(\xi)$. By matching the inner limit of
the late-order behavior [which is given in (78)] with the outer limit of the
local expansion, we obtain $\Lambda_{1-}$. We show the computational portion
of this procedure in Appendix A.2 and obtain
$c_{\epsilon}^{2}\Lambda_{1-}/\delta^{3/2}\approx 38.41$.

###  4.3 Calculations of the Remainders

As in Section 3, we truncate the asymptotic series after $N$ terms. This
yields

$\displaystyle u(\xi)=\sum_{j=0}^{N-1}\eta^{2j}u_{j}(\xi)+S_{N}(\xi)\,,\quad
v(\xi)=\sum_{j=0}^{N-1}\eta^{2j}v_{j}(\xi)+R_{N}(\xi)\,,$ (79)

where $S_{N}$ and $R_{N}$ are the remainder terms that we obtain by truncating
the series. The optimal truncation point (which we find using the same method
as in Section 3.3) is $N_{\mathrm{opt}}=|\chi|/2\eta+\omega$, where we choose
$\omega\in[0,1)$ to ensure that $N_{\mathrm{opt}}$ is an integer. Note that
$N_{\mathrm{opt}}\to\infty$ as $\eta\to 0$.

We insert the truncated series expressions from (79) into the governing
equation (58). Using the recursion relation (67) and the late-order ansatz
(68), we obtain

$\displaystyle
c_{\epsilon}^{2}\eta^{2}R_{N}^{\prime\prime}(\xi)-c_{\epsilon}^{2}\chi^{\prime}(\xi)^{2}R_{N}(\xi)\sim-\frac{\eta^{2N}\chi^{\prime}(\xi)^{2}V(\xi)\Gamma(2N+3/2)}{\chi(\xi)^{2N+3/2}}\quad\mathrm{as}\quad\eta\rightarrow
0$ (80)

after some algebra.

The right-hand side of (80) is exponentially small and one can neglect it,
except in the neighborhood of the curve $\mathrm{Im}(\chi)=0$, which is the
Stokes curve. To capture the behavior of the remainder $R_{N}$ near the Stokes
curve, we write it using the same adapted WKB ansatz as in (40). That is,

$R_{N}(\xi)\sim
A(\xi)\mathrm{e}^{-\chi/\eta}\quad\mathrm{as}\quad\eta\rightarrow 0\,,$ (81)

where $A(\xi)$ is a Stokes switching parameter that varies rapidly near the
Stokes curve and is constant outside of the rapidly varying neighborhood.
Inserting (81) into (80) yields

$\displaystyle-2A^{\prime}\chi^{\prime}\mathrm{e}^{-\chi/\eta}\sim-\frac{\eta^{2N-1}(\chi^{\prime})^{2}\Gamma(2N+3/2)}{\chi^{2N+3/2}}\,.$
(82)

We transform (82) to treat $\chi$ as the independent variable and write
$\chi=r\mathrm{e}^{i\theta}$ in polar coordinates. We fix $r$, and we consider
only variations in the angular direction. Using the optimal truncation
$N=N_{\mathrm{opt}}$ and simplifying using Stirling’s formula gives

$\displaystyle\frac{\,\mathrm{d}A}{\,\mathrm{d}\theta}\sim\frac{\mathrm{i}}{\eta^{2}}\sqrt{\frac{\pi\rho}{2}}\exp\left(\frac{\rho}{\eta}\left(\mathrm{e}^{\mathrm{i}\theta}-1\right)-\frac{\mathrm{i}\theta
r}{\eta}-\mathrm{i}\theta(1/2+2\omega)\right)\,.$ (83)

The right-hand side of (83) is exponentially small in $\eta$, except in the
neighborhood of $\theta=0$. Defining an inner region
$\theta=\eta^{1/2}\bar{\theta}$, we find that

$\displaystyle\frac{\,\mathrm{d}A}{\,\mathrm{d}\bar{\theta}}\sim\mathrm{i}\sqrt{\frac{\pi\rho}{2\eta^{3}}}\mathrm{e}^{-r\bar{\theta}^{2}/2}\,.$
(84)

By integrating (84), we see that the behavior of $A$ as the Stokes curve is
crossed is

$\displaystyle
A\sim\mathrm{i}\sqrt{\frac{\pi}{2\eta^{3}}}\int_{-\infty}^{\sqrt{\rho}\bar{\theta}}\mathrm{e}^{-s^{2}/2}\,\mathrm{d}s\,.$
(85)

Evaluating the integral in (85) and using the form of $R_{N}$ from (81), we
find that the exponentially small contribution from $\xi_{1-}$ is

$[R_{N}]_{-}^{+}\sim\frac{\mathrm{i}\pi\Lambda}{\eta^{3/2}\sqrt{\chi^{\prime}}}\mathrm{e}^{-\chi/\eta}\,.$
(86)

To capture the change across the Stokes curve, we also need to include the
contribution from $\xi_{0-}$ that is conjugate to $\xi_{1-}$; this
contribution is given by the complex conjugate of (86). Adding (86) and its
complex conjugate, we find that the change in the exponentially small
contribution as the Stokes curve is crossed from left to right is

$[R_{N}]_{-}^{+}\sim\frac{\mathrm{i}\pi\Lambda}{\eta^{3/2}\sqrt{\chi^{\prime}}}\mathrm{e}^{-\chi/\eta}+\mathrm{c.c.}\,,$
(87)

where the complex conjugate c.c. indicates the contribution from the
singularity at $\xi_{0-}$.

The overall exponentially small contribution to the asymptotic behavior of
$v(\xi)$ 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). In the
limit $\eta\rightarrow 0$, the exponentially small terms have the asymptotic
expression

$\displaystyle
v_{\exp}\sim\frac{\mathrm{i}\pi}{\eta^{3/2}}\mathrm{e}^{-\mathrm{Re}(\chi_{1-}(\xi))/\eta}$
$\displaystyle\frac{\Lambda_{1-}\mathrm{e}^{-\mathrm{i}\mathrm{Im}(\chi_{1-}(\xi))/\eta}}{\sqrt{\chi_{1-}^{\prime}(\xi)}}$
$\displaystyle+$
$\displaystyle\frac{\mathrm{i}\pi}{\eta^{3/2}}\mathrm{e}^{-\mathrm{Re}(\chi_{1+}(\xi))/\eta}\frac{\Lambda_{1+}\mathrm{e}^{-\mathrm{i}\mathrm{Im}(\chi_{1+}(\xi))/\eta}}{\sqrt{\chi_{1+}^{\prime}(\xi)}}+\mathrm{c.c}\,,$
(88)

where $\chi_{1-}$ is the singulant that is associated with $\xi=\xi_{1-}$ and
$\chi_{1+}$ is the singulant that is associated with $\xi=\xi_{1+}$. Using the
relation $\mathrm{Re}(\chi_{1+})=\mathrm{Re}(\chi_{1-})$ from (71), the
asymptotic expression (88) becomes

$v_{\exp}\sim\frac{\mathrm{i}\pi}{\eta^{3/2}}\mathrm{e}^{-\mathrm{Re}(\chi_{1-}(\xi))/\eta}\left[\frac{\Lambda_{1-}\mathrm{e}^{-\mathrm{i}\mathrm{Im}(\chi_{1-}(\xi))/\eta}}{\sqrt{\chi_{1-}^{\prime}(\xi)}}+\frac{\Lambda_{1+}\mathrm{e}^{-\mathrm{i}\mathrm{Im}(\chi_{1+}(\xi))/\eta}}{\sqrt{\chi_{1+}^{\prime}(\xi)}}\right]+\mathrm{c.c}\,.$
(89)

In Appendix A.2, we show that $\Lambda_{1+}=\mathrm{i}\Lambda_{1-}$. Far
behind the central solitary wave, $|u_{0}(\xi-1)-u_{0}(\xi+1)|\to 0$ because
the leading-order traveling wave is exponentially localized. Consequently,
$\chi_{1,\pm}^{\prime}(\xi)\sim\mathrm{i}\sqrt{3}\delta^{1/4}/c_{\epsilon}$ in
the limit $\xi\rightarrow-\infty$, as indicated in (69). This yields a
convenient trigonometric simplification:

$\displaystyle
v_{\exp}\sim\frac{4\pi\Lambda_{1-}}{\eta^{3/2}\sqrt{\chi_{1-}^{\prime}(\xi)}}\mathrm{e}^{-\mathrm{Re}(\chi_{1-}(\xi))/\eta}\cos\bigg{(}\frac{\mathrm{Im}(\chi_{1-}(\xi)-\chi_{1+}(\xi))}{2\eta}$
$\displaystyle+\frac{\pi}{4}\bigg{)}$ $\displaystyle\times\cos$
$\displaystyle\bigg{(}\frac{\mathrm{Im}(\chi_{1+}(\xi)+\chi_{1-}(\xi))}{2\eta}-\frac{\pi}{4}\bigg{)}\,.$
(90)

From the integral expression (70) for the singulant $\chi$, we see that
$\mathrm{Re}(\chi_{1,\pm})$ and
$\mathrm{Im}(\chi_{1-})-\mathrm{Im}(\chi_{1+})$ do not depend on $\xi$.
Consequently, we write the amplitude of the far-field waves as

$v_{\mathrm{amp}}\sim\frac{4\pi
c_{\epsilon}^{1/2}\Lambda_{1-}}{3^{1/4}\delta^{1/8}\eta^{3/2}}\mathrm{e}^{-\mathrm{Re}(\chi_{1-}(\xi))/\eta}\cos\bigg{(}\frac{\mathrm{Im}(\chi_{1-}(\xi)-\chi_{1+}(\xi))}{2\eta}+\frac{\pi}{4}\bigg{)}\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (91)

###  4.4 Comparison of our Asymptotic and Computational Results

As in Section 3.4, we compare our asymptotic results for the far-field
amplitude to numerical simulations. We again use the velocity Verlet algorithm
and again simulate the equations of motion on a large periodic domain with
relative coordinates. For diatomic Hertzian chains, the relative coordinate
system is given in terms of $r(n,t)$, where

$\displaystyle r(2n,t)=v(2n+1,t)-u(2n,t)\,,\quad
r(2n+1,t)=u(2n+2,t)-v(2n+1,t)\,.$ (92)

In this coordinate system, the equations of motion (57)–(58) become

$\displaystyle\ddot{r}(2n,t)$
$\displaystyle=\bigg{(}1+\frac{1}{\eta^{2}}\bigg{)}[\delta-r(2n,t)]^{3/2}_{+}-[\delta-r(2n-1,t)]^{3/2}_{+}-\frac{1}{\eta^{2}}[\delta-r(2n+1,t)]^{3/2}_{+}\,,$
$\displaystyle\ddot{r}(2n+1,t)$
$\displaystyle=\bigg{(}1+\frac{1}{\eta^{2}}\bigg{)}[\delta-r(2n+1,t)]^{3/2}_{+}-\frac{1}{\eta^{2}}[\delta-r(2n,t)]^{3/2}_{+}-[\delta-r(2n+2,t)]^{3/2}_{+}\,.$

As with the woordpile chain, the domain includes $M=2^{10}$ particles, with
indices $n\in\\{-M/2+1,\ldots,M/2\\}$. At each time step, we multiply the
solution by the windowing function in (54). The initial condition is given by
the leading-order solution (62).

We calculate the amplitude of the far-field waves for precompressions of
$\delta=5$, $\delta=7.5$, $\delta=10$, and $\delta=15$ with long-wavelength
parameters of $\epsilon=0.4$, $\epsilon=0.5$, and $\epsilon=0.6$ for a range
of values of the mass-ratio parameter $\eta$. As with woodpile chains, we run
each simulation for a sufficiently long time so that the amplitude of the far-
field oscillations appears to reach a constant value. In Fig. 8, we present
the results of our computations and compare them to our asymptotic results
from (91).

(a) $\delta=3$ (b) $\delta=5$ (c) $\delta=10$ (d) $\delta=15$

Figure 8: Comparison of our asymptotic results and numerical simulations of
the far-field wave amplitude of nanoptera in a diatomic Hertzian chain with
precompressions of (a) $\delta=3$, (b) $\delta=5$, (c) $\delta=10$, and (d)
$\delta=15$. For each value of $\delta$, we show our computations for long-
wavelength parameters of $\epsilon=0.4$, $\epsilon=0.5$, and $\epsilon=0.6$.
Our asymptotic approximation of the far-field amplitude is accurate for a wide
range of values of the mass-ratio parameter $\eta$, with certain important
exceptions. For fixed values of $\delta$ and $\epsilon$, two oscillatory wave
trains interfere destructively for specific values of $\eta$. Our asymptotic
approximation is capable of predicting the destructive interference that
occurs for $\eta\approx 0.32$ and $\eta\approx 0.22$, and the accuracy of our
approximation becomes progressively better for progressively smaller values of
$\eta$. However, the approximation does not predict the destructive
interference that occurs for $\eta\approx 0.25$ and $\eta\approx 0.19$. Our
approximation for the anti-resonance conditions becomes less accurate for
progressively smaller $\epsilon$.

We make two observations from Fig. 8. First, our asymptotic approximation
accurately predicts the amplitude of the far-field waves for a wide range of
values of $\eta$. Second, our asymptotic calculation does not always identify
the values of $\eta$ for which destructive interference causes the far-field
oscillations to vanish. In each panel of Fig. 8, we observe four values of
$\eta$ (these are $\eta\approx 0.32$, $\eta\approx 0.25$, $\eta\approx 0.22$,
and $\eta\approx 0.19$) at which the far-field waves in our numerical
computations vanish, corresponding to localized solitary waves with no
exponentially small oscillations behind the leading-order solitary wave. From
(91), our asymptotic analysis predicts these values to be solutions of

$\cos\bigg{(}\frac{\mathrm{Im}(\chi_{1-}(\xi)-\chi_{1+}(\xi))}{2\eta}+\frac{\pi}{4}\bigg{)}=0\,,$
(93)

which simplifies to solving

$\displaystyle
2\pi\eta\left(K+\frac{1}{4}\right)=\frac{\sqrt{3}\mathrm{i}}{c_{\epsilon}}\int_{\xi_{1-}}^{\xi_{1+}}\bigg{(}\delta+\frac{(u_{0}(\xi-1)-u_{0}(\xi+1))}{2}\bigg{)}^{1/4}\,\mathrm{d}s$
(94)

for integer values of $K$. The expression (94) is an orthogonality condition
[57, 56, 93, 44, 1]. Values of $\eta$ that satisfy (94) correspond to dips in
the asymptotically predicted amplitude in Fig. 8. These values yield localized
solitary waves. We see that the condition (94) is capable of predicting the
localized solitary waves that arise for $\eta\approx 0.32$ and $\eta\approx
0.22$. However, it fails to predict the localized solitary waves that arise
for $\eta\approx 0.25$ and $\eta\approx 0.19$. In fact, our exponential
asymptotic analysis appears to systematically miss every second point (i.e.,
if it captures one point, then it misses the next point, and vice versa) at
which the oscillations vanish. Additionally, for progressively smaller values
of $\epsilon$, the orthogonality condition (94) becomes progressively less
accurate at predicting the values of $\eta$ that give localized solitary
waves.

Systematically missing half of these solitary-wave solutions appears to be a
limitation of the method that we used to obtain our asymptotic expressions.
There are several approximations in our analysis. Notably, we analytically
continued the long-wavelength KdV approximation from the real axis into the
complex plane. The orthogonality condition (94) depends on the behavior of
$u_{0}(\xi)$ on a contour that connects the singularities at $\xi_{1-}$ and
$\xi_{1+}$ (see Fig. 6). These singularities are far from the real axis, and
their distance from the real axis increases as the long-wavelength parameter
$\epsilon$ decreases. Consequently, it is possible that errors in this
approximation — and, in particular, in the locations of the singular points —
become progressively larger for progressively smaller $\epsilon$.

Detecting destructive interference involves determining configurations for
which two oscillations with very short wavelengths precisely cancel each
other. The individual Stokes contributions have wavelengths that are of size
$\mathcal{O}(\eta)$ as $\eta\rightarrow 0$. Errors in the positions of the
singularities on this very short scale are capable of significantly changing
the values of $\eta$ for which cancellation occurs. It is possible that
approximation errors in the leading-order solution lead to shifts in the
singularities in the analytically-continued solution. This may disrupt the
precise cancellation that is necessary for the manifestation of localized
solitary-wave solutions, causing our asymptotic approximation to miss some of
them.

The issue of missing solitary-wave solutions did not arise in a recent study
of a diatomic FPUT chain [56], despite the use of a similar long-wavelength
approximation to continue a leading-order solution into the complex plane. It
is not apparent why the method worked effectively for diatomic FPUT chains but
has failed to identify analogous behavior in our diatomic Hertzian chain. This
is an interesting open question that merits further study.

##  5 Conclusions and Discussion

In the present study, we derived asymptotic approximations for traveling waves
in two particle chains that are singular perturbations of a monoatomic
particle chain with precompression. Specifically, we considered a woodpile
chain and a diatomic Hertzian chain with a small ratio between the masses of
its two types of particles. The traveling-wave solutions of these systems are
nanoptera, which consist of solitary wave along with non-decaying oscillations
of exponentially small amplitude.

A monoatomic Hertzian chain supports traveling solitary-wave solutions, which
one can approximate as soliton solutions of the KdV equation under suitable
assumptions. We used this monoatomic solitary wave as the leading-order
behavior in both woodpile chains and diatomic Hertzian chains. In both
systems, we found that non-decaying oscillations appear in the wake of a
primary traveling wave. These oscillations arise from the Stokes phenomenon,
in which exponentially small contributions to an analytically-continued
solution switch on or off as Stokes curves are crossed in the complex plane.

By employing exponential asymptotic analysis, we obtained an asymptotic form
for these exponentially small contributions, and we used this form to
determine the amplitude of the non-decaying far-field waves that appear behind
the leading-order wave. For woodpile chains, the solution has one Stokes
curve, which produces one exponentially small oscillation in the region behind
the wave front. We determined an asymptotic approximation of this behavior,
and we compared our asymptotic results with numerical computations. We saw
that the asymptotic prediction accurately approximates the wave behavior in
the far field. However, the asymptotic approximation of the far-field wave
amplitude becomes progressively less accurate for progressively smaller values
of the precompression parameter $\delta$. It is possible that this is a
consequence of the inadequacy of examining a weakly nonlinear regime of the
woodpile chain as nonlinear effects become more significant, leading to errors
in the description of exponentially small effects. It would be interesting to
study woodpile chains without precompression to investigate whether it is
possible to analyze the behavior of the far-field oscillations in the strongly
nonlinear regime of woodpile systems using exponential asymptotic methods. We
discuss the feasibility 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 which the mass ratio between light and heavy particles is
small. We found that there are two Stokes curves in the analytically-continued
traveling-wave solution. These produce two exponentially small oscillations —
which have identical amplitude but different phases — behind the wave front.
There exists a set of values of the mass-ratio parameter $\eta$ at which the
oscillations interfere destructively to produce a localized solitary wave.

Our comparison of the asymptotic results and numerical computations in
diatomic Hertzian chains shows that our exponential asymptotic analysis
accurately approximates the amplitude of the far-field waves for a wide range
of values of $\eta$. However, our computations also reveal that our asymptotic
analysis is unable to detect the wave cancellation that occurs for several
values of $\eta$. Furthermore, the predicted mass ratio values that produce
wave cancellation are less accurate than those in analyses of other particle
chains (such as the exponential asymptotic calculations from [57, 56]),
particularly for small values of the long-wavelength parameter $\epsilon$. It
is possible that this inaccuracy arises from our use of a long-wavelength
approximation for the leading-order solution; this approximation may deviate
from the exact solution away from the real axis. This is an important
question, as the Stokes phenomenon depends on the behavior of singularities in
the analytically-continued solution. It is worthwhile to study whether using
an alternative method of approximating the leading-order solution can capture
all of the values of $\eta$ that cause the oscillations to cancel.

Our work also opens the following interesting questions. In the present study,
we used a similar asymptotic approach as the one in [57, 56] and evaluated our
asymptotic approximations by comparing them to numerical computations. It
would be interesting to develop a rigorous existence proof of nanopteron
solutions in both diatomic Hertzian chains and woodpile chains. One possible
approach for attempting this is by adapting the Beale-ansatz method of [26,
40].

As we discussed in Section 1, one-sided nanopteron solutions in both woodpile
and diatomic Hertzian chains are metastable because the oscillations slowly
draw energy away from the wave front. This indicates that the wave cannot
persist indefinitely; instead, it must eventually decay. This property is not
present in the leading-order approximation, nor it is apparent in the form of
the exponentially small oscillations in the analysis either in the present
study or in [57, 56]. The long-time decay must arise from interactions between
the 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 or not one can
detect such interactions by using an appropriate long-time rescaling or by
calculating higher-order corrections to the oscillatory wave train.

In our analysis, we considered a weakly nonlinear regime for both woodpile and
diatomic Hertzian chains. In principle, one can apply the exponential
asymptotic method that was developed in [17] with a few additional
complications to study strongly nonlinear systems, such as a singularly
perturbed Hertzian system without precompression. However, there are two
significant challenges that need to be overcome to undertake such an analysis.
The first challenge is that one no longer has access to a closed-form leading-
order solution to analytically continue in a straightforward fashion. Perhaps
it is possible to use the leading-order approximation that was derived in
[83], but it is not apparent whether this approximation is valid if it is
analytically continued away from the real axis. This is an important point to
clarify because the Stokes phenomenon is caused by singularities in the
complex plane. The second challenge relates to the form of the potential in
(2). This potential is not smooth at any point at which the argument in the
brackets is $0$. In the present study, the precompression parameter
$\delta_{0}$ is larger than or equal to the relative displacement of two
adjacent particles, so the argument in the brackets is always non-negative and
the potential in (2) is smooth. Without precompression, however, the argument
in the brackets can switch between positive and negative signs, and the
potential is not smooth everywhere. An important direction for future studies
is to explore the effects of such a non-smooth potential on exponential
asymptotic analysis. If these challenges are overcome, it would be useful to
compare the ensuing results with those of the uncompressed woodpile [98] and
diatomic Hertzian [44] chains.

##  6 Acknowledgement

CJL and GD were supported by Australian Research Council Discovery Project
DP190101190.

##  Appendix A Determining the Prefactor Constants

We determine the prefactor constants for both the woodpile chain and the
diatomic Hertzian chain. Assuming that the power series diverges in a
factorial-over-power fashion, it must cease to be asymptotic in a narrow
region in the neighborhood of the singularities in the complex plane. Inside
this region, the earlier terms in the series are not larger asymptotically
than the later series terms in the limit $\eta\rightarrow 0$. Describing the
solution behavior near one of these singularities requires obtaining a local
expansion of the solution near the singular point. Once we have obtained this
local expansion, we calculate the prefactor constants in the late-order terms
using asymptotic matching to ensure that the local expansions are consistent
with the power-series behavior.

###  A.1 Woodpile Chain

To determine the value of $\Lambda_{+}$, we match the late-order expansion in
the outer region with the local solution in an inner region near the
singularity at $\xi=\xi_{0+}$ by using Van Dyke’s matching principle.

As $\xi\rightarrow\xi_{0+}$, we find that

$\displaystyle u_{0}(\xi)$
$\displaystyle\sim-\frac{1}{\alpha-1}\frac{\delta}{\xi-\xi_{0+}}+\mathcal{O}(\xi-\xi_{0+})\,,\quad$
$\displaystyle u_{0}(\xi+1)$
$\displaystyle\sim\frac{1}{\alpha-1}\delta\epsilon\,\mathrm{coth}(\epsilon)+\mathcal{O}(\xi-\xi_{0+})\,,$
(95) $\displaystyle v_{0}(\xi)$
$\displaystyle\sim-\frac{1}{\alpha-1}\frac{\delta}{\xi-\xi_{0+}}+\mathcal{O}(\xi-\xi_{0+})\,,\quad$
$\displaystyle u_{0}(\xi-1)$
$\displaystyle\sim-\frac{1}{\alpha-1}\delta\epsilon\,\mathrm{coth}(\epsilon)+\mathcal{O}(\xi-\xi_{0+})\,.$
(96)

To locate the relevant inner region, we need to determine where the validity
of the late-order term ansatz breaks down. From the form of the late-order
ansatz (28), we see that this occurs for $\eta^{2}\chi^{-2}=\mathcal{O}(1)$ as
$\eta\rightarrow 0$. That is, it occurs for
$\eta^{2}(\xi-\xi_{0+})^{-2}=\mathcal{O}(1)$; this corresponds to the inner
scaling $\xi-\xi_{0+}=\eta\overline{\xi}$. From asymptotic balancing, the
appropriate rescaled inner variables are

$u(\xi)=-\frac{1}{\alpha-1}\frac{\delta}{\eta\overline{\xi}}+\hat{u}(\overline{\xi})\,,\quad
u(\xi+1)=\hat{u}(\overline{\xi}+\eta^{-1})\,,\quad
u(\xi-1)=\hat{u}(\overline{\xi}-\eta^{-1})\,,\quad
v(\xi)=-\frac{1}{\alpha-1}\frac{\delta}{\eta\overline{\xi}}+\frac{\hat{v}(\overline{\xi})}{\eta}\,.$
(97)

Retaining the leading-order terms as $\eta\rightarrow 0$, the rescaled inner
equation gives

$-\frac{1}{\alpha-1}\frac{2\delta}{\overline{\xi}^{3}}+\frac{\mathrm{d}^{2}\hat{v}(\overline{\xi})}{\mathrm{d}\overline{\xi}^{2}}=-\frac{k}{c_{\epsilon}^{2}}\hat{v}(\overline{\xi})\,.$
(98)

We express $\hat{v}$ as a power series

$\hat{v}(\overline{\xi})\sim\sum_{j=1}^{\infty}\frac{v_{n}}{\overline{\xi}^{2j+1}}\quad\mathrm{as}\quad\overline{\xi}\rightarrow
0\,,$ (99)

and we note that we include the leading-order singularity as part of the
rescaling process (97). Substituting (99) into (98) yields

$-\frac{1}{\alpha-1}\frac{2\delta}{\overline{\xi}^{3}}+\sum_{j=1}^{\infty}\frac{(2j+1)(2j+2)v_{j}}{\overline{\xi}^{2j+3}}=-\frac{k}{c_{\epsilon}^{2}}\sum_{j=1}^{\infty}\frac{v_{j}}{\overline{\xi}^{2j+1}}\,.$
(100)

By matching orders of $\overline{\xi}$, we obtain the recurrence relation

$v_{1}=2\delta
c_{\epsilon}^{2}/(k(\alpha-1))\,,\quad(2j+2)(2j+1)c_{\epsilon}^{2}v_{j}=-kv_{j+1}\,.$
(101)

Solving the recurrence relation (101) gives

$v_{j}=\frac{1}{\alpha-1}\delta\left(-1\right)^{j+1}\left(\frac{c_{\epsilon}^{2}}{k}\right)^{j}\Gamma(2j+1)\,.$
(102)

By comparing the series expression (102) with the inner limit of the late-
order ansatz, we obtain

$\Lambda_{+}=\lim_{j\rightarrow\infty}\frac{v_{j}(\mathrm{i}\sqrt{k}/c_{\epsilon})^{2j+1}}{\Gamma(2j+1)}=-\frac{1}{\alpha-1}\frac{\mathrm{i}\delta\sqrt{k}}{c_{\epsilon}}\,.$
(103)

###  A.2 Diatomic Hertzian Chain

To determine $\Lambda_{1-}$, we match the late-order term expansion in the
outer region with a local expansion in an inner region near $\xi=\xi_{1-}$. We
perform this analysis in the neighborhood of $\xi=\xi_{1-}$. From (68), we see
that the factorial-over-power ansatz breaks down when
$\eta^{2}\chi^{-2}=\mathcal{O}(1)$ as $\eta\to 0$. We introduce an inner
scaling $\xi-\xi_{1-}=\eta^{4/3}\bar{\xi}$, which gives the rescaled inner
variables

$\displaystyle
u(\xi+1)=-\frac{8\delta}{\eta^{4/3}\bar{\xi}}+\hat{u}(\bar{\xi}+\eta^{-4/3})\,,\quad
u(\xi-1)=\hat{u}(\bar{\xi}-\eta^{-4/3})\,,\quad
v(\xi)=-\frac{4\delta}{\eta^{4/3}\bar{\xi}}+\frac{\hat{v}(\bar{\xi})}{\eta^{4/3}}\,.$
(104)

Retaining the leading-order terms as $\eta\to 0$, the rescaled inner equation
gives

$\displaystyle-\frac{8\delta}{\bar{\xi}^{3}}+\frac{\,\mathrm{d}^{2}\hat{v}(\bar{\xi})}{\,\mathrm{d}{\bar{\xi}}^{2}}=\frac{1}{c_{\epsilon}^{2}}\left[\bigg{(}\frac{4\delta}{\bar{\xi}}-\hat{v}(\xi)\bigg{)}^{3/2}-\bigg{(}\frac{4\delta}{\bar{\xi}}+\hat{v}(\xi)\bigg{)}^{3/2}\right]\,.$
(105)

We express $\hat{v}$ in terms of a power series

$\displaystyle\hat{v}(\bar{\xi})\sim\sum_{j=1}^{\infty}\frac{a_{j}}{\bar{\xi}^{\frac{3j}{2}+1}}\quad\mathrm{as}\quad\bar{\xi}\to
0\,.$ (106)

Matching the inner expansion (106) with the outer ansatz (78) yields

$\displaystyle\Lambda_{1-}=\lim_{j\to\infty}\frac{a_{j}\mathrm{i}^{2j+2}2^{5j+4}\delta^{j/2+1/2}}{3^{j+1/2}c_{\epsilon}^{2j+2}\Gamma(2j+3/2)}\,.$
(107)

We numerically compute values for $a_{j}$ by substituting the series
expression (106) into (105) and solving a recurrence relation for $a_{j}$ in
terms of the previous series coefficients. We obtain an approximation
$\Lambda_{\mathrm{approx}}$ for $\Lambda_{1-}$ by evaluating the right-hand
side of (107) for large values of $j$. When we take $j$ to be sufficiently
large, we obtain an accurate approximation of the exact value of
$\Lambda_{1-}$.

In Fig. 9, we show the behavior of $\Lambda_{\mathrm{approx}}$ as we increase
$j$; we observe that our approximation is converging. For sufficiently large
values of $j$, we compute that

$\Lambda_{1-}=38.41\,\delta^{3/2}/c_{\epsilon}^{2}\,.$ (108)

Figure 9: Determining $\Lambda_{\mathrm{approx}}$ using (107) for
progressively larger values of $j$. The dashed line indicates the value to
which the rescaled $\Lambda_{\mathrm{approx}}$ appears to converge as we
increase $j$.

Following the same approach as above yields

$\Lambda_{1+}=\mathrm{i}\Lambda_{1-}\,.$ (109)

## References

  * [1] G. L. Alfimov, A. S. Korobeinikov, C. J. Lustri, and D. E. Pelinovsky.  Standing lattice solitons in the discrete 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. Ávalos and S. Sen.  How solitary waves collide in discrete granular alignments.  Phys. Rev. E, 79:046607, 2009.
  * [4] E. Ávalos and S. Sen.  Granular chain between asymmetric boundaries and the quasiequilibrium state.  Phys. Rev. E, 89:053202, 2014.
  * [5] E. Ávalos, D. Sun, R. L. Doney, and S. Sen.  Sustained strong fluctuations in a nonlinear chain at acoustic vacuum: 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ÉS, 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, The Netherlands, 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 in one-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 $\phi^{4}$ breather.  Nonlinearity, 3(1):177–195, 1990.
  * [13] J. P. Boyd.  Weakly non-local solitons for capillary-gravity waves: Fifth-degree Korteweg–de Vries equation.  Physica D, 48:129–146, 1991.
  * [14] J. P. Boyd.  Weakly Nonlocal Solitary Waves and Beyond-All-Orders Asymptotics: Generalized Solitons and Hyperasymptotic Perturbation Theory, volume 442 of Mathematics and Its Applications.  Kluwer Publishers, 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 nonlinear ordinary 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 granular crystals.  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 a composite 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, S. Sen, and P. G. Kevrekidis.  On the generation and propagation of solitary waves in integrable and nonintegrable 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ët 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 diatomic Fermi–Pasta–Ulam–Tsingou problem.  arXiv preprint arXiv:2010.02107, 2020.
  * [28] F. Fraternali, G. Carpentieri, A. Amendola, R. E. Skelton, and V. F. Nesterenko.  Multiscale tunability of solitary wave dynamics in tensegrity metamaterials.  Appl. Phys. Lett., 105(20):201903, 2014.
  * [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 FPUT lattices.  Math. Eng., 1:419–433, 2019.
  * [32] J. D. Goddard.  Nonlinear elasticity and pressure-dependent wave speeds in granular media.  Proc. Roy. Soc. A, 430(1878):105–131, 1990.
  * [33] R. Grimshaw and N. Joshi.  Weakly nonlocal solitary waves in a singularly perturbed Korteweg–de Vries equation.  SIAM J. Appl. Math., 55(1):124–135, 1995.
  * [34] U. Harbola, A. Rosas, A. H. Romero, and K. Lindenberg.  Pulse propagation in randomly decorated chains.  Phys. Rev. E, 82:011306, 2010.
  * [35] E. B. Herbold and V. F.Nesterenko.  Propagation of rarefaction pulses in discrete materials with strain-softening behavior.  Phys. Rev. Lett., 110:144101, 2013.
  * [36] E. B. Herbold, J. Kim, V. F. Nesterenko, S. Y. Wang, and C. Daraio.  Pulse propagation in a linear and nonlinear diatomic periodic chain: Effects of acoustic frequency band-gap.  Acta Mech., 205:85–103, 2009.
  * [37] H. Hertz.  Üeber die berührung fester elastischer körper.  J. Reine Angew. Math., 92:156–171, 1881.
  * [38] E. J. Hinch.  Perturbation Methods.  Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, UK, 1991.
  * [39] E. J. Hinch and S. Saint-Jean.  The fragmentation of a line of balls by an impact.  Proc. Roy. Soc. A, 455:3201–3220, 1999.
  * [40] A. Hoffman and J. D. Wright.  Nanopteron solutions of diatomic Fermi–Pasta–Ulam–Tsingou lattices with small mass-ratio.  Physica D, 358:33–59, 2017.
  * [41] 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 diatomic granular crystals.  Europhys. Lett., 101(4):44003, 2013.
  * [42] J. K. Hunter and J. Scheurle.  Existence of perturbed solitary wave solutions to a model equation for water waves.  Physica D, 32(2):253–268, 1988.
  * [43] G. Iooss and K. Kirchgässner.  Travelling waves in a chain of coupled nonlinear oscillators.  Commun. Math. Phys., 211(2):439–464, 2000.
  * [44] K. R. Jayaprakash, Y. Starosvetsky, and A. F. Vakakis.  New family of solitary waves in granular dimer chains with no precompression.  Phys. Rev. E, 83(3):036606, 2011.
  * [45] K. R. Jayaprakash, Y. Starosvetsky, A. F. Vakakis, and O. V. Gendelman.  Nonlinear resonances leading to strong pulse attenuation in granular dimer chains.  J. Nonlinear Sci., 23(3):363–392, 2013.
  * [46] 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.
  * [47] S. Job, F. Melo, A. Sokolow, and S. Sen.  How Hertzian solitary waves interact with boundaries in a 1D granular medium.  Phys. Rev. Lett., 94:178002, 2005.
  * [48] E. Kim, R. Chaunsali, H. Xu, J. Jaworski, J. Yang, P. G. Kevrekidis, and A. F. Vakakis.  Nonlinear low-to-high-frequency energy cascades in diatomic granular crystals.  Phys. Rev. E, 92:062201, 2015.
  * [49] E. Kim, F. Li, C. Chong, G. Theocharis, J. Yang, and P. G. Kevrekidis.  Highly nonlinear wave propagation in elastic woodpile periodic structures.  Phys. Rev. Lett., 114(11):118002, 2015.
  * [50] E. Kim, A. J. Martinez, S. E. Phenisee, P.G. Kevrekidis, M. A. Porter, and J. Yang.  Direct measurement of superdiffusive energy transport in disordered granular chain.  Nat. Commun., 9:640, 2018.
  * [51] E. Kim and J. Yang.  Wave propagation in single column woodpile phononic crystals: Formation of tunable band gaps.  J. Mech. Phys. Solids, 71:33–45, 2014.
  * [52] D. J. Korteweg and G. de Vries.  XLI. On the change of form of long waves advancing in a rectangular canal, and on a new type of long stationary waves.  Lond. Edinb. Dubl. Phil. Mag., 39(240):422–443, 1895.
  * [53] A. N. Lazaridi and V. F. Nesterenko.  Observation of a new type of solitary waves in a one-dimensional granular medium.  J. Appl. Mech. Tech. Phy., 26:405–408, 1985.
  * [54] L. Liu, G. James, P. Kevrekidis, and A. Vainchtein.  Breathers in a locally resonant granular chain with precompression.  Physica D, 331:27–47, 2016.
  * [55] L. Liu, G. James, P. Kevrekidis, and A. Vainchtein.  Strongly nonlinear waves in locally resonant granular chains.  Nonlinearity, 29(11):3496–3527, 2016.
  * [56] C. J. Lustri.  Nanoptera and Stokes curves in the 2-periodic Fermi–Pasta–Ulam–Tsingou equation.  Physica D, 402:132239, 2020.
  * [57] C. J. Lustri and M. A. Porter.  Nanoptera in a period-2 Toda chain.  SIAM J. Appl. Dyn. Syst., 17(2):1182–1212, 2018.
  * [58] F. S. Manciu and S. Sen.  Secondary solitary wave formation in systems with generalized Hertz interactions.  Phys. Rev. E, 66:016616, 2002.
  * [59] M. Manciu, S. Sen, and A. J. Hurd.  Impulse propagation in dissipative and disordered chains with power-law repulsive potentials.  Physica D, 157(3):226–240, 2001.
  * [60] M. Manjunath, A. P. Awasthi, and P. H. Geubelle.  Wave propagation in random granular chains.  Phys. Rev. E, 85:031308, 2012.
  * [61] M. Manjunath, A. P. Awasthi, and P. H. Geubelle.  Family of plane solitary waves in dimer granular crystals.  Phys. Rev. E, 90:032209, 2014.
  * [62] A. J. Martínez, P. G. Kevrekidis, and M. A. Porter.  Superdiffusive transport and energy localization in disordered granular crystals.  Phys. Rev. E, 93:022902, 2016.
  * [63] A. J. Martínez, M. A. Porter, and P. G. Kevrekidis.  Quasiperiodic granular chains and hofstadter butterflies.  Phil. Trans. Roy. Soc. A, 376:20170139, 2018.
  * [64] A. J. Martínez, H. Yasuda, E. Kim, P. G. Kevrekidis, M. A. Porter, and J. Yang.  Scattering of waves by impurities in precompressed granular chains.  Phys. Rev. E, 93:052224, 2016.
  * [65] A. Molinari and C. Daraio.  Stationary shocks in periodic highly nonlinear granular chains.  Phys. Rev. E, 80:056602, 2009.
  * [66] V. F. Nesterenko.  Propagation of nonlinear compression pulses in granular media.  J. Appl. Mech. Tech. Phy., 24:733–743, 1983.
  * [67] V. F. Nesterenko.  High Rate Deformation of Heterogeneous Materials.  Nauka, Novosibirsk, Russia, 1992.
  * [68] V. F. Nesterenko.  Nonlinear waves in “sonic vacuum”.  Fiz. Goreniya Vzryva, 28:121–123, 1992.
  * [69] V. F. Nesterenko.  Examples of “sonic vacuum”.  Fiz. Goreniya Vzryva, 29:132–134, 1993.
  * [70] V. F. Nesterenko.  Solitary waves in discrete medium with anomalous compressibility.  Fiz. Goreniya Vzryva, 29:134–136, 1993.
  * [71] V. F. Nesterenko.  Dynamics of Heterogeneous Materials.  Springer-Verlag, Heidelberg, Germany, 2001.
  * [72] V. F. Nesterenko.  Waves in strongly nonlinear discrete systems.  Phil. Trans. Roy. Soc. A, 376:20170130, 2018.
  * [73] V. F. Nesterenko, C. Daraio, E. B. Herbold, and S. Jin.  Anomalous wave reflection at the interface of two strongly nonlinear granular media.  Phys. Rev. Lett., 95(15):158702, 2005.
  * [74] Y. Okada, S. Watanabe, and H. Tanaca.  Solitary wave in periodic nonlinear lattice.  J. Phys. Soc. Jpn., 59(8):2647–2658, 1990.
  * [75] A. B. Olde Daalhuis, S. J. Chapman, J. R. King, J. R. Ockendon, and R. H. Tew.  Stokes phenomenon and matched asymptotic expansions.  SIAM J. Appl. Math., 55(6):1469–1483, 1995.
  * [76] L. Ponson, N. Boechler, Y. M. Lai, M. A. Porter, P. G. Kevrekidis, and C. Daraio.  Nonlinear waves in disordered diatomic granular chains.  Phys. Rev. E, 82:021301, 2010.
  * [77] M. A. Porter, C. Daraio, E. B. Herbold, I. Szelengowicz, and P. G. Kevrekidis.  Highly nonlinear solitary waves in periodic dimer granular chains.  Phys. Rev. E, 77:015601, 2008.
  * [78] M. A. Porter, C. Daraio, I. Szelengowicz, E. B. Herbold, and P. G. Kevrekidis.  Highly nonlinear solitary waves in heterogeneous periodic granular media.  Physica D, 238(6):666–676, 2009.
  * [79] M. A. Porter, P.G. Kevrekidis, and C. Daraio.  Granular crystals: Nonlinear dynamics meets materials engineering.  Phys. Today, 68(11):44–50, 2015.
  * [80] 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.
  * [81] M. Przedborski, T. A. Harroun, and S. Sen.  Granular chains with soft boundaries: Slowing the transition to quasiequilibrium.  Phys. Rev. E, 91:042207, 2015.
  * [82] M. Przedborski, S. Sen, and T. A. Harroun.  The equilibrium phase in heterogeneous Hertzian chains.  J. Stat. Mech., 2017(12):123204, 2017.
  * [83] S. Sen, J. Hong, J. Bang, E. Ávalos, and R. Doney.  Solitary waves in the granular chain.  Phys. Rep., 462(2):21–66, 2008.
  * [84] S. Sen, T. R. Krishna Mohan, and J. M. M. Pfannes.  The quasi-equilibrium phase in nonlinear 1D systems.  Physica A, 342(1):336–343, 2004.
  * [85] S. Sen and M. Manciu.  Solitary wave dynamics in generalized hertz chains: An improved solution of the equation of motion.  Phys. Rev. E, 64:056605, 2001.
  * [86] S. Sen, M. Manciu, and J. D. Wright.  Solitonlike pulses in perturbed and driven Hertzian chains and their possible applications in detecting buried impurities.  Phys. Rev. E, 57:2386–2397, 1998.
  * [87] Y. Shen, P. G. Kevrekidis, S. Sen, and A. Hoffman.  Characterizing traveling-wave collisions in granular chains starting from integrable limits: The case of the Korteweg–de Vries equation and the Toda lattice.  Phys. Rev. E, 90:022905, 2014.
  * [88] D. A. Spence.  Self similar solutions to adhesive contact problems with incremental loading.  Proc. Roy. Soc. A, 305(1480):55–80, 1968.
  * [89] G. G. Stokes.  On the discontinuity of arbitrary constants which appear in divergent developments.  Trans. Cam. Phil. Soc., 10:106–128, 1864.
  * [90] D. Sun, C. Daraio, and S. Sen.  Nonlinear repulsive force between two solids with axial symmetry.  Phys. Rev. E, 83:066605, 2011.
  * [91] Y. Tabata.  Stable solitary wave in diatomic Toda lattice.  J. Phys. Soc. Jpn., 65(12):3689–3691, 1996.
  * [92] G. Theocharis, N. Boechler, P. G. Kevrekidis, S. Job, M. A. Porter, and C. Daraio.  Intrinsic energy localization through discrete gap breathers in one-dimensional diatomic granular crystals.  Phys. Rev. E, 82:056604, 2010.
  * [93] A. Vainchtein, Y. Starosvetsky, J. D. Wright, and R. Perline.  Solitary waves in diatomic chains.  Phys. Rev. E, 93(4):042210, 2016.
  * [94] L. Vergara.  Scattering of solitary waves from interfaces in granular media.  Phys. Rev. Lett., 95:108002, 2005.
  * [95] L. Vergara.  Delayed scattering of solitary waves from interfaces in a granular container.  Phys. Rev. E, 73:066623, 2006.
  * [96] L. Verlet.  Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules.  Phys. Rev., 159:98–103, 1967.
  * [97] 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.
  * [98] H. Xu, P. G. Kevrekidis, and A. Stefanov.  Traveling waves and their tails in locally resonant granular systems.  J. Phys. A Math. Theor., 48(19):195204, 2015.
  * [99] H. Yasuda, C. Chong, J. Yang, and P. G. Kevrekidis.  Emergence of dispersive shocks and rarefaction waves in power-law contact models.  Phys. Rev. E, 95:062216, 2017.