Nanoptera and Stokes Curves in the 2-Periodic Fermi-Pasta-Ulam-Tsingou Equation

TL;DR

This paper investigates nonlocal solitary waves called nanoptera in a period-2 FPUT lattice, using exponential asymptotics to analyze wave trains and identify conditions for localized solutions, validated by numerical simulations.

Contribution

It introduces an exponential asymptotic method to analyze nanoptera and Stokes curves in the 2-periodic FPUT equation, revealing conditions for localized solitary waves.

Findings

Identification of Stokes curves where oscillations appear

Conditions for mass ratios that eliminate oscillations

Validation of asymptotic predictions with numerical simulations

Abstract

This work presents asymptotic solutions to a singularly-perturbed, period-2 FPUT lattice and uses exponential asymptotics to examine `nanoptera', which are nonlocal solitary waves with constant-amplitude, exponentially small wave trains which appear behind the wave front. Using an exponential asymptotic approach, this work isolates the exponentially small oscillations, and demonstrates that they appear as special curves in the analytically-continued solution, known as `Stokes curves' are crossed. By isolating these the asymptotic form of these oscillations, it is shown that there are special mass ratios which cause the oscillations to vanish, producing localized solitary-wave solutions. The asymptotic predictions are validated through comparison with numerical simulations.

Tables

Table 1. Table 1: Comparison between δ 𝛿 \delta values for nanopteron-free solutions computed with ϵ = 0.25 italic-ϵ 0.25 \epsilon=0.25 , and the asymptotic prediction obtained from ( 58 ) for corresponding values of K 𝐾 K . The asymptotic prediction provides a good approximation for the numerically-obtained solutions. As in the amplitude predictions, this is not as accurate as the predictions for the diatomic Toda lattice from [ 37 ] , even for larger values of K 𝐾 K (and hence smaller values of δ 𝛿 \delta ).

$K$	Computed $\delta$	Asymptotic $\delta$
2	0.3922	0.4258
3	0.2809	0.3041
4	0.2248	0.2365
5	0.1932	0.1935
6	0.1677	0.1638

Equations

$$m(j)\ddot{x}(j,t)=F(x(j+1,t)-x(j,t))-F(x(j,t)-x(j-1,t)),$$

$$F(r)=r+r^2.$$

$$g_r\sim \frac{G\Gamma (r+\gamma )}{{\chi }^{r+\gamma }}\mathrm{as}r\to \infty ,$$

$$g_{\exp }\sim \mathcal{S}G{\mathrm{e}}^{-\chi /\delta },$$

$$r(j,t)=x(j+1,t)-x(j,t),p(j,t)=\dot{x}(j,t),$$

$$\dot{r}(j,t)=p(j+1,t)-p(j,t),m_j\dot{p}(j,t)=F(r(j,t))-F(r(j-1,t)),$$

$$(r(j,t),p(j,t))=3{ϵ}^2{\mathrm{sech}}^2(\beta ϵ(j-c_{ϵ}t))(1,-c_0)+\mathcal{O}({ϵ}^{5/2}),$$

$$\beta =\sqrt{\frac{3(m_1^2+2m_1{m}_2+m_2^2)}{2(m_1^2-m_1{m}_2+m_2^2)}},c_{ϵ}=(1+{ϵ}^2)\sqrt{\frac{2}{m_1+m_2}}.$$

$$x(j,t)=\frac{3ϵ}{\beta }\tanh (\beta ϵ(j-c_{ϵ}t))+\mathcal{O}({ϵ}^{3/2}).$$

$$c_{ϵ}^2{y}^{\prime \prime }(\xi )$$

$${\delta }^2{c}_{ϵ}^2{z}^{\prime \prime }(\xi )$$

$$c_{ϵ}^2{y}^{\prime \prime }(\xi )$$

$$z(\xi )$$

$$y(\xi )=\sum _{r=0}^{\infty }{\delta }^{2r}{y}_r(\xi ),z(\xi )=\sum _{r=0}^{\infty }{\delta }^{2r}{z}_r(\xi ),$$

$$c_{ϵ}^2{y}_0^{\prime \prime }(\xi )$$

$$0$$

$$y_0(\xi )=\sqrt{6}ϵ\tanh \left(\sqrt{\frac{3}{2}}ϵ\xi \right)+\mathcal{O}({ϵ}^{3/2}),z_0(\xi )=\frac{1}{2}(y(\xi +1)+y(\xi -1)),$$

$${\xi }_{s,N,\pm }=\frac{(2N-1)\pi \mathrm{i}}{\sqrt{6}ϵ}\pm 1,N\in \mathbb{Z},$$

$$z_0(\xi )\sim (\xi -{\xi }_{s,N,\pm }{)}^{-1}\mathrm{as}\xi \to {\xi }_{s,N,\pm },$$

$$c_{ϵ}^2{y}_r^{\prime \prime }(\xi )=$$

$$+(z_1(\xi +1)-y_1(\xi ))(z_{r-1}(\xi +1)-y_{r-1}(\xi ))F^{\prime \prime }(z_0(\xi +1)-y_0(\xi ))$$

$$-(y_1(\xi )-z_{r-1}(\xi -1))(y_1(\xi )-z_{r-1}(\xi -1))F^{\prime \prime }(y_0(\xi )-z_0(\xi -1))+\dots ,$$

$$c_{ϵ}^2{z}_{r-1}^{\prime \prime }(\xi )=$$

$$+(y_1(\xi +1)-z_1(\xi ))(y_{r-1}(\xi +1)-z_{r-1}(\xi ))F^{\prime \prime }(y_0(\xi +1)-z_0(\xi ))$$

$$-(z_1(\xi )-y_1(\xi -1))(z_{r-1}(\xi )-y_{r-1}(\xi -1))F^{\prime \prime }(z_0(\xi )-y_0(\xi -1))+\dots ,$$

$$y_r(\xi )\sim \frac{Y(\xi )\Gamma (2r+\alpha )}{\chi (\xi {)}^{2r+\alpha }},z_r(\xi )\sim \frac{Z(\xi )\Gamma (2r+\beta )}{\chi (\xi {)}^{2r+\beta }},\mathrm{as}r\to \infty ,$$

$$\frac{c_{ϵ}^2({\chi }^{\prime }(\xi ){)}^{2}Z(\xi )\Gamma (2r+\alpha )}{\chi (\xi {)}^{2r+\alpha }}-\frac{2c_{ϵ}^2{\chi }^{\prime }(\xi )Z^{\prime }(\xi )\Gamma (2r+\alpha -1)}{\chi (\xi {)}^{2r+\alpha -1}}-\frac{c_{ϵ}^2{\chi }^{\prime \prime }(\xi )Z(\xi )\Gamma (2r+\alpha -1)}{\chi (\xi {)}^{2r+\alpha -1}}$$

$$=$$

$$+\frac{Y(\xi +1)\Gamma (2r+\alpha -2)}{\chi (\xi +1{)}^{2r+\alpha -2}}(y_1(\xi +1)-z_1(\xi ))F^{\prime }(\frac{1}{2}(y_0(\xi +1)-y_0(\xi -1)))$$

$$+\frac{Y(\xi -1)\Gamma (2r+\alpha -2)}{\chi (\xi -1{)}^{2r+\alpha -2}}(z_1(\xi )-y_1(\xi -1))F^{\prime }(\frac{1}{2}(y_0(\xi +1)-y_0(\xi -1)))+\dots ,$$

$$c_{ϵ}^2({\chi }^{\prime }{)}^2=-2(1+y_0(\xi +1)-y_0(\xi -1)).$$

