# Exponential asymptotics of woodpile chain nanoptera using numerical analytic
continuation

Guo Deng111Electronic address<EMAIL_ADDRESS>Department of Mathematics
and Statistics, 12 Wally’s Walk, Macquarie University, New South Wales 2109,
Australia  Christopher J. Lustri222Corresponding Author. Electronic address:
<EMAIL_ADDRESS>Department of Mathematics and Statistics, 12
Wally’s Walk, Macquarie University, New South Wales 2109, Australia

###### Abstract

Travelling waves in woodpile chains are typically nanoptera, which are
composed of a central solitary wave and exponentially small oscillations.
These oscillations have been studied using exponential asymptotic methods,
which typically require an explicit form for the leading-order behaviour. For
many nonlinear systems, such as granular woodpile chains, it is not possible
to calculate the leading-order solution explicitly. We show that accurate
asymptotic approximations can be obtained using numerical approximation in
place of the exact leading-order behaviour. We calculate the oscillation
behaviour for Toda woodpile chains, and compare the results to exponential
asymptotics based on tanh-fitting, Padé approximants, and the adaptive
Antoulas-Anderson (AAA) method. The AAA method is shown to produce the most
accurate predictions of the amplitude of the oscillations and the mass ratios
for which the oscillations vanish. This method is then applied to study
granular woodpile chains, including chains with Hertzian interactions – this
method is able to calculate behaviour that could not be accurately
approximated in previous studies.

##  1 Introduction

###  1.1 Exponential asymptotics

Exponential asymptotics are techniques used to study behaviour in
exponentially small behaviour in singularly-perturbed systems. Important
behaviour on exponentially small scales arises in models from quantum
mechanics [57, 45, 59, 35, 40], optics [80, 21, 48, 58, 63], crystal growth
[16, 47], fluid mechanics [19, 52, 53, 20, 54], particle chains [55, 51, 27,
28] and cosmology [2].

Classic asymptotic power series cannot describe these exponentially small
terms, as they are smaller than any algebraic power of the asymptotic
parameter. Exponential asymptotic methods, first developed in [11, 12], make
use of the idea that truncating an asymptotic power series optimally produces
an exponentially small remainder. This remainder term can be studied in order
to calculate the asymptotic behaviour of the solution on exponentially small
scales. The method used in the paper, developed in [17, 64], are described in
more detail in Section 2.

Exponentially small asymptotic contributions in singularly-perturbed problems
typically exhibit behaviour known as Stokes’ phenomenon. This describes rapid
changes in the exponentially small contribution that occur as curves in the
complex plane, known as Stokes curves, are crossed.

Exponential asymptotic techniques typically require the explicit calculation
of the leading-order solution in the asymptotic limit. Stokes curves originate
at singularities in the leading-order solution, and this solution also affects
the form of the exponentially-small asymptotic behaviour that is switched as
these curves are crossed. For many nonlinear problems, such as [51, 27, 28,
18], it is impossible to calculate the leading-order behaviour explicitly. The
leading-order behaviour can be calculated numerically, but this numerical
solution must be analytically continued to determine the location of singular
points in the complex plane.

An example of a physical system for which the leading-order solution cannot be
calculated exactly is the granular woodpile chain studied in [28]. A tanh-
fitting method based on the work of [69] was used to approximate the leading-
order behaviour. This method produced accurate results for large values of the
interaction exponent but failed to capture important behaviour for smaller
values, including the important case of Hertzian particle interactions.

In this paper we investigate how numerical approximations to the leading-order
behaviour affect the accuracy of using exponential asymptotic methods to study
woodpile chains. We perform an exponential asymptotic analysis on nonlinear
waves in a Toda woodpile lattice, which has an exact leading-order solution.
We then compare three possible approximation methods for determining the
leading-order behaviour in the complex plane against the results of the exact
analysis: tanh-fitting [69], Padé approximants [7], and the adaptive Antoulas-
Anderson (AAA) method [60].

We then apply the most effective of these methods (AAA approximation) to
produce descriptions of exponentially small behaviour in granular woodpile
chains, including the Hertzian woodpile chains that could not be accurately
calculated in [28]. We show that this inaccuracy occured because tanh-fitting
cannot accurately identify the location of subdominant singularities. This
problem is avoided using a method based on AAA approximation for the leading-
order behaviour, which is able to accurately identify subdominant singularity
locations and describe their interaction effects.

###  1.2 Woodpile chains

The motion of homogeneous particle chains, where every particle is identical,
is governed by the system of differential–difference equations

$m\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}$, dots denote temporal derivatives, $\phi$ is the
interaction potential and $\phi^{\prime}$ is the derivative of $\phi$. It is
known that solitary wave solutions exist in lattices with superquadratic
potentials [34, 70]. Among lattices with superquadratic potentials, the Toda
lattice [74, 71, 72, 73, 75] and the granular chain [61, 62, 67] have been of
particular interest due to their theoretical and practical importance.

The Toda lattice describes particle interactions with interaction potential

$\phi(r)=\frac{a}{b}(\mathrm{e}^{-br}-1)+ar,$ (2)

where $a$ and $b$ are constants. The Toda lattice is a completely integrable
Hamiltonian system, allowing for the use of analytical tools to obtain and
study explicit solutions for the chain behaviour [33, 32, 41, 78, 24].

Granular chains describe the interaction of physical beads, and the
interaction potential is given by

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

where $\alpha>1$, $c$ is constant, and $\Delta$ is the equilibrium overlap of
adjacent particles due to precompression. Figure 1(a) illustrates a granular
chain with precompression. The interaction potential (3) is only nonzero
between particles that are in physical contact, and it is zero for
neighbouring particles that are not in contact.

The granular chain is important from a practical point of view, as it arises
in various applications. Previous theoretical, numerical, and experimental
studies have explored the generation [26, 49, 22, 23, 42], propagation [26,
49, 22, 23, 42], interaction [56, 44, 4, 6, 5, 25], and long-time dynamics [6,
5, 68, 65, 66] of solitary waves in granular chains.

A woodpile chain is composed of orthogonally stacked slender rigid cylinders,
illustrated in Figure 1(b). The elastic deformation in the direction
perpendicular to the stack direction is modelled by internal resonators, where
the mass and coupling constant of these resonators are determined by system
properties, such as the mass, shape, and material of the cylinders. For a
stack of identical cylinders, each resonator has the identical mass and
coupling constant. The stack of identical cylinders can be modelled as a
homogeneous chain with a prescribed interaction potential, where each chain
particle is also connected to an external particle by a linear spring. We
denote the mass of particles in the monatomic chain by $m_{1}$, the mass of
external particles by $m_{2}$, the mass ratio by $\eta^{2}=m_{2}/m_{1}$ and
the spring constant by $k$.

We consider woodpile chains with interaction potential along the direction of
the stack given by Toda (2) and power law (3) interactions. The configuration
of this idealized model of woodpile chains is shown in Figure 1(c). For the
remainder of this study, we use the term “woodpile chain” to describe this
idealized model.

(a) Precompressed particle chain (b) Physical woodpile (c) Idealized model

Figure 1: (a) Schematic illustration of a granular chain with precompression.
The chain consists of identical, aligned spheres with radius $R$, with a
compressive force $F_{0}$ applied at both ends of the chain to compress the
equilibrium particle positions. (b) The physical configuration of orthogonally
stacked rods known as a “woodpile chain”. (c) An idealized mathematical model
of the physical configuration in (b). In this model, the interaction between
heavy particles is given by (2) or (3), and light spherical particles
(sometimes called “resonators”) are attached to each heavy particle by a
spring.

The governing equations of an idealized woodpile chain are

$\displaystyle m_{1}\ddot{u}(n,t)=$
$\displaystyle\phi^{\prime}(x(n+1,t)-x(n,t))$
$\displaystyle-\phi^{\prime}(x(n,t)-x(n-1,t))-k[u(n,t)-v(n,t)]\,,$ (4)
$\displaystyle m_{2}\ddot{v}(n,t)=$ $\displaystyle k[u(n,t)-v(n,t)]\,,$ (5)

where $u(n,t)$ and $v(n,t)$, respectively, denote the displacement of the
$n^{\mathrm{th}}$ particle of mass $m_{1}$ and $m_{2}$ with respect to their
equilibrium positions at time $t$.

Previous studies [27, 28, 81, 46] have shown that traveling wave solutions in
woodpile chains are often not localized solitary waves, but are instead a type
of nonlinear wave known as a nanopteron. A nanopteron is the sum of a solitary
wave a train of exponentially small oscillations that extend indefinitely in
one or both directions of the central wave without decaying. These
oscillations typically appear across Stokes curves, and must be studied using
exponential asymptotic methods. In Figure 2, we illustrate examples of a
solitary wave and one- and two-sided nanoptera.

(a) A solitary wave (b) A one-sided nanopteron (c) A two-sided nanopteron

Figure 2: Comparison of the profiles of (a) a standard solitary wave, (b) a
one-sided nanopteron, and (c) a two-sided nanopteron that each propagates at
speed $c$. 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.

###  1.3 Rational approximation

Rational approximation methods approximate a function as a ratio of a
numerator and denominator polynomial. These methods are particularly accurate
for functions which contain singularities in or near the approximation domain.
Any approximation method that will be used for exponential asymptotics must be
able to accurately describe singular behaviour; rational approximation methods
are ideally suited for this purpose.

Historically, the most widely-used rational approximation method is Padé
approximation [7]. The approximation is written in terms of the power series
coefficients of the target function at a point. This method is well-suited to
problems where the power series coefficients can be obtained exactly, and
where we wish to approximate the solution accurately at a particular point.
For many problems we wish to study, we wish to approximate the solution over a
wider domain, and calculation of the power series coefficients requires the
numerical approximation of higher-order derivatives of the solution. Due to
these limitations, we apply the Padé approximation only to Toda woodpile
chains, for which the power series coefficients can be calculated exactly, for
purposes of comparison.

More recently-developed rational approximation methods, such as the Antoulas-
Anderson method [3], AAA method [60], vector fitting [39], rational Krylov
fitting (RKFIT) method [9, 10], and iterative rational Krylov algorithm (IRKA)
[8, 38], avoid some of the challenges faced when applying the Padé
approximation. These methods typically fit a rational approximation to a set
of data points distributed over a larger domain, with the pole locations
determined by some least-squares optimization process. The discussion in [60]
contains a discussion regarding the accuracy and speed of these methods. For
the purposes of this study, we will apply the AAA algorithm.

The Antoulas-Anderson method involves interpolating every pair of points to
determine a parametrization of all minimal-degree rational functions. The AAA
method [60] adapts this method to solve a least-squares minimization problem
to approximate the solution on a set of support points across the domain. The
AAA method has been applied to solve nonlinear eigenvalue problems [50],
represent conformal maps [76, 36], calculate rational minimax approximations
[31], and compute fractional diffusion [43]. In this paper, we will apply a
method based on AAA approximation to obtain leading-order approximations to
both Toda and woodpile nonlinear chains. The Padé and AAA approximation
methods are described in Section 3.