$$\chi =\pm \frac{\mathrm{i}\sqrt{2}}{c_{ϵ}}{\int }_{{\xi }_s}^{\xi }\sqrt{1+y_0(s+1)-y_0(s-1)}\mathrm{d}s.$$

$$\mathrm{Re}(\chi )=\pm \frac{\mathrm{i}\sqrt{2}}{c_{ϵ}}{\int }_{{\mathcal{C}}_1}\sqrt{1+y_0(s+1)-y_0(s-1)}\mathrm{d}s,\mathrm{Im}(\chi )=\pm \frac{\sqrt{2}}{c_{ϵ}}{\int }_{{\mathcal{C}}_2}\sqrt{1+y_0(s+1)-y_0(s-1)}\mathrm{d}s.$$

$$c_{ϵ}^2({\chi }^{\prime }{)}^2\sim -4(\xi -{\xi }_s{)}^{-1}\mathrm{as}\xi \to {\xi }_s.$$

$$\chi \sim \frac{4\mathrm{i}}{c_{ϵ}}(\xi -{\xi }_s{)}^{1/2}\mathrm{as}\xi \to {\xi }_s.$$

$$2{\chi }^{\prime }{Z}^{\prime }+{\chi }^{\prime \prime }Z=0.$$

Full text

Keywords: solitary waves, exponential asymptotics, nanoptera, FPUT lattice

AMS Subject Classifications: 34E15, 35Q51, 34C15, 37K60

1 Introduction

In this study, we demonstrate the important role played by Stokes phenomenon in the behaviour of travelling waves propagating through a diatomic Fermi-Pasta-Ulam-Tsingou (FPUT) lattice. The classical FPUT lattice contains an infinitely long line of masses connected to their neighbours by identical springs. After non-dimensionalization, this system can be represented by

[TABLE]

where $x(j,t)$ represents the position of the $j$th particle at time $t$, a dot refers to differentiation with respect to time, and $F$ is an interaction force given by

[TABLE]

A diatomic lattice has $m(j)=m_{1}$ if $j$ is odd, and $m(j)=m_{2}$ if $j$ is even, where $m_{1}\neq m_{2}$. This system has a speed of sound given by $c_{0}=\sqrt{2/(m_{1}+m_{2})}$. In this study, we are interested in the behaviour of supersonic travelling waves (with speed $c_{\epsilon}>c_{0}$) in diatomic systems with small mass ratio. We therefore set $\delta$ such that $m_{2}/m_{1}=\delta^{2}$ with $0<\delta\ll 1$.

The behaviour of travelling waves in monoatomic FPUT lattice has been treated comprehensively in, for example, [23, 24, 25]. These studies determined that the travelling wave solution is a localized solitary wave that is a regular perturbation away from a solution to the Kortweg-de Vries (KdV) equation in the long-wave asymptotic limit. These studies led to the natural question of whether the results could be generalized to polyatomic FPUT lattices, and how travelling waves would propagate through such systems.

In a later analysis, Gaison et al. [26] found that the equations governing travelling waves in diatomic FPUT systems can be written as a singular, rather than regular, perturbation to a solution to the KdV equation. The KdV solution is a long-wave approximation that depends on a small parameter, and the solution is accurate up to algebraic order in this parameter. We will later use this KdV solution as a leading-order solution to the diatomic FPUT lattice equation (1) in the small mass ratio limit. A number of other rigorous and formal studies of the behaviour of periodic lattice systems have been undertaken in the past, including [17, 43, 44, 46, 53].

Faver & Wright [21] and Hoffman & Wright [28] used a Beale ansatz (introduced in [3]) to rigorously prove the existence of supersonic travelling wave solutions in the diatomic FPUT lattice, which are composed of the sum of an exponentially localized solitary wave and periodic oscillations with exponentially small amplitude. Solutions of this form are known as ‘nanopteron’ solutions. Boyd introduced the term ‘nanopteron’ in [11] to describe weakly nonlocal solitary waves, which approximately satisfy the classical definition of a solitary wave with the distinction that a nanopteron wave is not exponentially localized, but rather tends to a small-amplitude oscillation on either one or both sides of the leading-order solitary wave, illustrated in Figure 1. Iooss & Kirchgässner [30] used related methods to identify nanopteron solutions in a chain of nonlinear oscillators with FPUT interactions and a potential function.

Typical nanopteron solutions contain a leading-order travelling wave, and a train of exponentially small oscillations on one or both sides of this front. When these solutions were first identified, it was commonly assumed that waves of this form were typically steady, and that these waves could propagate indefinitely in the same manner as classical solitary waves. This idea was contradicted by the work of Boyd [12], who studied a fifth-order Korteweg-de Vries (KdV) and concluded that only symmetric two-sided nanoptera could truly propagate in this manner. One-sided nanoptera, or two-sided nanoptera with unequal amplitude on either side, must necessarily radiate energy and can therefore not propagate indefinitely. This argument was made explicitly for the fifth-order KdV in [4], and the analysis therein can equally be applied to diatomic lattice systems.

This is consistent with numerical studies of one-sided nanopteron solutions in [27, 31, 50], which showed numerical examples of one-sided nanoptera in diatomic FPUT, Hertzian, and Toda chains respectively, each of which radiated energy and eventually decayed. Notably, metastable nanopteron solutions to the diatomic FPUT lattice equation were directly computed in the numerical study [27]. This study presented the numerical technique used in Section 5 to validate the formal asymptotic analysis.

One-sided nanoptera are therefore not truly steady in these particle systems, but instead decay on a slow timescale in a fashion that is not visible in the leading-order approximation to the wavefront or the exponentially small oscillations. Giardetti et al. [27] conjecture that this decay occurs in the diatomic FPUT system on a timescale that is exponentially large in the small mass ratio parameter. Hence, these solutions are not true travelling wave solutions, and are instead described as “metastable” or “quasi-stable” solutions.

The detailed analysis in this paper will compute the behaviour of one-sided metastable solutions; however, it is straightforward to adapt the analysis to obtain a form for symmetric two-sided solutions, and this will also be noted at appropriate points in the analysis. The term “nanopteron” will be used to describe both stable waves with symmetric oscillations, as well as metastable solutions with waves on one side. In the latter case, the leading-order asymptotics provide a useful approximation to the wave behaviour in the quasi-steady regime.

A number of other studies have been performed on diatomic FPUT (or related) lattices. Qin [45] proved the existence of periodic supersonic wavetrains in such a lattice, which correspond to the far-field periodic oscillations far from the travelling wave front in the present study. Faver later used a Beale ansatz to prove the existence of nanopteron solutions in FPUT lattices with 2-periodic spring potential (2), rather than the 2-periodic mass considered here.

Previous mathematical and experimental studies of lattices with periodic masses have revealed dynamics which do not arise in uniform lattices [31, 32, 39, 44, 48, 50]. In particular, diatomic lattices permit new families of solitary waves which can exist only for discrete values of the ratio between the masses of the two particles in a diatomic unit [31, 32, 50]. In [31], these ratios are described as anti-resonances. Asymptotic descriptions of these mass ratios for the diatomic Toda lattice were first obtained in [50] using a matched asymptotic expansion technique, and explained in [37] as destructive interference between two distinct exponentially small wavetrains in the solution. This study exploited the asymptotic form of the exponentially small terms to find a simplified asymptotic expression for the mass ratios that produce cancellation.

More generally, exponential asymptotic methods have been used to discuss the behaviour of a wide range of systems of with discrete dynamics; this includes both differential-difference and difference equations. An early of exponential asymptotics being applied to a discrete problem is found in [49], which uses Borel transform techniques to study exponentially small seperatrix splitting in the discrete Hénon map. A number of subsequent studies used Borel transform methods to identify exponentially localized waves in discrete systems, including the study of radiationless travelling kinks in $\phi^{4}$ equations [42] and sliding solitons in the discrete nonlinear Schrödinger (NLS) equation [41]. Radiationless travelling waves were also identified in saturable NLS systems [38], FPUT oscillator chains [30, 29] and a nonlocal sine-Gordon model for Josephson structures [2]. Recently [1] used similar exponential asymptotic techniques to the present study to compute exponentially small oscillations in soliton solutions to the discrete NLS equation, and compare these solutions to those of the continuous second-order NLS that obtained in the continuum limit; this analysis found that the localized solutions of the continuous system were not the same as the transparent points (or approximately localized solutions) of the discrete equation. Similar techniques were used in [36] to study dislocations in the discrete Frenkel-Kontorova model.

In these previous studies, exponential asymptotic techniques were used to identify localized solutions (otherwise known as “sliding solutions” or “radiationless solutions”, as in [2, 42, 41]). Each system possesses travelling leading-order solutions that, in general, correspond to solutions with exponentially small oscillatory tails. By directly calculating the behaviour of these oscillations, the authors were able to identify special parameter sets in which they disappear entirely, producing a purely localized solitary wave. This can occur either because the amplitudes are proportional to a constant (known as “Stokes constants”) which takes zero value for certain parameter choices, as in [38, 42, 41], or because there are multiple oscillatory contributions that interfere destructively for certain parameter choices, as in [1, 37]. These wave-free solutions are true solitary waves which do not decay due to radiation being lost to the oscillatory tail.

This paper uses an adaptation of the technique developed in [37] based on [16, 40] for studying travelling waves in a diatomic Toda lattice with small mass ratio. This previous work found that travelling waves in the diatomic Toda are nanoptera, and that the exponentially small oscillations appear as curves in the complex plane known as ‘Stokes curves’, discussed in Section 1.1. In the present study, we show that corresponding behaviour is found in supersonic travelling wave solutions for the diatomic FPUT lattice with small mass ratio. We will also show that the special nanopteron-free asymptotic solutions – true solitary waves – computed in [37, 50] are also present in the diatomic FPUT lattice, and are caused by destructive interference between two distinct oscillatory contributions.

In the present work, we do not obtain a rigorous existence proof, and we instead take a similar approach to [37, 50] and compare the results of our formal analysis to computational results. In this study, we use exponential asymptotic techniques developed by [16, 40] to provide a mathematical description of nanopteron solutions in the period-2 FPUT lattice.

The remainder of this paper is organized as follows. The method of exponential asymptotics is introduced in Section 1.1. The equations for diatomic FPUT lattices are shown in Section 2, the long wave approximation to the solution is provided. In Section 3, the solutions to the FPUT equations are expanded as asymptotic power series in the small mass ratio. The leading-order behaviour of this series is found in Section 3.1, and an ansatz is applied to determine the late-order behaviour in Section 3.2, with some technical details found in Appendix A. A full exponential asymptotic analysis is performed in Section 4. This involves finding the Stokes structure of the solution in 4.1, and determining the asymptotic form of the exponentially small oscillations in Section 4.2. In Section 5, the asymptotic results are compared to numerical simulations, and it is noted that nanopteron-free solutions are apparent in both the numerical and asymptotic results. In Section 6, an asymptotic expression is obtained for the special mass ratios that produce nanopteron-free solutions. The paper concludes in Section 7.

1.1 Exponential Asymptotics

We examine the asymptotic behavior of exponentially small, non-decaying waves that appear in the wake of a solitary-wave front in diatomic FPUT lattices. However, determining the behavior of terms that are exponentially small compared to the leading-order solution in the $\epsilon\rightarrow 0$ asymptotic limit is impossible using classical asymptotic series expansions, because the exponentially small contribution is necessarily smaller than any power of the small parameter $\delta$. Therefore, we apply specialized techniques, known as ‘exponential asymptotics’, to determine behavior on this scale [13, 14]. This section contains a brief outline of the method that will be used in the present study, adapted from the description in [37].

As we uncover this exponentially small behavior, we will see that the analytic continuation of the asymptotic solution includes curves known as ‘Stokes curves’ [47]. These curves are related to the behavior of exponentially small contributions to the solutions. As a Stokes curve is crossed, the exponentially small contribution experiences a smooth, rapid change in value in the neighborhood of the curve. In many problems, including the present investigation, the exponentially small contribution to the solution appears only on one side of a Stokes curve.

The central idea of exponential asymptotic methods is that a divergent asymptotic series may be truncated to approximate the exact solution. Furthermore, one can choose the truncation point to minimize the error between the approximation and the exact solution; this is known as ‘optimal truncation;. When a divergent series is truncated optimally, the approximation error is generally exponentially small in the asymptotic limit [14]. The problem may then be rescaled to directly determine this approximation error, allowing the exponentially small component of the solution to be determined in the absence of the asymptotic series itself. This idea was introduced by Berry [7, 8], and it was employed in [6, 9] to determine the position of Stokes curves, and associated switching behaviour in special functions such as the Airy function.

In the present paper, we apply an exponential asymptotic method developed by Olde Daalhuis et al. [40] for linear differential equations and extended by Chapman et al. [16] to nonlinear ordinary differential equations. Here we provide a brief outline of the process; see the above papers for a more detailed explanation of the methodology.

The first step in exponential asymptotic analysis is to express the solution as an asymptotic power series. In many singular perturbation problems, including the problem considered in the present study, the asymptotic series solution diverges. For a more detailed discussion of asymptotic series divergence, see [15, 19]. Optimally truncating these divergent series typically requires a general form for the asymptotic series coefficients, and it is frequently algebraically intractable to obtain such a general form. In practice, however, one does not require the exact form of the series coefficients. Instead, one needs only the so-called ‘late-order terms’, or asymptotic expressions for the $r$th series coefficient in the $r\rightarrow\infty$ limit.

Dingle [19] noted that the terms of the divergent asymptotic power series of a singularly perturbed system are typically obtained by repeated differentiations, and therefore diverge in a predictable factorial-over-power fashion. Noting this observation, Chapman et al. [16] proposed writing an ansatz for the late-order terms that is capable of describing this form of late-order term behavior. One can write an ansatz for the $r$th term of a divergent asymptotic series (denoted $g_{r}$) as $r\rightarrow\infty$ using the form

[TABLE]

where $G$, $\gamma$, and $\chi$ are functions that do not depend on $r$ but are free to vary with independent variables (and hence with $z$). The ‘singulant’ $\chi$ equals [math] at values of $z$, denoted by $z=z_{s}$, at which the leading-order behavior $g_{0}$ is singular. This ensures that the late-order ansatz for $g_{r}$ also has a singularity at $z=z_{s}$ and that the singularity increases in strength as $r$ increases. In practice, $\gamma$ typically takes constant value, and we will assume this to be the case in the present study. One can then use the ansatz (3) to optimally truncate an asymptotic expansion. The method developed in Olde Daalhuis et al. [40] involves substituting the resulting truncated series expression into the original problem to obtain an equation for the exponentially small remainder term.

The exponentially small contribution to the asymptotic solution that one obtains using the above method, denoted $g_{\mathrm{exp}}$, generally takes the form in the limit $\delta\rightarrow 0$ given by

[TABLE]

where the ‘Stokes multiplier’ $\mathcal{S}$ varies rapidly from [math] to a nonzero value as one crosses a Stokes curve. This behavior is known as ‘Stokes switching’, and it occurs along curves at which the switching exponential is maximally subdominant and hence where the singulant $\chi$ is real and positive [9]. The variation is smooth, and it occurs in a neighborhood of width $\mathcal{O}(\sqrt{\delta})$ that contains the Stokes curve.