##  2 Exponential asymptotic method

The methodology used in this study is very similar to that of [55, 51, 27,
28]. The explanation of the methodology is similar to the explanation
contained in these previous studies.

We study the oscillatory tails of nanoptera in woodpile chains. We denote the
mass ratio between light particles and heavy particles, with mass $m_{2}$ and
$m_{1}$ respectively, by $\eta^{2}=m_{2}/m_{1}$. This system is singularly
perturbed in the limit $\eta\to 0$. The oscillatory tails are exponentially
small in $\eta$.

Exponential asymptotic methods were developed in [11, 12], and applied to
study Stokes’ phenomenon in special functions. The key observation was that if
an asymptotic series is truncated to minimise the error, the remainder term is
exponentially small. By rescaling the problem to study the remainder term, it
is possible to directly compute the exponentially small behaviour. These
techniques were further developed in [14] to further reduce the truncation
error. See [13] for a summary of these results.

The method applied in this study was developed in [17, 64]. We write the
solution $u$ to the singularly-perturbed governing equation as the power
series

$u\sim\sum_{j=0}^{\infty}\eta^{rj}u_{j}\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (6)

where $r$ is the number of times that $u_{j-1}$ must be differentiated to
obtain $u_{j}$. The leading-order solution, $u_{0}$, is obtained by setting
the small parameter $\eta$ to be equal to zero, and solving the resultant
equation.

We obtain an expression for $u_{j}$ by substituting the series (6) into the
governing equation and matching terms at powers of $\eta$. To determine
$u_{j}$ we need to differentiate earlier terms in the power series. If the
leading-order solution $u_{0}$ contains singular points, the repeated
differentiation guarantees that later terms in the series exhibit asymptotic
behaviour known as “factorial-over-power divergence” [29].

To capture the factorial-over-power divergence, [17] propose an asymptotic
ansatz for the terms $u_{j}$ in the limit $j\rightarrow\infty$, known as
“late-order terms”,

$u_{j}\sim\frac{U\Gamma(rj+\gamma)}{\chi^{rj+\gamma}}\quad\mathrm{as}\quad
j\to\infty\,.$ (7)

The parameter $\gamma$ is constant and $G$ and $\chi$ are functions of any
independent variables but are independent of $j$. The function $\chi$ is known
as the “singulant”. It is equal to $0$ at each singularity of the leading-
order solution $u_{0}$, ensuring that each $u_{j}$ is singular at the same
locations as the leading-order solution.

By substituting (6) and (7) into the governing equation and matching orders of
$\eta$, we can find $\chi$ and $G$. Through a local analysis of the solution
near singular points, we can obtain $\gamma$. Stokes curves originate at
singular points of $u_{0}$, which satisfy $\chi=0$, and follow curves on which
$\chi$ is real and positive.

The series (6) is then truncated, with the heuristic from [15] used to
determine the optimal truncation point $N_{\mathrm{opt}}$. This gives

$u=\sum_{j=0}^{N_{\mathrm{opt}}-1}\eta^{rj}u_{j}+u_{\exp}\,,$ (8)

where $u_{\mathrm{exp}}$ is the exponentially small truncation remainder. By
substituting (8) into the governing equation, we can obtain a new equation for
the exponentially small remainder term. Following [64], we can apply matched
asymptotic expansions in the neighbourhood of the Stokes curve to find that
the remainder has the form

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

where $\mathcal{S}$ is known as the Stokes multiplier; it is a function of the
independent variables in the problem. Away from the Stokes curve, the Stokes
multiplier takes constant value, and (9) reduces to the standard WKB ansatz.
In the neighborhood of the Stokes curve of width $\mathcal{O}(\sqrt{\eta})$
with $\eta\to 0$, the Stokes multiplier $\mathcal{S}$ undergoes a rapid
change, which is known as Stokes switching.

The leading-order behaviour $u_{0}$ plays several important roles in this
process. Firstly, the terms $u_{j}$ are singular at singular points of $u_{0}$
in the complex plane with the singularity strength increasing with $j$,
leading to factorial-over-power divergence. Secondly, Stokes curves in the
solution originate at singular points of $u_{0}$. Finally, from (9) we see
that the exponential scaling of $u_{\mathrm{exp}}$ is determined by the
appearance of $\chi$ in the exponent, and $\chi$ satisfies the condition that
$\chi=0$ at singular points of $u_{0}$.

The reliance of exponential asymptotic methods on the explicit calculation of
$u_{0}$ causes challenges when studying systems for which $u_{0}$ cannot be
calculated exactly. The purpose of this study is to determine whether
numerical approximations to the leading-order behaviour – in particular, the
numerical analytic continuation of a simulated solitary wave solution – can be
used to produce an approximation, denoted $\hat{u}_{0}$, that is able to
replicate exponential asymptotic results.

##  3 Rational approximation methods

In this section, we outline two methods for analytic continuation using
rational function approximation that we will use as part of a numerical
algorithm for finding the leading-order approximation $\hat{u}_{0}$.

###  3.1 Padé approximation

The $[N,M]$ Padé approximant, denoted $[N,M](x)$ to a function with power
series given by

$\displaystyle F(x)=\sum_{j=0}^{\infty}f_{j}x^{j}$ (10)

is written explicitly in terms of the power series coefficients $f_{j}$ in
(10) as $[N,M](x)=P(x)/Q(x)$, where $P(x)$ and $Q(x)$ are the determinants of
$N\times N$ matrix expressions given in [7]. The $[N,N](\xi)$ rational
approximation has the same power series as $F(x)$ up to order $N$. Increasing
the dimension of the matrix $N$, and therefore the accuracy of the
approximation, introduces more poles into the rational approximation.

For woodpile chains with Toda particle interactions, we can explicitly
determine the leading-order solution exactly, and it is therefore
straightforward to obtain $f_{j}$ exactly. For woodpile chains with power-law
interactions, we cannot find an exact solution. Instead, we have a numerical
leading-order solution. To obtain $f_{j}$, we need to calculate higher-order
terms in the power series numerically. This is impractical, as taking
differentiation magnifies accumulated numerical inaccuracy.

###  3.2 The AAA approximation

The Padé approximation determines the numerator and denominator polynomials by
evaluating the derivative of the function at a single point. We also consider
exponential asymptotics using the AAA method, which is based on fitting the
numerator and denominator polynomials at points sampled over a domain. We
summarise the first presentation of this method, from [60].

The AAA method approximates the function $f(z)$ as a rational expression,

$\displaystyle
f(z)\approx\frac{n(z)}{d(z)}={\sum_{j=1}^{m}\frac{w_{j}f_{j}}{z-z_{j}}}\,\Bigg{/}\,{\sum_{j=1}^{m}\frac{w_{j}}{z-z_{j}}}.$
(11)

The weights $w_{j}$ and the pole locations $z_{j}$ are found using an
iterative method. The $m^{\mathrm{th}}$ iteration takes as input the values of
$z_{j}$ for $j=1,\ldots,m$. The first part of the iteration process produces
the values of $w_{j}$ for $j=1,\ldots,m$, and the second part produces
$z_{m+1}$.

The first part of the process involves solving a linear least-squares problem
over a set of points in the restricted domain
$Z^{(m)}=Z\setminus\\{z_{1},\ldots,z_{m}\\}$; these points are labelled
$Z_{i}^{(m)}$. The use of the restricted domain ensures that the rational
approximation does not contain poles on the set. The function $f$ is
interpolated, giving $f_{1}=f(z_{1}),\dots,f_{m}=f(z_{m})$ at
$z_{1},\dots,z_{m}$. The weight vector $w=(w_{1},w_{2},\dots,w_{m})^{T}$ is
chosen to solve the least-squares problem of minimizing
$\mathinner{\\!\left\lVert fd-n\right\rVert}_{Z^{(m)}}$ subject to
$\mathinner{\\!\left\lVert w\right\rVert}_{m}=1$, where
$\mathinner{\\!\left\lVert\cdot\right\rVert}_{Z^{(m)}}$ is the discrete 2-norm
over $Z^{(m)}$ and $\mathinner{\\!\left\lVert\cdot\right\rVert}_{m}$. is the
discrete 2-norm on $m$-vectors. This is equivalent to solving the matrix
problem

$\displaystyle\mathrm{minimize}\,\,\,\,\mathinner{\\!\left\lVert
A^{(m)}w\right\rVert}_{M-m},\quad\mathinner{\\!\left\lVert
w\right\rVert}_{m}=1,$ (12)

where $A$ is a $(M-m)\times m$ matrix given by

$\displaystyle
A^{(m)}=\begin{bmatrix}\frac{F_{1}^{(m)}-f_{1}}{Z_{1}^{(m)}-z_{1}}&\dots&\frac{F_{1}^{(m)}-f_{m}}{Z_{1}^{(m)}-z_{m}}\\\
\vdots&\ddots&\vdots\\\
\frac{F_{M-m}^{(m)}-f_{1}}{Z_{M-m}^{(m)}-z_{1}}&\dots&\frac{F_{M-m}^{(m)}-f_{m}}{Z_{M-m}^{(m)}-z_{m}}\\\
\end{bmatrix}.$ (13)

This minimization problem can be solved efficiently using the singular value
decomposition (SVD), and taking $w$ as the final right singular vector in a
reduced SVD.

Finally, the value of $z_{m+1}$ is determined for the next iteration. It is
found by choosing the value of $z_{m+1}$ that maximizes the quantity

$\displaystyle\sum_{j=1}^{m}\frac{w_{j}f(z_{m+1})}{z_{m+1}-z_{j}}-\sum_{j=1}^{m}\frac{w_{j}f_{j}}{z_{m+1}-z_{j}}.$
(14)

This process is repeated until the $L^{\infty}$ norm of the difference between
the input function and the rational approximation is smaller than a specified
tolerance.

One potential obstacle in applying rational approximation methods is that
these methods produce meromorphic functions; all singularities in the
approximation are isolated poles. If the leading-order behaviour contains
branch points, the solution will not be accurately approximated in some local
vicinity of the pole. In previous studies of AAA approximation [37, 76, 77],
it was noted that the presence of branch points in the solution creates an
exponential clustering of poles (and zeroes) in the approximation around the
branch points of the target function, which accurately approximate the effects
of a branch cut. A goal of this study is to determine whether this pole
accumulation is capable of accurately predicting exponentially small behaviour
in Toda woodpile chains, whose leading-order solution contains logarithmic
branch cuts.

##  4 Woodpile chains with Toda interactions