This primary advantages of this approach are that it it reveals the important role played by the Stokes curves in the asymptotic solution, and it does not require the computation of terms beyond the leading-order expression to obtain the form of exponentially small correction terms. The latter of these advantages makes the technique particularly useful for the many nonlinear problems for which obtaining even these low-order correction terms is intractable. See the review article [14] or monograph [13] for more details on exponential asymptotics and their application to nonlocal solitary waves, [7, 8, 15] for examples of previous studies of exponential asymptotics, and [16, 40] for more details on the particular methodology that we apply in the present paper.

2 Diatomic FPUT Equation

We consider the diatomic FPUT lattice equations given in (1)–(2), with $m(j)=m_{1}$ for $j$ odd, and $m(j)=m_{2}$ for $j$ even. We express the mass ratio as $\delta=\sqrt{m_{2}/m_{1}}$, and set $0<\delta\ll 1$.

Long waves with small amplitudes in diatomic FPUT lattices were previously studied in [26] in terms of the offset $r(j,t)$ and the particle velocity $p(j,t)$, defined by

[TABLE]

giving the governing equations as

[TABLE]

where $m_{j}=m_{1}$ if $j$ is odd, and $m_{j}=m_{2}$ if $j$ is even. In [26], a long-wave solution was found in terms of a small parameter $\epsilon$ that specifies the amplitude and width of the travelling wave. The long-wave solution was given by

[TABLE]

in the limit that $\epsilon\rightarrow 0$, where

[TABLE]

We see that small $\epsilon$ corresponds to a long, small-amplitude travelling wave moving at slightly above the speed of sound in the system, denoted $c_{0}$. By integrating the particle velocities, we obtain

[TABLE]

This approximate travelling wave solution will play an important role in subsequent analysis.

While representing the system in terms of the offset $r$ and particle velocity $p$ is algebraically simpler than the system (1), it is convenient for the present analysis to seperate the motion of the heavy and light particles. For clarity in subsequent analysis, we will denote the position of the heavier particles as $y(j,t)$ where $j$ takes odd values, and the lighter particles as $z(j,t)$ where $j$ takes even values. It is always possible to remove the larger mass through non-dimensionalizing the system, so without any loss of generality, we set $m_{1}=1$ and $m_{2}=\delta^{2}$ in all subsequent analysis. We are interested in travelling wave solutions, and hence convert to a moving frame with velocity $c_{\epsilon}$ parameterized by $\xi$, where $\xi=j-c_{\epsilon}t$. We write $y(j,t)=y(\xi)$ and $z(j,t)=z(\xi)$, and the governing equations become

[TABLE]

where a dash denotes a derivative with respect to $\xi$.

Before considering the behaviour of this system for small $\delta$, we must determine whether the $\delta=0$ system contains exponentially small oscillations in the far field due to the discrete nature of the system itself, such as those seen for the discretized Korteweg-de Vries (KdV) equation in [33]. It is straightforward to see that when $\delta=0$, the governing equations become

[TABLE]

This shows that the position of the heavy particles is governed by a scaled monoatomic FPUT equation, while the lighter particles occupy the average position of their nearest neighbours. From previous work including [23, 24, 25], we know that travelling waves in the monoatomic FPUT lattice are localized solitary waves. Consequently, the long-wave approximation to the $\delta=0$ system does not introduce exponentially small far-field oscillations. This confirms that any far-field oscillations in solutions to the full system (10)–(11) are introduced by the singular perturbation caused by the small mass ratio parameter $\delta$.

3 Series Expansion

We expand the dependent variables as a Taylor series in $\delta$, giving

[TABLE]

where the coefficients $y_{r}$ and $z_{r}$ depend on the long-wave parameter $\epsilon$, but not the mass ratio $\delta$. We may determine equations for the coefficients by applying these series expressions to (10)–(11), and matching in the limit that $\delta\rightarrow 0$.

3.1 Leading-order series terms

At leading order in the limit that $\delta\rightarrow 0$, this gives

[TABLE]

This is equivalent to the full diatomic FPUT system (1) with $\delta=0$. Consequently, we may construct an approximate travelling wave solution by setting $m_{1}=1$ and $m_{2}=0$ in (7). This gives the long-wave solution in the limit that $\epsilon\rightarrow 0$ as

[TABLE]

where $y_{0}(\xi)$ describes the behaviour of heavy particles with $j$ odd, and $z_{0}(\xi)$ describes the behaviour of light particles with $j$ even. There is a significant difference at this stage between the present analysis of the diatomic FPUT lattice, and the analysis of the diatomic Toda lattice from [33]; in the analysis of the Toda lattice, the leading order solution is known exactly, whereas for the FPUT lattice, it is approximated up to $\mathcal{O}(\epsilon^{3/2})$. Consequently, we have two small parameters in the system, which introduces an extra source of error into the approximation.

Finally, for the purposes of subsequent analysis, it is important to observe that the expression for $z_{0}(\xi)$ is singular at $\xi=\xi_{s,N,\pm}$ where

[TABLE]

and that the singular behaviour is given by

[TABLE]

for any choice of $N\in\mathbb{Z}$ and sign.

We need only consider the contributions associated with singularities nearest to the real axis, or $N=1$ and $N=0$. In all subsequent analysis, we will consider the exponentially small oscillations cased by the singularity in $z_{0}(\xi)$ at $\xi=\xi_{s,1,-}$. For simplicity of notation, we will denote this particular choice as $\xi=\xi_{s}$ in subsequent analysis, giving $\xi_{s}={\pi\mathrm{i}}/{\sqrt{6}\epsilon}-1$.

The contributions associated with singularities at $\xi={\xi}_{s,0,-}$, $\xi_{s,1,+}$, and ${\xi}_{s,0,+}$ will also contribute to the asymptotic form of the far-field oscillations, and the corresponding results will be stated after the conclusion of the detailed analysis for $\xi=\xi_{s}$. The four singularities are illustrated in Figure 2.

3.2 Late-order series terms

The recursion relation is obtained by applying the series expression (14) to the governing equations (10)–(11) and matching orders of $\delta$, giving

[TABLE]

where the omitted terms are products containing $y_{r-k}$ and $z_{r-k}$ with $k>1$. These terms will be smaller than those retained in the limit that $r\rightarrow\infty$, which contain $y_{r}$, $z_{r}$, $y_{r-1}$ and $z_{r-1}$. This omission is discussed in more detail at the end of this section, after the late-order ansatz defined in (22) has been applied.

In principle this recursion relation could be repeatedly applied in order to obtain terms in the series up to arbitrarily large values of $r$, given all previous terms in the series. This process would involve solving the algebraic equation (21) for $z_{r}$, and then solving the differential-difference equation (3.2) to obtain $y_{r}$,. Repeating this process allows for the calculation of series terms up to arbitrary order; however, this is technically challenging, and it will not reveal the presence of exponentially small oscillations in the far-field, as such oscillations are typically exponentially small in the singularly perturbed limit ($\delta\rightarrow 0$).

Instead, we must follow [37] and determine the asymptotic form of the series terms in the limit that $r\rightarrow\infty$, known as the late-order terms. We see that obtaining the value of $y_{r}$ and $z_{r}$ for large $r$ requires two differentiations of the term $z_{r-1}$. Hence, we can follow the process devised in Chapman et al. [16] and apply a factorial-over-power late-order ansatz to approximate these terms as $r\rightarrow\infty$.