The scaled governing equations for a woodpile chain governed by a Toda
interaction potential (2), with $a=1$ and $b=1$, are given by

$\displaystyle\ddot{u}(n,t)$
$\displaystyle=\mathrm{e}^{-[u(n,t)-u(n-1,t)]}-\mathrm{e}^{-[u(n+1,t)-u(n,t)]}-k[{u}(n,t)-{v}(n,t)]\,,$
(15) $\displaystyle\eta^{2}\ddot{{v}}(n,t)$
$\displaystyle={k}[{u}(n,t)-{v}(n,t)]\,.$ (16)

In the limit $\eta\rightarrow 0$ we expand $u(x,t)$ and $v(x,t)$ as a power
series in $\eta^{2}$, giving

$\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)\,.$ (17)

Substituting the series expression (17) into (16) and matching at the leading
order in the limit $\eta\to 0$, we find that $u_{0}(n,t)=v_{0}(n,t)$, and

$\displaystyle\ddot{u}_{0}(n,t)$
$\displaystyle=\mathrm{e}^{-[u_{0}(n,t)-u_{0}(n-1,t)]}-\mathrm{e}^{-[u_{0}(n+1,t)-u_{0}(n,t)]}\,.$
(18)

We introduce a co-moving frame $\xi=n-ct$, where $c$ is the solitary wave
velocity. In the co-moving frame, (18) can be written as

$\displaystyle c^{2}{u}^{\prime\prime}_{0}(\xi)$
$\displaystyle=\mathrm{e}^{-[u_{0}(\xi)-u_{0}(\xi-1)]}-\mathrm{e}^{-[u_{0}(\xi+1)-u_{0}(\xi)]}\,,$
(19)

where prime denotes derivative with respect to $\xi$.

Equation (19) has an exact solitary-wave solution, which we present in (21).
This exact solution will be applied to determine the form of the exponentially
small oscillations in the nanopteron solutions for the woodpile system using
exponential asymptotics. We will then use different numerical approaches to
approximate $u_{0}$, giving an approximate leading-order solution (denoted
$\hat{u}_{0}$), and compare the output of exponential asymptotic analyses
based on these approximate solutions.

The leading-order solution is singular at points in the complex plane, and
these singularities occur in complex-conjugate pairs. We denote the
singularity locations as $\xi_{s}$ and $\xi_{s}^{*}$, where
$\mathrm{Im}(\xi_{s})>0$. Each singularity pair is connected by a Stokes
curve, which switches on an exponentially small contribution.

Assume that $u_{0}$ is singular at $\xi=\xi_{s}$, and that the singularity
behaviour is locally either logarithmic, so that
$u_{0}\sim\mu\log(\xi-\xi_{s})$ as $\xi\to\xi_{s}$ for some $\mu$, or a pole
or branch point, so that $u_{0}\sim\mu(\xi-\xi_{s})^{\nu}$ as $\xi\to\xi_{s}$
for some $\mu$ and $\nu$. This local behaviour is sufficient to determine the
exponentially small contribution that appears across the Stokes curve; this
analysis is presented in Appendix A. The exponentially small solution
contribution to $v(\xi)$, denoted $v_{\mathrm{exp}}$, is given by

$v_{\mathrm{exp}}\sim\frac{2|V|\pi}{\eta^{\nu}}\exp\left(-\frac{\sqrt{k}\mathrm{Im}(\xi_{s})}{c\eta}\right)\cos\left(\frac{\sqrt{k}(\xi-\mathrm{Re}(\xi_{s}))}{c\eta}-\phi\right)\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (20)

where $V=\mu$ if the singularity is logarithmic, and
$V=\mathrm{i}\mu\sqrt{k}/c$ if the solution is a pole or branch point. In the
case of a logarithmic singularity, we take $\nu=0$. The quantity $\phi$ is
given by $\pi/2$ if the singularity is logarithmic, and the phase of $V$ if
the singularity is a pole or branch cut. The term $u_{\mathrm{exp}}$ can be
obtained directly from $v_{\mathrm{exp}}$, but we omit these details here.

This contribution is present for $\mathrm{Re}(\xi)>\mathrm{Re}(\xi_{s})$,
which corresponds to the region behind the leading-order solitary wave. The
amplitude of the oscillations decays exponentially as $\mathrm{Im}(\xi_{s})$
increases. The dominant oscillations are therefore those associated with the
singularity pairs that are nearest to the real axis, which are exponentially
larger in the small $\eta$ limit than singularity pairs which are further
away.

###  4.1 Exact leading-order solution

(a) Leading-order solution $u_{0}(\xi)$ for $\kappa=1$ (b)
$\mathrm{Re}(u_{0})$ (c) $\mathrm{Im}(u_{0})$ (d) Stokes Curves

(e) Comparison of Results

Figure 3: (a) The exact leading-order wave $u_{0}$ as a function $\xi$ for
$\kappa=1$. (b)–(c) The behaviour of $u_{0}$ for $\xi\in\mathbb{C}$ in the
upper half-plane. (d) The Stokes curves that connect the singularity pairs. As
Stokes curves are crossed from right to left, exponentially small oscillations
appear in the solution. In ① there are no oscillations. In ②, there is one
oscillatory contribution. In ③, both oscillatory contributions are present.
When these oscillations have opposite phase, they cancel. Region ③ contains no
oscillations in this case, and the solution is a solitary wave. (e) Comparison
between asymptotic predictions and simulations of the oscillation amplitudes
for $\kappa=1$ and $\kappa=2$ respectively, with $k=1$. The asymptotic
predictions accurately predict the oscillation amplitude and the values of
$\eta$ at which the oscillations disappear.

Equation (18) is a scaled version of the integrable Toda chain and may be
solved exactly, as in [74]. Recalling that $u_{0}=v_{0}$, the leading-order
solution to the system (15)–(16) is given by

$\displaystyle u_{0}(n,t)=\log\left(\frac{\cosh(\kappa(n-1)-\beta
t)}{\cosh(\kappa n-\beta t)}\right)\,,\quad u_{0}(n,t)=v_{0}(n,t)\,,$ (21)

where $\kappa$ is the soliton parameter, $\beta=\sinh\kappa$ and the velocity
of the soliton is given by $c=\sinh\kappa/\kappa$. We write this in terms of
an offset co-moving frame, given by $\xi=n-ct-1/2$, where the offset is
introduced to make the wave symmetric about $\xi=0$. The solution in this
frame is given by

$\displaystyle
u_{0}(\xi)=\log\left(\frac{\cosh(\kappa(\xi-1/2))}{\cosh(\kappa(\xi+1/2))}\right)\,,\quad
u_{0}(\xi)=v_{0}(\xi)\,.$ (22)

The leading-order solution $u_{0}$ is shown in Figure 3(a).

If we analytically continue the leading-order solutions $u_{0}$ and $v_{0}$ so
that $\xi\in\mathbb{C}$, we see that these functions are singular at a set of
points described by

$\displaystyle\xi_{1,\pm,j}=-\tfrac{1}{2}\pm\tfrac{\mathrm{i}\pi}{2\kappa}(2j+1)\quad\mathrm{and}\quad\xi_{2,\pm,j}=\tfrac{1}{2}\pm\tfrac{\mathrm{i}\pi}{2\kappa}(2j+1),$
(23)

where $j$ is a non-negative integer. The behaviour of $u_{0}$ for
$\xi\in\mathbb{C}$ in the upper half-plane is shown in Figure 3(b)–(c). These
contour plots show the real and imaginary part of $u_{0}(\xi)$, and depict the
two singularities nearest to the real axis. Recall that the exponentially
small behaviour is dominated exponentially as $\eta\to 0$ by contributions
from the singularities that are nearest to the real axis.

These singularities are connected by Stokes curves, shown in Figure 3(d).
There are two exponentially small contributions; one appears on the left-hand
side of the curve $\mathrm{Re}(\xi)=-1/2$, and one appears on the left-hand
side of the curve $\mathrm{Re}(\xi)=1/2$. The oscillatory behaviour in the
nanopteron solution is therefore a superposition of the two sets of
exponentially small oscillations.

From (20), we determine that the exponentially small oscillations generated by
the singularity pair located at $\xi=\xi_{2,\pm,j}$ have the form

$v_{\mathrm{exp}}\sim 2\pi\exp\left(-\frac{\sqrt{k}\pi(2j+1)}{2\kappa
c\eta}\right)\sin\left(\frac{\sqrt{k}(\xi-1/2)}{c\eta}\right)\quad\mathrm{as}\quad\eta\rightarrow
0\,.$ (24)

The exponentially small oscillations generated by singularity pair located at
$\xi=\xi_{1,\pm,j}$ have the form from (24), with an offset in $\xi$ of $+1/2$
rather than $-1/2$. From these two expressions, it is clear that the dominant
oscillations are those obtained by choosing $j=0$. Taking the sum of these
oscillations and simplifying them allows us to determine the amplitude of the
combined oscillatory behaviour, given by

$\mathrm{Amplitude}\sim 4\pi\exp\left(-\frac{\sqrt{k}\pi(2j+1)}{2\kappa
c\eta}\right)\sin\left(\frac{\sqrt{k}}{2c\eta}\right).$ (25)

A comparison between the asymptotic prediction of the amplitude and numerical
results calculated using the numerical scheme from [27, 28] is shown in Figure
3(e). The asymptotic predictions show strong agreement with the numerical
simulations, especially for small values of $\eta$. The asymptotic analysis is
able to accurately predict values of the mass ratio at which the oscillations
disappear entirely (which appear in Figure 3(e) as sharp negative cusps);
these correspond to solitary-wave solutions.

###  4.2 Tanh-fitted leading-order solution

(a) Leading-order approximation $\hat{u}_{0}(\xi)$ for $\kappa=1$ (b)
$\mathrm{Re}(\hat{u}_{0})$ (c) $\mathrm{Im}(\hat{u}_{0})$ (d) Leading-order
approximation $\hat{u}_{0}(\xi)$ for $\kappa=2$

(e) $\mathrm{Re}(\hat{u}_{0})$ (f) $\mathrm{Im}(\hat{u}_{0})$

Figure 4: A comparison between the tanh-fitted leading order solution
$\hat{u}_{0}$ and the exact leading order solution $u_{0}$ for (a) $\kappa=1$,
and (d) $\kappa=2$. The analytic continuation of $\hat{u}_{0}$ is shown in
(b)–(c) for $\kappa=1$ and (e)–(f) for $\kappa=2$. For $\kappa=1$,
$\hat{u}_{0}$ contains one dominant singularity pair. For $\kappa=2$,
$\hat{u}_{0}$ contains two equally dominant singularity pairs. Hence,
$\hat{u}_{0}$ cannot produce solutions for which the oscillations cancel if
$\kappa=1$, but this is possible for $\kappa=2$.

In [69], the authors note that solitary waves often have a kink shape, and
proposed a method of approximating solitary waves by fitting a tanh function
to the solution. The authors assume solitary waves in Hertz chains can be
approximated as $u_{0}\approx\hat{u}_{0}$, where $\hat{u}_{0}$ has the form

$\displaystyle\hat{u}_{0}(\xi)=-\frac{A}{2}\tanh(f(\xi)),\quad
f(\xi)=\sum_{n=0}^{N}C_{2n+1}\xi^{2n+1},\quad\xi=n-c_{A}t,$ (26)

where $A$ is the amplitude of the solitary wave, $c_{A}$ is the velocity of
the solitary wave with amplitude $A$, $N$ is a positive integer, $\xi$ defines
the co-moving frame and the coefficients $C_{2n+1}$ are determined by
numerically approximating higher-order derivatives of the leading-order
behaviour and equating the results with an analytic expression for the
corresponding derivatives of (26). This approximation is appropriate a
solitary wave profile for which $u_{0}(\xi)\to\mp A/2$ as $\xi\to\pm\infty$.
It is clear from Figure 3(a) that the Toda solitary wave falls into this
category.

For $\kappa=1$, we obtain $C_{1}=0.92423$, $C_{3}=0.020876$ and a negative
value for $C_{5}$. For $\kappa=2$, we obtain $C_{1}=1.5232$, $C_{3}=0.32505$,
and $C_{5}<0$. If negative coefficients are included in the approximation, the
profile behaves correctly in the vicinity of $\xi=0$, but the wave profile is
not monotonic and can tend to incorrect values as $|\xi|\to\infty$. There is a
limit to the approximation accuracy that can be obtained using this method,
and it cannot be continued indefinitely. Hence, we only consider
approximations using coefficients up to $C_{3}$.

We show a comparison the approximation $\hat{u}_{0}$ and the true leading-
order solution $u_{0}$ in Figure 4(a) and (d) for $\kappa=1$ and $\kappa=2$.
The behaviour of $\hat{u}_{0}$ in the complex plane is presented in Figure
4(b)–(c) for $\kappa=1$ and (e)–(f) for $\kappa=2$. These figures show a
difference between the two cases. For $\kappa=1$ there is a single dominant
singularity pair, located at $\xi=\pm 1.84035\mathrm{i}$. For $\kappa=2$,
there are two equally dominant singularity pairs, located at $\xi=-0.4923\pm
1.2817\mathrm{i}$ and $\xi=0.4923\pm 1.2817\mathrm{i}$.

In fact, there exists a critical value $\kappa_{c}\approx 1.56$ for $\kappa$.
If $\kappa<\kappa_{c}$, then $\hat{u}_{0}$ contains one singularity pair
closest to the real axis; when $\kappa>\kappa_{c}$, the single singularity
pair splits into two distinct pairs. Oscillations can only cancel if there are
two contributions of the same amplitude. Hence, the oscillations can only
vanish if $\kappa$ exceeds this critical value. From the results in Figure 25,
we know this is not a true feature of the asymptotic approximation, but rather
an inaccuracy introduced by the approximation method.

Using (20) and the local behaviour of $\hat{u}_{0}$ near singularities, the
amplitude of the oscillations for $\kappa=1$ is given by

$\displaystyle\mathrm{Amplitude}\sim\frac{1.4043\pi\sqrt{k}}{\eta
c_{A}}\exp\left(-\frac{1.84035\sqrt{k}}{c_{A}\eta}\right).$ (27)

For $\kappa=2$, the amplitude of the oscillations is given by

$\displaystyle\mathrm{Amplitude}\sim\frac{1.6120\pi\sqrt{k}}{\eta
c_{A}}\exp\left(-\frac{1.2817\sqrt{k}}{c_{A}\eta}\right)\cos\left(\frac{0.4923\sqrt{k}}{c_{A}\eta}-1.4435\right).$
(28)

We show a comparison between the amplitude based on the approximation
$\hat{u}_{0}$ and that based on the exact leading-order solution $u_{0}$ in
Figure 5. For $\kappa=1$, which is below the critical value, the asymptotic
prediction cannot identify values of $\eta$ at the oscillations vanish, and
the approximation underestimates the amplitude. For $\kappa=2$, the asymptotic
prediction identifies the values of $\eta$ at which the oscillations vanish,
but the amplitude deviates significantly from numerical simulations.

Figure 5: Comparison of asymptotic predictions based on the approximate
leading-order $\hat{u}_{0}$ and predictions based on $u_{0}$. For $\kappa=1$,
there are no values of $\eta$ which result the oscillations vanishing, while
for $\kappa=2$, the predicted values of $\eta$ at which the oscillations
vanish are accurate. In both cases, the amplitude is not accurately predicted
by the exponential asymptotic method.

###  4.3 Padé-approximated leading-order solution

Figure 6: Padé approximants with different values of $N$, which is the order
of the numerator and denominator polynomials. The size of the region in which
the approximant agrees well with the exact solution increases as $N$
increases. For $N\geq 8$, the approximant is accurate in the region containing
the solitary wave, although the approximants generated using $N=8$, $16$ do
differ visibly from the exact solution for $|\xi|>5$.

In this section we use the Padé approximant described in (3.1) to approximate
the leading-order solution. As we are approximating a function which tends to
constant values at $\pm\infty$, we expect that the degree of the polynomial in
the numerator will be close to that in the denominator, so we take $N=M$.

In Figure 6 we show that the size of the region in which the Padé approximant
accurately predicts the solution behaviour increases with $N$. We wish to
select an approximation that accurately describes the behaviour of the
solution in the region where the wave varies. From Figure 6, we see that
choosing $N=8$ is accurately describes the solitary wave in $\xi\in[-5,5]$. We
could choose higher values of $N$, but the subsequent analysis does not change
significantly.

For $\kappa=1$, the Padé approximant is used to generate $\hat{u}_{0}$, given
by

$\displaystyle\hat{u}_{0}(\xi)=[8,8](\xi)=-\frac{0.9242\xi+0.5441\xi^{3}+0.1132\xi^{5}+0.005663\xi^{7}}{1+0.8509\xi^{2}+0.2742\xi^{4}+0.03295\xi^{6}+0.0004223\xi^{8}}.$
(29)

This approximation is shown in Figure 7(a), and the analytic continuation is
shown in Figure 7(b)–(c). The solution contains two singularity pairs closest
to the real axis, located at $\xi=-0.4559\pm 1.6651\mathrm{i}$ and
$\xi=0.4559\pm 1.6651\mathrm{i}$. The solution also contains subdominant
pairs, one of which can be seen in Figure 7(b)–(c). These pairs will generate
Stokes curves, shown in Figure 7(d). The oscillations will be dominated by
contributions from the dominant pair, but the remaining poles will also switch
on behaviour that does not necessarily cancel. The amplitude will nearly
cancel at particular values of $\eta$, but the cancellation will not be
perfect.

(a) Exact leading-order solution $u_{0}(\xi)$ for $\kappa=1$ (b)
$\mathrm{Re}(\hat{u}_{0})$ (c) $\mathrm{Im}(\hat{u}_{0})$ (d) Stokes Curves

(e) Comparison of Results

Figure 7: (a) Approximate leading-order wave $\hat{u}_{0}$ for $\kappa=1$.
(b)–(c) The behaviour of $\hat{u}_{0}$ for $\xi\in\mathbb{C}$ in the upper
half-plane. (d) The Stokes curves that connect the singularity pairs. In ①
there are no oscillations. In ②, there is one oscillatory contribution. In ③,
a second oscillatory contribution appears which has different amplitude to the
first. In ④, a third oscillatory contribution is switched on. When the first
and third oscillations are in opposite phase, they cancel. There are still
oscillations present due to the second Stokes curve, but they are smaller than
the contributions from the dominant poles. (e) Comparison between asymptotic
predictions for $u_{0}$ and $\hat{u}_{0}$ for $\kappa=1$ and $\kappa=2$. The
asymptotic predictions accurately predict both the amplitude and the values of
$\eta$ at which the dominant oscillations vanish.

Using (20), the amplitude of the oscillations generated by $\hat{u}_{0}$ is

$\displaystyle\mathrm{Amplitude}\sim\frac{1.311\pi\sqrt{k}}{\eta
c_{A}}\exp\left(-\frac{1.6651\sqrt{k}}{c_{A}\eta}\right)\cos\left(\frac{0.4559\sqrt{k}}{c_{A}\eta}+2.0766\right).$
(30)

This expression does not contain the exponentially subdominant oscillations
generated by the subdominant pole pairs. Hence, we don’t expect the amplitude
to truly disappear at values of $\eta$ for which this expression is equal to
zero, as the subdominant contributions will not cancel at these points.

For $\kappa=2$, the Padé approximant is used to generate $\hat{u}_{0}$, given
by

$\displaystyle\hat{u}_{0}(\xi)=[8,8](\xi)=-\frac{3.04638\xi+259.567\xi^{3}+182.411\xi^{5}+160.653\xi^{7}}{85.7651\xi^{2}+107.787\xi^{4}+102.879\xi^{6}+25.8326\xi^{8}}.$
(31)

The singularities closest to the real axis are located at $\xi=-0.5205\pm
0.08846\mathrm{i}$ and $\xi=0.5205\pm 0.8846\mathrm{i}$. The amplitude of the
oscillations is given by (20) as

$\displaystyle\mathrm{Amplitude}\sim\frac{1.1690\pi\sqrt{k}}{\eta
c_{A}}\exp\left(-\frac{0.8846\sqrt{k}}{c_{A}\eta}\right)\cos\left(\frac{0.5205\sqrt{k}}{c_{A}\eta}+1.4161\right).$
(32)

In Figure 7(e), the oscillation amplitude generated by $\hat{u}_{0}$ for
$\kappa=1$ and $\kappa=2$ are compared with the oscillation amplitude
generated by the exact solution $u_{0}$. These amplitudes were generated using
all singularity pairs in the solution, rather than just the dominant pairs. In
each case, the approximation accurately captures the asymptotic behaviour
generated by the exact leading-order. It is able to correctly predict both the
amplitude of the oscillations, and the values of $\eta$ at which the
oscillations appear to cancel. This is not true cancellation, as there are
still exponentially subdominant oscillations in the solution; however, these
contributions are so small that they are not visibly apparent in the figure.

These results indicate that using rational approximation to perform numerical
analytic continuation produces approximations that are useful for exponential
asymptotics. Computing the Padé approximation requires finding Taylor
coefficients of the leading-order solution. For the woodpile chain with Toda
interactions, these coefficients are straightforward to obtain. For a woodpile
chain with power-law interactions, these coefficients can only be obtained
numerically by approximating higher-order derivatives at the point $\xi=0$
numerically.

###  4.4 AAA-approximated leading-order solution