The required ansatz takes the form

[TABLE]

where we assume that $\alpha$ and $\beta$ take constant value.222This assumption may be omitted, as it was in [16]. In this case, the late-order equation derived by substituting the ansatz into the governing equation, given in this study in (23), contains extra terms that are smaller than the leading-order behaviour in (23) as $r\rightarrow\infty$, but larger than the first correction terms. Matching at this order gives the result that $\gamma$ is constant. An identical balancing argument to the equivalent analysis for the diatomic Toda lattice from [37] gives $\alpha$ = $\beta-2$.

At this stage we note that for sufficiently large $r$, the terms of the asymptotic series (14) diverge in a “factorial-over-power” fashion. This divergence is captured by the form of the late-order ansatzes in (22), where the gamma function grows faster as $r\rightarrow\infty$ than the algebraic decay caused by the increasing power of the singulant in the denominator. Consequently, $y_{r}\gg y_{r-k}$ and $z_{r}\gg z_{r-k}$ for $k>0$ as $r\rightarrow\infty$, with the size being controlled by the argument of the gamma function for asymptotic matching purposes.

Substituting the late-order terms (22) into the recursion relation (21) and keeping only the first two orders in the large $r$ limit gives

[TABLE]

where we use the leading-order relationship between $y_{0}(\xi)$ and $z_{0}(\xi)$, and the omitted terms are smaller than those retained as $r\rightarrow\infty$. The leading order terms in this expression grow as $\Gamma(2r+\alpha)/{\chi^{2r+\alpha}}$ in the limit that $r\rightarrow\infty$, and are therefore the same size in this limit as $z_{r+1}$ (noting that $\beta=\alpha-2$). These factorial-over-power terms of this size are therefore $\mathcal{O}(z_{r+1})$ in the limit that $r\rightarrow\infty$.

The first correction terms in this expression grow as $\Gamma(2r+\alpha-1)/{\chi^{2r+\alpha-1}}$ in the limit that $r\rightarrow\infty$, and are therefore smaller in this limit as than $z_{r+1}$, but larger than $z_{r}$. In fact, these terms correspond to the size the late-order ansatz would take for $z_{r+1/2}$, by straightforward substitution. These factorial-over-power terms of this size are therefore $\mathcal{O}(z_{r+1/2})$ in the limit that $r\rightarrow\infty$, even though this is a purely algebraic construction, rather than a term in the asymptotic series (14).

The terms that were omitted from (23) all grow as $\Gamma(2r+\alpha-4)/{\chi^{2r+\alpha-4}}$ in the limit that $r\rightarrow\infty$. Using the notation from before, they are are $\mathcal{O}(z_{r-1})$, and therefore smaller than the terms that were retained in the limit that $r\rightarrow\infty$.

3.2.1 Calculating $\chi$

Matching (23) at $\mathcal{O}(z_{r+1})$ as $r\rightarrow\infty$, we find that

[TABLE]

Integrating this, and recalling that $\chi=0$ at the singularity location $\xi=\xi_{s}$ gives

[TABLE]

The integral contour is depicted in Figure 3. While any contour may be chosen, it is helpful to divide the contour into a vertical component $\mathcal{C}_{1}$ and a horizontal component $\mathcal{C}_{2}$. The real contribution to this integral arises by integrating down $\mathcal{C}_{1}$, while the imaginary contribution is caused by integrating along $\mathcal{C}_{2}$.

Consequently, this allows us to obtain the real and imaginary parts of the singulant,

[TABLE]

We note that $\mathrm{Re}(\chi)$ is constant for real $\xi$, which causes the far field oscillations to have constant amplitude.

For Stokes switching to occur, we recall from Section 1.1 that $\mathrm{Re}(\chi)>0$, which corresponds to the positive choice of sign. Hence, we restrict our subsequent analysis to this sign choice.

It will be important in the subsequent analysis to know the local behaviour of the singulant near the singular point. It is straightforward to show by direct computation that

[TABLE]

The explicit local expansions for $y_{0}(\xi+1)$ and $y_{0}(\xi-1)$ are given in (59) and (60). The dominant contribution to the right-hand side of (27) comes from the leading-order local behaviour of $y_{0}(\xi+1)$, which is independent of $\epsilon$. This local expression may now be rearranged to give

[TABLE]

3.2.2 Calculating $Z$ and $\alpha$

Matching (23) at $\mathcal{O}(z_{r-1/2})$, after some simplification,

[TABLE]

This expression may be integrated to give

[TABLE]

where $\Lambda$ is an arbitrary constant that one may determine by considering an inner expansion of the solution in the neighbourhood of $\xi=\xi_{s}$, and matching the outer limit of this expansion with the inner limit of the late-order ansatz (22). Performing this inner analysis requires the value of $\alpha$ from (22).

To determine $\alpha$, we use (28) and (30) to find a local expression for $Z$. Direct computation gives

[TABLE]

By comparing this expression with the local form of the late-order ansatz, we see that

[TABLE]

In order for this equation to be consistent with the initial behaviour, which has singularity strength of one at $\xi=\xi_{s}$, we require that $\alpha/2-1/4=1$, which gives $\alpha=5/2$.

Using this value of $\alpha$, as well as local expressions for $\chi$ and $Z$ from (28) and (31) respectively, we may perform a local analysis to determine $z_{r}$ in the neighbourhood of $\xi=\xi_{s}$, and apply asymptotic matching in order to determine the unknown $\Lambda$. The detailed local analysis is presented in Appendix A, and yields $\Lambda={8\sqrt{\pi\mathrm{i}}}/{c_{\epsilon}}$.

We have now fully determined the late-order ansatz form for $z_{r}$, given in (22), and know that

[TABLE]

where $\chi$ is given by (25) with the positive choice of sign.

4 Exponential Asymptotics

The main idea of exponential asymptotics is that one can truncate a divergent asymptotic series optimally at some term number $N_{\mathrm{opt}}$ such that the remainder is exponentially small in size. The truncated series can be expressed as

[TABLE]

where $y_{\mathrm{exp}}$ and $z_{\mathrm{exp}}$ are exponentially small in the limit $\delta\rightarrow 0$.

4.1 Stokes Structure

We recall that there are four exponentially small asymptotic contributions, associated with the singularities depicted in Figure 2. We denote the corresponding singulants as $\chi_{1,-}$, $\chi_{0,-}$, $\chi_{1,+}$ and $\chi_{0,+}$.

Recall from Section 1.1 that Stokes curves are curves in the complex plane that correspond to $\mathrm{Im}(\chi)=0$ and $\mathrm{Re}(\chi)>0$. This means that the Stokes structure may be fully determined from the form of (25), as well as the corresponding singulant equation for the three remaining singularity contributions of interest. In each case, the Stokes curves extend vertically from the singularities of $z_{0}(\xi)$. This is depicted in Figure 4.

We prescribe that the behaviour preceding the wave front is undisturbed, and hence conclude that $z_{\mathrm{exp}}$ is zero as $\xi\rightarrow-\infty$. We therefore conclude that the four exponentially small contributions to $z_{\mathrm{exp}}$ are present to the right of the corresponding Stokes curves, and therefore are only present in the wake of the leading-order travelling wave.

4.2 Remainder Calculations

Finally, we may determine the form of the contributions that are switched on as the Stokes curves are crossed. In order to accomplish this, we truncate the asymptotic series after $N$ terms, giving

[TABLE]

where $S_{N}$ and $R_{N}$ are the remainders obtained by truncating the series after $N$ terms. If $N$ is chosen optimally, these remainder terms are exponentially small [7]. We may then analyse these remainder terms in a neighbourhood of the Stokes curves in order to determine the exponentially small asymptotic contributions to the solution behaviour obtained as the Stokes curves are crossed.