Motivated by the accurate results produced by Padé approximation, we formulate
a method based on rational approximation for $u_{0}$ using the AAA algorithm;
an implementation of this algorithm is included in the Chebfun Matlab package
[30]. The first step is to compute the leading-order solitary wave $u_{0}$
numerically by exciting a monatomic Toda lattice with a velocity impulse. The
second step is to use the AAA algorithm to obtain a rational approximation for
the leading-order wave, $\hat{u}_{0}$.

####  4.4.1 Numerical Leading-Order Solution

For the numerical step, we use the velocity Verlet algorithm [79, 1] to
predict the particle behaviour. This is a symplectic integrator that is
designed to conserve the energy of a system. 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\,,$ (33)

where the time step is given by $\Delta t=t_{n}-t_{n-1}$, 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, the acceleration
$\boldsymbol{a}_{n}$ can be obtained directly from $\boldsymbol{x}_{n}$. Each
simulation contains a large number of particles on a periodic domain. The
domain information and the size of the time steps is presented in Table 1.

Table 1: Parameters for simulations of leading-order solitary wave in the Toda woodpile chain. $\kappa$ | Number of particles | $\Delta t$ | Velocity impulse
---|---|---|---
1 | 300 | 0.01 | $\boldsymbol{v}_{61}(0)=2.3525$
2 | 300 | 0.01 | $\boldsymbol{v}_{61}(0)=7.2560$

Initially, all particles in the chain are at rest. At the initial time step, a
velocity impulse is applied to the particle at $n=61$. The velocity impulse
generates a solitary wave and dispersive waves. The dispersive waves travel
slower than the solitary wave, eventually seperating entirely. This is
illustrated in Figure 8, which shows the motion of the particles at $n=160$,
$200$, and $240$. As the distance from the initial impulse grows, the
seperation between the solitary wave and the dispersive waves increases. At
$n=161$, the seperation is sufficient that we can isolate the numerical
approximation to the solitary wave without capturing any effects due to the
dispersive waves. We tune the magnitude of the velocity impulse to obtain
solitary waves with amplitude $4$ and $2$, corresponding to $\kappa=2$ and
$\kappa=1$. The parameters are given in Table 1.

Figure 8:  Displacement of particles $n=$ 160, 200, and 240 as a function of
time. The time $t$ is offset by a shift ($\delta$) so that the solitary wave
profiles are overlaid. The solitary wave profile in simulations is shifted up
by half the amplitude, as the particles start with zero displacement. The
solitary wave travels faster than the dispersive waves, and the two wave
contributions seperate. Eventually, the seperation is large enough that we can
isolate the solitary wave, without any measurable effects due to dispersive
waves.

(a) AAA-approximated $\hat{u}_{0}(\xi)$ for $\kappa=1$ (b)
$\mathrm{Re}(\hat{u}_{0})$ (c) $\mathrm{Im}(\hat{u}_{0})$ (d) AAA-approximated
$\hat{u}_{0}(\xi)$ for $\kappa=2$

(e) $\mathrm{Re}(\hat{u}_{0})$ (f) $\mathrm{Im}(\hat{u}_{0})$

Figure 9: (a) A comparison between the AAA approximated function $\hat{u}_{0}$
and the exact leading-order solution $u_{0}$ for $\kappa=1$. The analytic
continuation of $\hat{u}_{0}$ is shown in (b)–(c). The analytically-continued
approximation contains two dominant singularity pairs, with three additional
singularities located between them. Figures (d)–(f) show the same quantities
for $\kappa=2$. Again, the analytically-continued approximation contains two
dominant singularity pairs, with three additional singularities located
between them.

####  4.4.2 Exponential Asymptotic Analysis

We now apply the AAA algorithm to the numerical leading-order solitary wave to
obtain an explicit rational approximation $\hat{u}_{0}$. To do this, we
consider the solution behaviour within a finite domain, which we denote as
$[-d,d]$, and use this as the input to the AAA algorithm. The spatial
resolution, which we denote as $r$, is given by $\Delta t\cdot\xi$. For
$\kappa=1$, we select $d=4.7004$ and $r=0.0118$. For $\kappa=2$, we select
$d=3.6232$ and $r=0.0181$. The error tolerance of the AAA algorithm represents
the $L^{\infty}$ norm of the difference between the input function and the
rational approximation – reducing the size of the error tolerance has the
effect of increasing the order of the polynomials in the numerator and
denominator of the rational expression, and therefore the number of
singularity pairs in the approximation. For the purpose of this analysis, we
select an error tolerance of $10^{-13}$.

Figure 9(a) and (d) show a comparison between the approximated rational
function $\hat{u}_{0}$ and the exact leading-order solution $u_{0}$ for
$\kappa=1$ and $\kappa=2$ respectively. The analytic continuation of
$\hat{u}_{0}$ is shown in Figure 9(b)–(c) and (e)–(f) for $\kappa=1$ and
$\kappa=2$ respectively. In both cases, there two dominant singularity pairs,
with three more singularity pairs spread along a curve between the dominant
pairs. This suggests that leading-order solitary wave profiles with a similar
shape produce similar analytically-continued behaviour when approximated using
the AAA algorithm. This was not true in the tanh-fitted approximation (4), in
which the analytic continuation changed significantly between $\kappa=1$ and
$\kappa=2$.

We use the approximated function $\hat{u}_{0}$ as the leading-order solution
in an exponential asymptotic analysis. Each singularity pair generates a set
of oscillations with the form given by (20). The total oscillatory
contribution behind the leading-order wave is given by the superposition of
the oscillations generated by each pair. The results of this analysis are
presented in in Figure 10 for $\kappa=1$ and $2$. These predicted behaviour is
visually indistinguishable from the results obtained using $u_{0}$, and the
AAA-approximation method accurately predicts both the amplitude of the
oscillations and the values of $\eta$ at which the oscillations vanish.

Figure 10: Comparison between the asymptotic amplitude predicted using the
AAA-approximation for the leading-order, $\hat{u}_{0}$, and the exact leading-
order behaviour, $u_{0}$, for $\kappa=1$ and $\kappa=2$. The predicted
amplitudes show very strong agreement.

To demonstrate the effect of adjusting the tolerance, we repeated this method
with tolerances of $10^{-5}$ and $10^{-9}$. The results of this are shown in
Figure 11. If the tolerance is relaxed, the approximation becomes less
accurate. To quantify the accuracy of each simulation, we introduce an error
expression

$\displaystyle
R=\frac{\|\log_{10}(\mathrm{Amplitude}_{\mathrm{exact}})-\log_{10}(\mathrm{Amplitude}_{\mathrm{AAA}})\|}{\|\log_{10}(\mathrm{Amplitude}_{\mathrm{AAA}})\|},$
(34)

where $\|\textbf{f}\|^{2}=\sum_{i=1}^{N}f_{i}^{2}$ and
$\textbf{f}=(f_{1},f_{2},\dots f_{N})$. The amplitude vectors contain the
oscillation amplitude at 90001 different values of $\eta$, evenly spaced in
the domain $[0.03,0.12]$. This expression was chosen as a quantitative measure
of the visual fit of the logarithmic data. For AAA error tolerances of
$10^{-5}$, $10^{-9}$ and $10^{-13}$, the computed values of $R$ are $0.0695$,
$0.0232$ and $0.0104$ respectively. This indicates that the fit increases in
accuracy as the tolerance is reduced, as would be expected.

(a) Tolerance = $10^{-5}$ (b) Tolerance = $10^{-9}$ (c) Tolerance = $10^{-13}$

Figure 11: Comparison of oscillation amplitudes found using the AAA-
approximated $\hat{u}_{0}$ and the exact leading-order $u_{0}$ with error
tolerances in the AAA algorithm of (a) $10^{-5}$, (b) $10^{-9}$ and (c)
$10^{-13}$. The relative error $R$ is $0.0695$, $0.0232$ and $0.0104$ for (a),
(b) and (c) respectively.

To illustrate the effects of adjusting the sample domain and resolution, we
repeat this algorithm with different values for the sample resolution $r$ and
the domain $d$, and calculate the error $R$ in each case. A summary of these
calculations is presented graphically in Figure 12. Parameter choices which
produce an error of $R<0.025$ correspond to output which is visibly similar to
that given in Figure 11(b)–(c), while $R>0.05$ typically produces output which
is similar to that of Figure 11(a). Figure 12 shows that the error can be
reduced by increasing the width of the sample domain or decreasing the sample
resolution.

Figure 12: The relative error $R$ for the amplitude of the trailing
oscillations for different values of the sample domain $[-d,d]$ and the
resolution $r$. Increasing the width of the sample domain or decreasing the
sample resolution improves the accuracy of the asymptotic behaviour predicted
by the exact leading order solution.

In Figure 13, we compare the location of singularities in the analytic
continuation of (a) the exact solution, (b) the tanh-fitted approximation, (c)
the Padé approximant, and (d) the AAA approximation for the leading order
behaviour. We see from this image that the AAA approximation most accurately
predicts the location of the dominant singularity pairs. The location of the
dominant singularity pair strongly influences on both the amplitude and phase
of the oscillations, as seen in the asymptotic expression (20). This is
consistent with the observation that the behaviour of the oscillations that
was generated using a AAA approximation for the leading order behaviour agrees
very strongly with the asymptotic behaviour produced using the exact
solutions; furthermore, it outperforms the other tested methods.

Singularities in the exact solutions are logarithmic, while the singularities
in each of the approximations for the leading-order behaviour are simple
poles. This does not appear to have a significant impact on the woodpile
problem, due to the linear coupling between the two governing equations. This
will not necessarily be true for more complicated nonlinear problems – we will
discuss this in Section 6.

(a) Exact: $\mathrm{Re}(u_{0})$ (b) Tanh: $\mathrm{Re}(\hat{u}_{0})$ (c) Padé:
$\mathrm{Re}(\hat{u}_{0})$ (d) AAA: $\mathrm{Re}(\hat{u}_{0})$

Figure 13: The locations of singularities in the analytic continuation of the
(a) exact solution $u_{0}$, and (b)–(d) the three different approximated
solutions obtained in this section, with $\kappa=1$ and singular points
indicated by crosses. The amplitude of the oscillations is determined by
$\mathrm{Im}(\xi_{s})$, and the phase is determined by $\mathrm{Re}(\xi_{s})$.
In the exact leading-order behaviour, the dominant singularity pairs are
located in at $\xi=-1/2\pm\mathrm{i}\pi/2$ and $\xi=1/2\pm\mathrm{i}\pi/2$.
AAA approximation most accurately predicts the dominant singularity locations,
and is therefore able to most accurately reproduce the asymptotic predictions
made using the exact solution.

##  5 Woodpile chain with power-law interaction potential