In many places, the subsequent analysis is nearly identical to [37], and hence we omit some intermediate technical details. The reader may refer to [37] for a more detailed presentation of the asymptotic analysis contained in this section, or more general descriptions of the technique described in [16, 40].

To optimally truncate the series, we apply the heuristic described in [14], and determine the optimal truncation point by finding the point at which consecutive terms in the series are equal in size. This heuristic gives $N_{\mathrm{opt}}={|\chi|}/{2{\delta}}+\omega$, where we choose $\omega\in[0,1)$ in a way that ensures that $N_{\mathrm{opt}}$ is an integer. We note that $N_{\mathrm{opt}}\rightarrow\infty$ as $\delta\rightarrow 0$, as expected.

We now apply the truncated series expression to the governing equation. Using the recursion relation to simplify the result, we obtain as $\delta\rightarrow 0$ that

[TABLE]

where the terms that we omitted are smaller in magnitude than those that we retained by a factor of $\delta$ or more in the limit that $\delta\rightarrow 0$. Throughout the remainder of this section, the asymptotic limit under consideration is $\delta\rightarrow 0$, which will be omitted for simplicity.

Assuming the truncation point occurs after a sufficiently large number of terms, we may apply the expression for the singulant and the late-order term ansatz to (37) to obtain

[TABLE]

Outside of a region in the neighbourhood of the Stokes curve, the right-hand side of this expression is exponentially small. A Green-Liouville (or WKBJ) analysis outside of this neighbourhood shows that away from Stokes curves, the remainder takes the form $R_{N}\sim CV\mathrm{e}^{-\chi/{\delta}}$ as $\delta\rightarrow 0$, where $C$ is an arbitrary constant. In order to capture the variation in the neighbourhood of the Stokes curve, we write

[TABLE]

where $A(\xi)$ is a Stokes switching parameter that varies rapidly near the Stokes curve. Applying this expression to (38) and simplifying gives

[TABLE]

We write this using $\chi$ as an independent variable. Noting that $A^{\prime}(\xi)=\chi^{\prime}(\xi)\frac{\,\mathrm{d}A}{\,\mathrm{d}\chi}$ and rearranging gives

[TABLE]

We define polar coordinates ($\chi=r\mathrm{e}^{\mathrm{i}\theta}$) and consider only variation that occurs in the angular direction. Applying the optimal truncation $N=N_{\mathrm{opt}}$, and using an asymptotic expansion for the gamma function [20] subsequently gives

[TABLE]

The right-hand side of this expression is exponentially small in $\delta$, except in the neighborhood of $\theta=0$. We therefore define an inner region $\theta=\delta^{1/2}\bar{\theta}$ and thereby find that

[TABLE]

Consequently, by integration, we see that the behavior as the Stokes curve is crossed is

[TABLE]

Therefore, converting back to outer coordinates, as the Stokes curve is crossed from ${\theta}<0$ to ${\theta}>0$, we find that $\left[A\right]_{-}^{+}\sim\pi\mathrm{i}\delta^{3/2}$ as $\delta\rightarrow 0$, where we use the notation $[f]_{-}^{+}$ to describe the change in a function $f$ as the Stokes curve is crossed from $\theta<0$ to $\theta>0$. Because $R_{N}(\xi)$ must be zero ahead of the solitary wave (i.e., when $\mathrm{Im}(\chi)>0$), we find from (39) that behind the wave front, the exponentially small contribution is given in $\delta\rightarrow 0$ limit by

[TABLE]

where $\mathcal{S}$ varies rapidly from [math] to $1$ in the neighborhood of the Stokes curve as it is crossed.

We recall that there are four relevant exponentially small contributions, associated with each of the four singularities in Figure 2. The contributions associated with $\chi_{0,-}$ and $\chi_{0,+}$ are the complex conjugates of those associated with $\chi_{1,-}$ and $\chi_{1,+}$ respectively. We consequently find that the exponentially small contribution to the asymptotic behaviour in the wake of the leading-order solitary wave is given by

[TABLE]

where $\mathcal{S}_{1}$ varies from zero to one in the neighbourhood of the Stokes curve that follows $\mathrm{Re}(\xi)=-1$ as it is crossed from right to left, $\mathcal{S}_{2}$ varies from zero to one in the neighbourhood of the Stokes curve that follows $\mathrm{Re}(\xi)=-1$ as it is crossed from right to left, and c.c. denotes the complex conjugate. The behaviour of these coefficients is illustrated in Figure 5.

Noting that $\chi_{1,-}^{\prime}(\xi)=\chi_{1,+}^{\prime}(\xi)$, we simplify this to give

[TABLE]

where $\chi_{1,-}$ is given by the integral (25), and can be approximated in the long-wave limit using the value of $y_{0}(\xi)$ given in (17). The value of $\chi_{1,+}$ can be written as an integral and approximated in near-identical fashion.

Using the remainder equations (36)–(37), we can also show that

[TABLE]

Consequently, we have obtained an asymptotic description for the exponentially small far-field oscillations present in the travelling wave solution to (1) in the small mass ratio limit.

It is very straightforward to obtain the form of symmetric two-sided nanopteron solutions, which are steady rather than simply metastable. In this case, the exponentially small contributions are given by (47)–(48), with the change that $\mathcal{S}_{j}$ switches from $-1/2$ to $1/2$ as the corresponding Stokes lines are crossed, rather than zero to one. This produces symmetric oscillations on both sides of the wavefront, leading to a system in which the wave has both an energy source ahead of the wavefront, and an energy sink in its wake. These are asymptotic representations of true travelling waves that do not radiate energy as they propagate.

4.3 Far-field oscillations

It is possible to use (24) to simplify the form of $z_{\mathrm{exp}}$ in the far field, corresponding to $\xi\rightarrow-\infty$. We note that that $y_{0}(\xi+1)-y_{0}(\xi-1)$ decays exponentially in this limit, as the leading-order solution is known to be a true solitary wave. This implies that $\chi^{\prime}\sim{\mathrm{i}\sqrt{2}}/{c_{\epsilon}}$ as $\xi\rightarrow-\infty$.

Additionally, we see from contour deformation that $\mathrm{Re}(\chi_{1,-})=\mathrm{Re}(\chi_{1,+})$, and that this quantity is constant for both singulants. This is determined by deforming the contour for $\chi_{1,+}$ to that depicted in Figure 6. Integrating $\chi^{\prime}$ along contours $\mathcal{C}_{0}$ and $\mathcal{C}_{2}$ produces imaginary contributions, while integrating along $\mathcal{C}_{1}$ produces the same real contribution as $\chi_{1,-}$. This can be confirmed by comparing the contours in Figure 3 and Figure 6, and noting that the integrand is identical in both cases.

Using these simplifications, we find that in the far-field where both Stokes multipliers are active, the asymptotic behaviour of the system simplifies to

[TABLE]

Using the result that $\Lambda=8\sqrt{\pi\mathrm{i}}/c_{\epsilon}$ and simplifying gives

[TABLE]

Recall that $\chi^{\prime}\sim\mathrm{i}\sqrt{2}/c_{\epsilon}$ as $\xi\rightarrow-\infty$. We therefore write (50) in the convenient form

[TABLE]

where $\phi_{1,\pm}$ is a constant phase offset. It is straightforward to find a similar expression for $y_{\mathrm{exp}}$ in near-identical fashion.