In the previous section, we established that the exponentially small
oscillations present in nanoptera can be calculated using a AAA approximation
in place of the exact leading-order solitary wave. We now apply this method in
order to predict behaviour in a system where no exact leading-order solution
exists.

The governing equations for woodpile chains with power-law interaction
potential (3), with $c=1/(\alpha+1)$ for simplicity, are given by

$\displaystyle\ddot{u}(n,t)=$
$\displaystyle[\Delta+u(n-1,t)-u(n,t)]_{+}^{\alpha}$
$\displaystyle-[\Delta+u(n,t)-u(n+1,t)]_{+}^{\alpha}-k[{u}(n,t)-{v}(n,t)]\,,$
(35) $\displaystyle\eta^{2}\ddot{{v}}(n,t)=$
$\displaystyle{k}[{u}(n,t)-{v}(n,t)]\,.$ (36)

For this section, we consider woodpile chains with zero precompression, as in
[28]. We therefore set $\Delta=0$.

In the limit $\eta\rightarrow 0$ we expand $u(x,t)$ and $v(x,t)$ as asymptotic
power series in $\eta^{2}$, as in (17). We substutitue the series expression
(17) into (36) match at the leading order in the limit $\eta\to 0$, and
introduce the co-moving frame $\xi=n-ct$, where $c$ is the solitary wave
velocity. We find that $u_{0}(\xi)=v_{0}(\xi)$, and

$\displaystyle c^{2}{u}^{\prime\prime}_{0}(\xi)$
$\displaystyle=[u(\xi-1)-u(\xi,t)]_{+}^{\alpha}-[u(\xi)-u(\xi+1)]_{+}^{\alpha}\,,$
(37)

where prime denotes derivative with respect to $\xi$. In contrast to woodpile
chains with Toda interaction potential, $u_{0}(\xi)$ cannot be expressed
analytically.

Assuming we have some leading-order solution $u_{0}$, and we analytically
continue the expression so that $\xi\in\mathbb{C}$, we will again find
singularity pairs in the complex plane that generate Stokes curves. We label
the singularity locations as $\xi=\xi_{s}$ and $\xi_{s}^{*}$, where
$\mathrm{Im}(\xi_{s})>0$. We assume that the local behaviour of the solution
near the singularity is given by $u_{0}\sim\mu(\xi-\xi_{s})^{\nu}$ as
$\xi\to\xi_{s}$, for some $\mu$ and $\nu$. A similar exponential asymptotic
analysis to Section 4 shows that the exponentially small contribution has an
identical form to the oscillations in the Toda woodpile chain, presented in
(20).

In [28], the authors approximate the leading-order solution using tanh-
fitting, motivated by the work of [69]. Applying an exponential asymptotic
analysis allows the amplitude to be accurately predicted for a number of
different parameters; however, the asymptotic predictions differ significantly
from the numerical simulation for smaller interaction exponents, such as the
Hertzian interaction ($\alpha=1.5$). In this section, we will show that the
behaviour of the system can be accurately predicted if we instead use a AAA
approximation for the leading-order behaviour.

We obtain the leading-order solitary wave by applying a velocity impulse to a
monatomic particle chain as in Section 4.4. Unlike the Toda woodpile chain
simulations, we fix the size of the velocity impulse so that
$\boldsymbol{v}_{61}(0)=5$. We then determine the amplitude of the solitary
wave as part of the solution. The simulations were performed with identical
parameters to the Toda chain, from Table 1.

We then use the AAA algorithm to obtain a rational approximation for the
leading-order behaviour, which we denote $\hat{u}_{0}$. For the purpose of
this study, we use the parameters $d=4.6296$ and $r=0.0231$, and the AAA
tolerance is $10^{-13}$. In Figure 14(a), we compare $\hat{u}_{0}$ to the
numerical solution.

(a) Leading-order $\hat{u}_{0}(\xi)$ for $\alpha=1.5$ (b)
$\mathrm{Re}(\hat{u}_{0})$ (c) $\mathrm{Im}(\hat{u}_{0})$

Figure 14: (a) A comparison between the AAA approximated function
$\hat{u}_{0}$ and the exact leading-order solution $u_{0}$ for $\kappa=1$. The
analytic continuation of $\hat{u}_{0}$ is shown in (b)–(c).

We analytically continue $\hat{u}_{0}$, and determine the location of the
singularites in the complex plane. These singularities are all simple poles,
and are shown in Figure 15. We also illustrate the size of the residual
associated with each singularity. In each case, there are pairs which are
further from the real axis than the dominant singularity pair, but which have
much larger residuals. This means that they can still have a measurable effect
on the oscillatory behaviour.

This complicated pole structure explains why the tanh-fitting method was
unable to accurately capture details of the oscillation behaviour. The tanh-
fitting method did not accurately predict singularity locations that are
further from the real axis than the dominant singularity pair, and these
singularities have a significant affect on the solution behaviour due to their
large residuals.

(a) Singularities: $\alpha=1.5$ (b) Singularities: $\alpha=1.75$ (c)
Singularities: $\alpha=2$ (d) Singularities: $\alpha=2.5$

(e) Singularities: $\alpha=3$ (f) Singularities: $\alpha=3.5$

Figure 15: Singularities with positive imaginary part for
$\alpha=1.5,1.75,2,2.5,3,3.5$.

(a) $\alpha=1.5$ (b) $\alpha=1.75$ (c) $\alpha=2$ (d) $\alpha=2.5$

(e) $\alpha=3.5$ (f) $\alpha=3.5$

Figure 16: Comparison of the oscillation amplitudes predicted by asymptotics
based on AAA-approximated leading-order solutions (black), based on leading
order with tanh-fitted leading-order solutions (red, dashed) and based on
numerical simulations (blue dots). For larger values of $\alpha$, both the AAA
and tanh-fitted asymptotic predictions are accurate. For smaller values of
$\alpha$, predictions made based on AAA-approximation are more accurate than
tanh-fitting. The AAA-approximation method can accurately describe small
variations in the amplitude due to interactions between subdominant
contributions. This is particularly apparent at $\eta\approx 0.075$ for
$\alpha=1.5$.

To determine the behaviour of the exponentially small oscillations, we take
the sum of the contributions from every singularity pair, and then determine
the amplitude of this combined expressions. The contributions for each
singularity pair are given by (20). In Figure 16, we present the results for a
range of interaction exponents $\alpha$, and compare these predictions to
those made in [28] using the tanh-fitted method. For larger values of $\alpha$
both methods accurately predict the system behaviour, but for smaller values
of $\alpha$, the predictions made using AAA approximation are significantly
more accurate than those made using tanh-fitting.

Using AAA approximation to generate $\hat{u}_{0}$ also allows the asymptotic
behaviour that occurs due to subdominant poles. This can be seen most clearly
in Figure 16(a), which shows the results for $\alpha=1.5$. In addition to
accurately predicting the values of $\eta$ at which cancellation occurs, our
results also correctly describe the smaller variations in the amplitude that
can be seen between these values. These smaller variations occur due to the
cancellation of subdominant oscillatory contributions, generated by
singularity pairs that are further from the real axis than the dominant pairs.

##  6 Conclusion

We applied three approximation methods to study nanoptera in Toda woodpile
chains – tanh-fitting, Padé approximants, and the AAA approximation algorithm
– and compared the results with asymptotic predictions made using the exact
leading-order solution. We found that, while all three methods were capable of
describing aspects of the exponentially small behaviour, the AAA approximation
was most effective in determining two key features of the exponentially small
oscillations in a Toda woodpile chains: the amplitude of the oscillations, and
the values of the mass ratio at which the oscillations vanish.

Using the AAA algorithm to approximate the leading-order behaviour has several
advantages over tanh-fitting and Padé approximants. Approximations made using
tanh-fitting cannot be made arbitrarily precise, as we cannot take series
coefficients indefinitely. Padé approximants rely entirely on the behaviour of
the numerical solution at a point, which requires taking progressively higher
numerical derivatives to increase the accuracy of the rational approximation.
AAA approximation avoids this by fitting the solution over a specified domain,
rather than a single point. AAA approximation made the most accurate
predictions of the singularity locations, explaining its accurate prediction
of exponentially small behaviour.

We used the AAA algorithm to approximate the leading-order behaviour of
woodpile chains with power-law interaction. The subsequent asymptotic
predictions agree with numerics well for all values of the interaction
exponent $\alpha$. This is notable, as [28] studied this problem using tanh-
fitting methods, and found that the method did not accurately describe the
oscillations for moderate or small values of the interaction exponent,
including Hertzian interactions ($\alpha=1.5$).

For $\alpha=1.5$ our asymptotic prediction also captures details in the
asymptotic behaviour that were missed by the tanh-fitting method. By looking
at the pole locations, it is apparent why this occured. This method was able
to identify smaller variations in the amplitude that are caused by
interactions between subdominant poles that are not accurately identified
using tanh-fitting – in some cases, such as $\alpha=1.5$, these poles have
large residues, and therefore can visibly affect the solution behaviour.
Asymptotic behaviour generated by this method can include significant effects,
even if they are exponentially subdominant compared to the oscillation
amplitude.

The next step to establishing the viability of this method is to apply the AAA
approximation to determine the leading-order behaviour for more complicated
model systems such as diatomic particle chains [55, 51], which includes
additional challenges.The AAA approximation is effective at locating
singularities, which it represents as simple poles. In the case of branch
points, the approximation contains a cluster of simple poles that accurately
approximate the branch cut behaviour. The order of the singularity plays a
significant role in nonlinear problems, and it is likely that the local
approximation of $u_{0}$ near the singularity will need to be determined. Data
from [77] indicates that the manner in which the poles cluster near the
singular point depends on the order of the singularity. If the order of the
singularity can be determined from this clustering, the appropriate branch cut
could be included in the approximation function with the weight determined
through the least-squares process.

Evaluating the singulant for nonlinearly coupled chains will require the
evaluation of complex contour integrals of $u_{0}$. Figure 17 shows the
difference between the AAA approximation for the Toda chain and the exact
leading-order behaviour, which is small except in a region near the poles. The
singularity positions must be considered in prescribing an integral contour,
to control the accuracy of the approximation.

(a) $\log_{10}|\mathrm{Re}(u_{0}-\hat{u}_{0})|$ (b)
$\log_{10}|\mathrm{Im}(u_{0}-\hat{u}_{0})|$

Figure 17: Error of the AAA approximant $\hat{u}_{0}$ of the leading-order
behaviour $u_{0}$ for the Toda lattice with $\kappa=1$. The poles of the exact
solution are shown as black crosses. The AAA approximant closely approximates
the exact solution except in a region near the poles.

##  Appendix A Detailed exponential asymptotic analysis

The method described in this section is very similar to that of [55, 51, 27,
28], and we therefore omit much of the intermediate analysis.

###  A.1 Late-order terms