We see that $z_{\mathrm{exp}}$ contains two non-decaying wavetrains as $\xi\rightarrow-\infty$, far behind the wave front. These wave trains have high frequency and identical exponentially small amplitude in the limit that $\delta\rightarrow 0$. Calculating the phase offsets of the wavetrains requires direct evaluation of the integral (25), and plays an important role in finding the orthogonality condition in Section 6. Finding the orthogonality condition amounts to determining mass ratios that produce phase offsets for which the two wavetrains interfere destructively, eliminating the far-field oscillations.

In the case of steady symmetric two-sided nanoptera, $\mathcal{S}_{j}$ switches from $-1/2$ to $1/2$ as the Stokes lines are crossed. This means that the amplitude of the corresponding waves is half of the far-field amplitude for the metastable one-sided nanopteron, given in (51), and that the waves extend symmetrically in both directions.

5 Numerical comparisons

In order to determine the utility of the asymptotic description of far-field oscillations in (47), we compare the amplitude predicted by the asymptotic analysis with a numerical study of the diatomic FPUT travelling wave. While it may seem sensible to simulate the steady two-sided nanoptera, these waves do not appear in simulations with localized initial data. It is possible to produce such waves with carefully constructed numerical schemes, as in [12], in which symmetric oscillatory basis functions are used to capture the two-sided oscillations.

Fortunately the radiation of energy in the metastable one-sided system occurs on a sufficiently slow slow time-scale that the transient effects vanish substantially before the wave decay has any measurable effect on the amplitude of the oscillations. This time scale is conjectured in [27] to be exponentially large in the small mass ratio parameter. Consequently, we compute the one-sided metastable nanopteron solution with localized initial data, and determine the amplitude in this “quasi-steady” phase, in which transient effects have disappeared, but there is not yet any apparent decay in the amplitude of the oscillations.

The numerical method was implemented in Matlab using an implementation of the fourth-order Runge Kutta algorithm (RK4). Rather than computing the particle positions directly, the implementation computed the quantities $r(n,t)$ and $p(n,t)$ described in (5). The particle positions $y(n,t)$ and $z(n,t)$ were obtained by inverting the relationship between these quantities. This implementation was chosen because the leading-order behaviours $r_{0}(n,t)$ and $p_{0}(n,t)$ are both zero in the far field, which is computationally convenient. The initial conditions were chosen to be $r(n,0)=r_{0}(n,0)$ and $p(n,0)=p_{0}(n,0)$.

The domain was restricted to include $M=2^{12}$ particles with indices given by $-M/2+1\leq n\leq M/2$, with periodic boundary conditions. The initial condition was given by the leading-order travelling wave solution (17), and the time step was chosen to be $h=1/40$. In order to prevent interactions between the far-field oscillations and the leading-order travelling wave as it returns to its original position, a window function was applied to the solution. This involved multiplying $r(n,t)$ and $p(n,t)$ by a function $W(n-n_{\mathrm{front}}+M/8)$, where

[TABLE]

It is important to note that this windowing process necessarily causes the computed solution to not conserve energy. A more robust scheme would involve performing simulations on a very large domain, in which the edges are sufficiently far from the wave front that even a disturbance travelling at the speed of sound in the system, $c_{0}$, would be unable to reach the travelling wave within the simulated timespan. Through direct comparison with such systems, it is possible to conclude that windowing does not measurably affect the amplitude of the exponentially small oscillations.

The rationale that windowing will not disturb the main wave or its oscillations was given in [27], which cites the study of travelling waves in the FPUT lattice given in [24]. This study showed that a localized perturbation to a travelling wave in a monoatomic FPUT lattice generates radiation that travels slower than the wave itself. If one interprets the trailing window edge as an energy sink, this indicates that any effects in the solution caused by this energy sink should travel slower than the wave, and hence the edge of the windowed region.

The argument from [24] applies only to monoatomic lattices. It is possible that this argument be adapted to describe the effect of localized perturbations in diatomic chains; however, it is beyond the scope of this study. Instead, we rely on direct comparison between representative windowed and un-windowed simulations, which showed no discernible difference in amplitude during the quasi-steady phase when windowing was applied to the system.

This windowing applied in the numerical algorithm causes the system to fail to conserve energy; however, this did not cause apparent problems for the algorithm, which was still able to obtain the exponentially small oscillations after transient effects in the system dissipated. In order to eliminate transient effects, the computations were run for $0<t<t_{\mathrm{max}}$, where typical values of the maximum time $t_{\mathrm{max}}$ were chosen to be between $10^{5}$ and $10^{6}$, depending on the time required for transient effects to vanish. In each case, there was no measurable decay in the oscillation amplitudes for $t\leq t_{\mathrm{max}}$.

The results of this numerical computation are illustrated in Figure 7, in which they are compared to the asymptotic prediction (47). It was numerically challenging to determine wave amplitudes for values of $\delta$ smaller than those presented in Figure 7. We see that the asymptotic method does capture the qualitative behaviour predicted by the numerical simulations in the region available for comparison.

While it is challenging to extend the computations to smaller values of $\delta$, it is hoped that smaller values of $\delta$ would correspond to more accurate approximations, as in the case of the diatomic Toda lattice [37]; however, this could not be shown conclusively without a more accurate numerical study. It can be seen from the larger values of $\delta$ represented in Figure 7 (particular for $\epsilon=0.25$ and $\epsilon=0.3$) that the asymptotic approximation becomes qualitatively and quantitatively inaccurate for values of $\delta$ that are not particularly small, as is expected.

There is error introduced into the asymptotic solution by the fact that the leading-order solution (9) is an approximation that depends on a small parameter $\epsilon$, and we expect this to introduce an additional source of error to the asymptotic expansion. Naively, we would expect this error to vary straightforwardly in $\epsilon$, decreasing as $\epsilon\rightarrow 0$; however, from comparing the asymptotic and numeric results, it appears there is exists more complicated nonlinear interaction between the asymptotic error in $\epsilon$ and $\delta$.

Importantly, we see that both the numerical and asymptotic results predict the existence of mass ratios that cause the leading-order oscillations to cancel entirely, leading to truly localized solitary wave solutions. This corresponds to predictions made in [37] about the diatomic Toda lattice. We will study these points in more detail in Section 6.

Finally, by looking at individual numerical solutions, we see that sampling the very high frequency sinusoidal term at regular intervals often introduces slower periodic effects into the wave train that are not obvious from the asymptotic form given in (47). This is shown schematically in Figure 8 (a), representing an example $r(j,t)$ sampled at a particular time. In this schematic, the solid gray curve represents the full solution $r(\xi)=r(j-c_{\epsilon}t)$, where $\xi$ is a continuous variable, while the filled circles represent the sampled points for integer values of $j$ at fixed time $t=t_{s}$. The fast sinusoidally varying oscillations described in (47) are sampled by the discrete chain at regular intervals; as the sampling period is generally not an exact multiple of the sinusoid period, this can introduce slower periodic effects into the solution related to the interaction between the two periodicities.

The numerical wave train in Figure 8 (b) illustrates visible periodic effects introduced by sampling, while the wave train in Figure 8 (c) contains more complicated periodic effects. In both cases, it is still possible to calculate the amplitude of the underlying continuous waveform, although care must be taken in order to measure this quantity numerically. This sampling also explains the relationship between the results obtained in this study and the periodic effects seen in the numerically-computed wave trains in [27].

6 Orthogonality Condition

It is apparent in Figure 7 that there are particular values of $\delta$ which cause the far-field oscillation amplitude to vanish, due to destructive interference between the exponentially small oscillations. As the two wave trains in $z_{\mathrm{exp}}$ given in (47) have identical amplitude, it is possible to select parameter values so that the waves precisely cancel each other out. Such configurations produce genuinely localized solitary waves, even at exponentially small orders. Systems with these particular parameter values are said (for example, in [31, 32]) to satisfy an ‘anti-resonance’ condition. If this condition is satisfied, there are no oscillations generated in the wake of the leading order travelling wave. The wave therefore does not lose energy to an oscillatory wave train, and instead propagates without decaying.

This behaviour was seen in the diatomic Toda lattice [37], and it is possible to show that corresponding behaviour is present in asymptotic solutions the diatomic FPUT lattice. We note that in the wake of the wave front, where all oscillatory contributions are present, the exponentially small contribution is given by

[TABLE]

It is clear that destructive interference occurs if

[TABLE]

We refer to this condition as an ‘orthogonality condition’. We can simplify this condition by writing the singulant as

[TABLE]

This is a consequence of comparing the contour for $\chi_{1,+}$ shown in Figure 6 with the contour for $\chi_{1,-}$ shown in Figure 3. The orthogonality condition (54) now simplifies to give

[TABLE]

This gives the condition $\lambda/\delta=(2K+1)\pi\mathrm{i}$, for $K\in\mathbb{Z}$, or

[TABLE]

which can be approximated in the long-wave limit as

[TABLE]

While this integral does not have a convenient solution like the corresponding integral in the Toda problem in [37], it is straightforward to show that the integrand is real along the integral contour between the two singular points, and consequently that the result of the integral is real-valued. This therefore gives a set of real values for $\delta$ at which the far-field oscillations vanish.

In Table 1, asymptotic predictions of the values of $\delta$ that cause the far-field oscillations to cancel are compared with numerical predictions. These numerical predictions were found using the methods from Section 5, however the size of the time steps were decreased in size, and the maximum simulation time was increased in order to resolve the waves at smaller values of $\delta$ than those presented in Figure 7. This allowed us to determine the values of $\delta$ corresponding to $K=2,\ldots,6$.

It is clear from Table 1 that the asymptotic predictions of nanopteron-free values of $\delta$ may be accurately predicted using the asymptotic form of $z_{\mathrm{exp}}$. These values correspond to localized solitary waves that propagate without decaying due to energy radiating into the exponentially small oscillations in the wake of the wave front. Furthermore, we note that the analysis which produced these values does not depend on whether the original nanopteron solution under consideration is one-sided (and hence only metastable) or two-sided; both wave families produce truly localized waves for these parameter values.

7 Discussion and Conclusions

In this study, we considered travelling slightly supersonic wave behaviour in a diatomic FPUT lattice with small mass ratio. The asymptotic solutions of this system are nanoptera, or solitary-wave solutions that are not exponentially localized, but instead possess trains of oscillations in the far field behind the wave front. These oscillations are exponentially small, so their dynamics are invisible to ordinary asymptotic power-series approaches. The existence of these oscillations was previously proven in [21].

We demonstrated that far-field oscillations ‘switch’ across special curves in the complex plane known as Stokes curves, which originate at singularities in the analytic continuation of the leading-order behaviour. The far field oscillations present in the asymptotic wave behaviour are therefore a consequence of Stokes Phenomenon. We derived asymptotic forms for these exponentially small oscillations, given in (47) and (48), as well as a simplified expression for the behaviour of the oscillations far from the wave front 51. We compared the results of this analysis to numerical studies, and found that the asymptotic results were useful for predicting both the amplitude of the oscillations as well as the special mass ratios that produce localized solitary wave solutions.

One significant difference between the travelling wave solutions in the diatomic FPUT lattice as opposed to the previous analysis of the diatomic Toda lattice is that the leading order solution could not be determined exactly for the diatomic FPUT approximation, but was instead approximated using the long-wave limit approximation from [26]. This change introduced a second small parameter into the problem; however, with careful treatment, it was still possible to determine the leading-order behaviour of the late-order terms. This work establishes the effectiveness of exponential asymptotics for determining asymptotic nanopteron solutions in other lattice systems for which the leading-order behaviour must be approximated, such as the diatomic Hertizan lattice [44], the woodpile lattice [35], or lattices of resonant granular crystals [51].

By isolating the exponentially small terms, we found an asymptotic condition that predicted the wave trains would cancel entirely, given in (58). If this condition is satisfied, the solutions do not possess a wave train in the wake of the travelling wave front; instead, they consist of a localized solitary wave. Similar results were found using an exponential asymptotic analysis on the diatomic Toda lattice in [37], as well as other studies in which special choices of mass ratios and wave parameters produce localized solitary wave solutions [31, 34, 52]. Importantly, as in the diatomic Toda lattice from [37], this orthogonality condition arises as a consequence of Stokes Phenomenon, and the precise cancellation of exponentially small wave trains that appear as two different Stokes curves are crossed.

It is important to emphasise that the solution that we discussed in this paper is a formal solution asymptotic solution, and we did not supply a rigorous proof of these results. Faver & Wright [22] used rigorous estimates for the local leading-order solitary wave to prove the existence of the exponentially small oscillations in a nanopteron solution of FPUT systems. Typically, exponential asymptotic arguments can be made rigorous using the method of Borel transforms in order to replace the divergent tail of the asymptotic series with a quantity that can be bounded rigorously. This is beyond the scope of the present work; however, examples of rigorous exponential asymptotic bounding may be found in, for example, [5, 9, 10, 18] and others.

This paper considered one-sided nanopteron solutions in detail, while outlining the corresponding results for symmetric two-sided nanoptera. As noted, only the latter is an example of a true solitary wave. The one-sided solution is metastable, as it slowly loses energy to the far field oscillations which causes the wave to eventually decay on a long timescale. This decay is not visible in the computed asymptotic solution describing the leading-order wave or the small oscillations; a valuable direction of future study would be to determine the asymptotic order on which this decay becomes apparent. It would be of interest to determine whether the unsteady behaviour appears directly at subsequent orders of the formal asymptotic wave behaviour, or whether it requires careful application of multiple scale techniques in order to be captured asymptotically.

8 Acknowledgements

CJL thanks Prof. J. Douglas Wright for helpful discussions on the implementation of the numerical methods applied in Section 5, and Dr Justin Tzou for discussions about the manuscript. CJL is supported by ARC Discovery Project DP190101190.

Appendix A Determining $\Lambda$

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. Using Van Dyke’s matching principle we match the inner limit as $\xi\rightarrow\xi_{s}$ of the outer expansion with the outer limit of the inner expansion, determined below.

In the inner region near the singularity at $\xi=\xi_{s}$, we find that

[TABLE]

In order to locate the inner region, we must determine the region in which the inner analysis breaks down. From the form of the late-order ansatz, we find that this occurs at $\delta^{2}\chi^{-2}=\mathcal{O}(1)$ as $\delta\rightarrow 0$, or $\delta^{2}(\xi-\xi_{s})^{-1}=\mathcal{O}(1)$. This corresponds to the inner scaling $\xi-\xi_{s}=\delta^{2}\overline{\xi}$. The appropriate rescaled inner variables are given by

[TABLE]

Retaining the leading-order terms as $\delta\rightarrow 0$ in the rescaled inner equation gives

[TABLE]

We express $\hat{z}$ in terms of the local series

[TABLE]

as the term with power $j=1$ is already built into the inner form of $v(\xi)$. This gives

[TABLE]

From matching at leading order, we see that $a_{1}=-c_{\epsilon}^{2}/2$. At subsequent orders, we obtain the recurrence relation

[TABLE]

Noting the form of $a_{1}$, this gives

[TABLE]

By comparing the series expression with the inner limit of the late-order ansatz, we find using Stirling’s formula

[TABLE]