In the limit $\eta\rightarrow 0$ we write $u(x,t)$ and $v(x,t)$ as asymptotic
power series in $\eta^{2}$, giving the series expressions from (17). Writing
the governing equations (15)–(16) in terms of $\xi$ and matching at each order
of $\eta$ gives

$\displaystyle
c^{2}u_{j}^{\prime\prime}(\xi)=(u_{j}(\xi-1)-u_{j}(\xi))\exp{(u_{0}(\xi-1)-u_{0}(\xi))}$
$\displaystyle\quad-(u_{j}(\xi)-u_{j}(\xi+1))\exp{(u_{0}(\xi)-u_{0}(\xi+1))}-k[u_{j}(\xi)-v_{j}(\xi)]+\ldots\,,$
(38) $\displaystyle
c^{2}v_{j-1}^{\prime\prime}(\xi)=k[u_{j}(\xi)-v_{j}(\xi)]\,.$ (39)

We only retain terms containing $u_{j}$, $v_{j}$, and derivatives of $u_{j}$
and $v_{j-1}$, as these terms govern the late-order asymptotic behaviour.

We assume that $u_{0}$ contains a singularity at a point in the complex plane
$\xi=\xi_{s}$. We consider two cases: the case where the singularity in
$u_{0}$ is logarithmic, so that $u_{0}\sim\mu\log(\xi-\xi_{s})$ as
$\xi\to\xi_{s}$, and the case the singularity has asymptotic behaviour
$u_{0}\sim\mu(\xi-\xi_{s})^{\nu}$ as $\xi\to\xi_{s}$. In the case of the exact
solution (21), $\mu=1$ for the singularities at $\xi=\xi_{1,\pm,j}$ and
$\mu=-1$ for the singularities at $\xi=\xi_{2,\pm,j}$.

We pose a late-order ansatz that consists of sums of terms with the form

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

where $r$ is the number of times that $u_{j-1}$ and $v_{j-1}$ must be
differentiated to obtain $u_{j}$ and $v_{j}$, as illustrated in Section 2.
From (A.1)–(39) we have $r=2$. We set $\chi=0$ at $\xi=\xi_{s}$ so that late-
order terms are singular at the same location as the leading-order solution.
By substituting the late-order ansatz (40) into (A.1), we find that only
$\beta_{1}+2=\beta_{2}$ can produce a nontrivial asymptotic balance. This
implies that $u_{j}=\mathcal{O}(v_{j-1})$ as $j\rightarrow\infty$, and
therefore that $v_{j}\gg u_{j}$ as $j\to\infty$.

Applying the late-order ansatz (40) to (39) gives

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

where the omitted terms are $\mathcal{O}(v_{j-1})$ in the $j\rightarrow\infty$
limit. Matching terms at $\mathcal{O}(v_{j})$ in the $j\rightarrow\infty$
limit, we obtain $c^{2}(\chi^{\prime}(\xi))^{2}=-k$, implying
$\chi^{\prime}(\xi)=\pm\mathrm{i}\sqrt{k}/c$, and therefore

$\chi(\xi)=\pm\tfrac{\mathrm{i}}{c}\sqrt{k}(\xi-\xi_{s})\,.$ (42)

Stokes curves occur where $\mathrm{Im}(\chi)=0$ and $\mathrm{Re}(\chi)>0$
[14]. This corresponds to the positive sign choice for singularities in the
upper half plane, and the negative sign choice for those in the lower half
plane. We do not consider late-order terms associated with the remaining sign
choices, as they do not generate Stokes curves.

Matching terms at $\mathcal{O}(v^{\prime}_{j-1})$, we obtain the prefactor
equation $2V^{\prime}(\xi)\chi^{\prime}(\xi)=0$; $V$ is therefore constant.
The strength of the singularities must be consistent with the leading-order
behaviour in the limit that $\xi\to\xi_{s}$. For these expressions to be
consistent, we require $\beta_{2}=\nu$ and $\beta_{1}=\nu+2$.

We follow a procedure similar to that of [27, 28], and determine $V$ by
matching the outer expansion for the late-order terms (40) with an inner
solution in the neighbourhood of $\xi_{s}$. Through this analysis, we find
that $V=\mu$ if the singularity is logarithmic, and
$V=\mathrm{i}\mu\sqrt{k}/c$ if the singularity is a pole or a branch point.

###  A.2 Stokes switching

Truncating the series (17) after $N$ terms gives

$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)\,,$ (43)

where $S_{N}$ and $R_{N}$ are the remainder terms obtained by truncating the
series, and are exponentially small if we optimally truncate the series. We
denote the optimal truncation point as $N=N_{\mathrm{opt}}$. The heuristic in
[15] shows that $N_{\mathrm{opt}}$ is determined by finding the point where
consecutive terms in the series are equal in size, giving
$N_{\mathrm{opt}}=|\chi|/2\eta+\omega$, where $\omega\in[0,1)$ is chosen to
ensure that $N_{\mathrm{opt}}$ is an integer.

Inserting (43) into (15)–(16) gives governing equations for $R_{N}$ and
$S_{N}$ as $\eta\to 0$:

$\displaystyle
c^{2}S^{\prime\prime}_{N}(\xi)\sim\,[S_{N}(\xi+1)-S(\xi)]\mathrm{e}^{-[u_{0}(\xi+1)-u_{0}(\xi)]}$
$\displaystyle\quad\quad-[S_{N}(\xi)-S_{N}(\xi-1)]\mathrm{e}^{-[u_{0}(\xi)-u_{0}(\xi-1)]}-k[S_{N}(\xi)-R_{N}(\xi)]+\ldots,$
(44)
$\displaystyle\eta^{2}c^{2}R^{\prime\prime}(\xi)\sim\,k[S_{N}(\xi)-R_{N}(\xi)]-\eta^{2N}v^{\prime\prime}_{N-1}(\xi),$
(45)

where omitted terms are asymptotically subdominant as $\eta\to 0$. We also
find from a balancing argument that $R_{N}(\xi)\gg S_{N}(\xi)$ in this limit,
and we may therefore simplify this expression – in particular, the second
equation decouples from the first. Noting that $N_{\mathrm{opt}}$ is large as
$\eta\to 0$, we apply (40) for the late-order terms, giving

$\eta^{2}c^{2}R^{\prime\prime}(\xi)+k\chi^{\prime}(\xi)^{2}R_{N}(\xi)\sim-\eta^{2N}\frac{V\Gamma(2N+\beta_{2})}{\chi(\xi)^{2N+\beta_{2}}}\qquad\mathrm{as}\qquad\eta\to
0.$ (46)

The right-hand side is asymptotically subdominant compared to the balance on
the left-hand side except for a narrow region in the neighbourhood of the
Stokes curve. We use a WKB analysis to show that the remainder can be written
as $R_{N}\propto\mathrm{e}^{-\chi/\eta}$ as $\eta\to 0$ away from the Stokes
curve. To capture the behaviour of the solution near the Stokes curve, we
write the ansatz $R_{N}=A(\xi)\mathrm{e}^{-\chi/\eta}$ as $\eta\to 0$. By
writing $\chi=r\mathrm{e}^{\mathrm{i}\theta}$, the variation of $A$ in terms
of $\theta$ is given by

$\frac{\mathrm{d}A}{\mathrm{d}\theta}\sim\frac{V\chi^{\prime}\Gamma(2N+\beta_{2})}{2\chi^{2N+\beta_{2}}}\mathrm{e}^{\chi/\eta}\qquad\mathrm{as}\qquad\eta\to
0.$ (47)

The Stokes curve follows $\theta=0$. We can define an inner region near the
Stokes curve in terms of an inner variable $\hat{\theta}=\eta^{-1/2}\theta$,
and use asymptotic matching to determine the variation of $A$ in this region.
We find that $A=0$ on the left-hand side of the Stokes curve, and
$A\sim\mathrm{i}\pi V/\eta^{\nu}$ as $\xi\to 0$ on the right-hand side of the
Stokes curve.

The contribution given by the singularity at $\xi=\xi_{s}^{*}$ is the complex
conjugate of the contribution given by the singularity at $\xi=\xi_{s}$.
Hence, the exponentially small behaviour that appears on the left-hand side of
the Stokes curve connecting the singularity pair $\xi_{s}$ and $\xi_{s}^{*}$
is

$R_{N}\sim\frac{2|V|\pi}{\eta^{\nu}}\exp\left(-\frac{\sqrt{k}\,\mathrm{Im}(\xi_{s})}{c\eta}\right)\cos\left(\frac{\sqrt{k}(\xi-\mathrm{Re}(\xi_{s}))}{c\eta}-\phi\right)\quad\mathrm{as}\quad\eta\rightarrow
0\,,$ (48)

where $\nu=0$ if the singularity is logarithmic. The quantity $\phi$ is given
by $\pi/2$ if the singularity is logarithmic, and the phase of $V$ if it is a
pole or branch cut.

## Acknowledgments

The authors acknowledge the support of Australian Research Council Discovery
Project #190101190.

## References

  * [1] M. P. Allen and D. J. Tildesley.  Computer Simulation of Liquids.  Clarendon Press, Oxford, UK, 1987.
  * [2] N. Andersson and C. J. Howls.  The asymptotic quasinormal mode spectrum of non-rotating black holes.  Classical and Quant. Grav., 21:1623–1642, 2003.
  * [3] A. C. Antoulas and B. D. Q. Anderson.  On the Scalar Rational Interpolation Problem.  IMA J. Math. Control I., 3(2-3):61–88, 1986.
  * [4] E. Ávalos and S. Sen.  How solitary waves collide in discrete granular alignments.  Phys. Rev. E, 79:046607, 2009.
  * [5] E. Ávalos and S. Sen.  Granular chain between asymmetric boundaries and the quasiequilibrium state.  Phys. Rev. E, 89:053202, 2014.
  * [6] 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.
  * [7] G. A. Baker.  The Padé approximant method and some related generalizations.  In G. A. Baker and J. L. Gammel, editors, The Padé approximant in theoretical physics, chapter 1, pages 1–39. Academic Press, New York and London, 1970.
  * [8] C. Beattie and S. Gugercin.  Realization-independent $\mathcal{H}_{2}$-approximation.  In 2012 IEEE 51st IEEE Conference on Decision and Control (CDC), pages 4953–4958, 2012.
  * [9] M. Berljafa and S. Güttel.  Generalized rational Krylov decompositions with an application to rational approximation.  SIAM J. Matrix Anal. Appl., 36(2):894–916, 2015.
  * [10] M. Berljafa and S. Güttel.  The RKFIT algorithm for nonlinear rational approximation.  SIAM J. Sci. Comput., 39(5):A2049–A2071, 2017.
  * [11] M. V. Berry.  Stokes’ phenomenon; smoothing a Victorian discontinuity.  Pub. Math. de L’IHÉS, 68:211–221, 1988.
  * [12] M. V. Berry.  Uniform asymptotic smoothing of Stokes’s discontinuties.  Proc. Roy. Soc. Lond. A, 422(1862):7–21, 1989.
  * [13] 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.
  * [14] M. V. Berry and C. J. Howls.  Hyperasymptotics.  Proc. Roy. Soc. Lond. A, 430(1880):653–668, 1990.
  * [15] J. P. Boyd.  The devil’s invention: Asymptotic, superasymptotic and hyperasymptotic series.  Acta Appl. Math., 56(1):1–98, 1999.
  * [16] R. C. Brower, D. A. Kessler, J. Koplik, and H. Levine.  Geometrical approach to moving-interface dynamics.  Phys. Rev. Lett., 51:1111–1114, 1983.
  * [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] S. J. Chapman, P. H. Trinh, and T. P. Witelski.  Exponential asymptotics for thin film rupture.  SIAM J. Appl. Math., 73:232–253, 2013.
  * [19] S. J. Chapman and J. M. Vanden-Broeck.  Exponential asymptotics and capillary waves.  SIAM J. Appl. Math., 62(6):1872––1898, 2002.
  * [20] S. J. Chapman and J. M. Vanden-Broeck.  Exponential asymptotics and gravity waves.  J. Fluid Mech., 567:299–326, 2006.
  * [21] S.J. Chapman and G. Kozyreff.  Exponential asymptotics of localised patterns and snaking bifurcation diagrams.  Physica D, 238(3):319–354, 2009.
  * [22] C. Coste, E. Falcon, and S. Fauve.  Solitary waves in a chain of beads under Hertz contact.  Phys. Rev. E, 56:6104–6117, 1997.
  * [23] 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.
  * [24] P. Deift, S. Kamvissis, T. Kriecherbauer, and X. Zhou.  The Toda rarefaction problem.  Commun Pur. Appl. Math., 49(1):35–83, 1996.
  * [25] G. Deng, G. Biondini, and S. Sen.  Interactions of solitary waves in integrable and nonintegrable lattices.  Chaos, 30(4):043101, 2020.
  * [26] 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.
  * [27] G. Deng, C. J. Lustri, and M. A. Porter.  Nanoptera in weakly nonlinear woodpile chains and diatomic granular chains.  SIAM J. Appl. Dyn. Syst., 20(4):2412–2449, 2021.
  * [28] G. Deng and C.J. Lustri.  Nanoptera in nonlinear woodpile chains with zero precompression.  Physica D, 429:133053, 2022.
  * [29] R. B. Dingle.  Asymptotic Expansions: Their Derivation and Interpretation.  Academic Press, New York, NY, USA, 1973.
  * [30] T. A. Driscoll, N. Hale, and L. N. Trefethen, editors.  Chebfun Guide.  Pafnuty Publications, Oxford, 2014.
  * [31] S. I. Filip, Y. Nakatsukasa, L. N. Trefethen, and B. Beckermann.  Rational minimax approximation via adaptive barycentric representations.  SIAM J. Sci. Comput., 40(4):A2427–A2455, 2018.
  * [32] H. Flaschka.  On the Toda Lattice. II: Inverse-Scattering Solution.  Prog. Theor. Phys., 51(3):703–716, 1974.
  * [33] H. Flaschka.  The Toda lattice. II. Existence of integrals.  Phys. Rev. B, 9:1924–1925, 1974.
  * [34] G. Friesecke and J. A. D. Wattis.  Existence theorem for solitary waves on lattices.  Commun. Math. Phys., 161(2):391–418, 1994.
  * [35] S. Giller.  Topological expansion and exponential asymptotics in 1d quantum mechanics.  J. Phys. A-Math. Gen., 33(8):1543–1580, 2000.
  * [36] A. Gopal and L. N. Trefethen.  Representation of conformal maps by rational functions.  Numer. Math., 142:359–382, 2019.
  * [37] A. Gopal and L. N. Trefethen.  Solving laplace problems with corner singularities via rational functions.  SIAM J. on Numer. Anal., 57(5):2074–2094, 2019.
  * [38] S. Gugercin, A. C. Antoulas, and C. A. Beattie.  $\mathcal{H}_{2}$ model reduction for large-scale linear dynamical systems.  SIAM J. Matrix Anal. Appl., 30:609–638, 2008.
  * [39] B. Gustavsen and A. Semlyen.  Rational approximation of frequency domain responses by vector fitting.  IEEE Transactions on power delivery, 14(3):1052–1061, 1999.
  * [40] G. A. Hagedorn and A. Joye.  Semiclassical dynamics with exponentially small error estimates.  Commun. Math. Phys., 207:439–465, 1999.
  * [41] M. Hénon.  Integrals of the Toda lattice.  Phys. Rev. B, 9:1921–1923, 1974.
  * [42] E. J. Hinch and S. Saint-Jean.  The fragmentation of a line of balls by an impact.  P. R. Soc. A, 455:3201–3220, 1999.
  * [43] C. Hofreither.  A unified view of some numerical methods for fractional diffusion.  Comput. Math. Appl., 80(2):332–350, 2020.  Numerical Methods for Scientific Computations and Advanced Applications II.
  * [44] 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.
  * [45] A. Joye.  Exponential asymptotics in a singular limit for $n$-level scattering systems.  SIAM J. Math. Anal., 28(3):669–703, 1997.
  * [46] 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.
  * [47] Martin D. Kruskal and Harvey Segur.  Asymptotics beyond all orders in a model of crystal growth.  Stud. Appl. Math., 85(2):129–181, 1991.
  * [48] F. Lawless.  Exponential asymptotics and spectral theory for optical tunnelling.  PhD thesis, Dublin City University, 1993.
  * [49] 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.
  * [50] P. Lietaert, K. Meerbergen, J. Pérez, and B. Vandereycken.  Automatic rational approximation and linearization of nonlinear eigenvalue problems.  IMA J. Numer. Anal., 2021.
  * [51] C. J. Lustri.  Nanoptera and Stokes curves in the 2-periodic Fermi–Pasta–Ulam–Tsingou equation.  Physica D, 402:132239, 2020.
  * [52] C. J. Lustri and S. J. Chapman.  Steady gravity waves due to a submerged source.  J. Fluid Mech., 732:660–686, 2013.
  * [53] C. J. Lustri. and S. J. Chapman.  Unsteady flow over a submerged source with low Froude number.  Eur. J. Appl. Math., 25(5):655–680, 2014.
  * [54] C. J. Lustri, C. C. Green, and S. W. McCue.  Selection of a Hele-Shaw bubble via exponential asymptotics.  SIAM J. Appl. Math., 80:289–311, 2020.
  * [55] C. J. Lustri and M. A. Porter.  Nanoptera in a period-2 Toda chain.  SIAM J. Appl. Dyn. Syst., 17(2):1182–1212, 2018.
  * [56] F. S. Manciu and S. Sen.  Secondary solitary wave formation in systems with generalized Hertz interactions.  Phys. Rev. E, 66:016616, 2002.
  * [57] V. P. Maslov.  Global exponential asymptotics of solutions of the tunnel equations and problems of large deviations.  Matematicheskii Institut imeni Steklova Trudy, 163:150–180, 1984\.
  * [58] R. E. Meyer.  On exponential asymptotics for nonseparable wave equations I: Complex geometrical optics and connection.  SIAM J. Appl. Math., 51(6):1585–1601, 1991.
  * [59] R. E. Meyer.  On exponential asymptotics for nonseparable wave equations II: EBK quantization.  SIAM J. Appl. Math., 51(6):1602–1615, 1991.
  * [60] Y. Nakatsukasa, O. Sète, and L. N. Trefethen.  The AAA algorithm for rational approximation.  SIAM J. Sci. Comput., 40:1494–1522, 2018.
  * [61] V. F. Nesterenko.  Propagation of nonlinear compression pulses in granular media.  J. Appl. Mech. Tech. Phy., 24:733–743, 1983.
  * [62] V. F. Nesterenko.  Dynamics of Heterogeneous Materials.  Springer-Verlag, Heidelberg, Germany, 2001.
  * [63] S. Nixon and J. Yang.  Exponential asymptotics for solitons in PT-symmetric periodic potentials.  Stud. Appl. Math., 133:373–397, 2014.
  * [64] 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.
  * [65] M. Przedborski, T. A. Harroun, and S. Sen.  Granular chains with soft boundaries: Slowing the transition to quasiequilibrium.  Phys. Rev. E, 91:042207, 2015.
  * [66] M. Przedborski, S. Sen, and T. A. Harroun.  The equilibrium phase in heterogeneous Hertzian chains.  J. Stat. Mech., 2017(12):123204, 2017.
  * [67] S. Sen, J. Hong, J. Bang, E. Ávalos, and R. Doney.  Solitary waves in the granular chain.  Phys. Rep., 462(2):21–66, 2008.
  * [68] 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.  Proceedings of the VIII Latin American Workshop on Nonlinear Phenomena.
  * [69] 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.
  * [70] A. Stefanov and P. G. Kevrekidis.  On the existence of solitary traveling waves for generalized Hertzian chains.  J. Nonlinear Sci., 22:327–349, 2012.
  * [71] M. Toda.  Vibration of a chain with nonlinear interaction.  J. Phys. Soc. Jpn., 22(2):431–436, 1967.
  * [72] M. Toda.  Wave propagation in anharmonic lattices.  J. Phys. Soc. Jpn., 23(3):501–506, 1967.
  * [73] M. Toda.  Waves in nonlinear lattice.  Progress of Theoretical Physics Supplement Prog. Theor. Phys. Supp., 45:174–200, 1970.
  * [74] M. Toda.  Theory of Nonlinear Lattices.  Springer-Verlag, New York/Berlin, 1981.
  * [75] M. Toda and M. Wadati.  A soliton and two solitons in an exponential lattice and related equations.  J. Phys. Soc. Jpn., 34(1):18–25, 1973.
  * [76] L. N. Trefethen.  Numerical conformal mapping with rational functions.  Comput. Meth. Funct. Th., 20:369–387, 2020.
  * [77] L. N. Trefethen, Y. Nakatsukasa, and J. A. Weideman.  Exponential node clustering at singularities for rational approximation, quadrature, and PDEs.  Numer. Math., pages 1–28, 2021.
  * [78] S. Venakides, P. Deift, and R. Oba.  The Toda shock problem.  Commun. Pur. Appl. Math., 44(8-9):1171–1242, 1991.
  * [79] L. Verlet.  Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules.  Phys. Rev., 159:98–103, 1967.
  * [80] A. D. Wood.  Exponential asymptotics and spectral theory for curved optical waveguides.  In Asymptotics Beyond All Orders, pages 317–326. Springer, 1991\.
  * [81] 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.