# A simple model of radiating solitary waves

J. Douglas Wright

###### Abstract.

To understand an oft-observed but poorly understood phenomenon in which a
solitary wave in a dispersive equation slowly deteriorates due a persistent
emission of radiation (i.e. a “radiating solitary wave”), we propose a bare-
bones model which captures many essential features and which we are capable of
analyzing completely by way of the Laplace transform. We find that wave
amplitude decreases at an exponential rate but with a decay constant that is
(in many cases) small beyond all orders of the frequency.

##  1\. Introduction

In the articles [27, 32, 14], simulations of solitary waves in spatially
heterogeneous variants of the Fermi-Pasta-Ulam-Tsingou (FPUT) and Toda
lattices demonstrated that such waves do not propagate without change of form
but instead continuously emit a small trailing ripple. The systems conserve
energy and consequently the solitary waves experience a commensurate and
extremely slow attenuation in amplitude. See Figure 1 for a representative
depiction.

Figure 1. An example of a “radiating solitary wave” in a diatomic FPUT
lattice. The wave is propagating to the right and the horizontal axis is in a
moving reference frame. Figure taken from [14].

A similar phenomena occurs in simulations of mass-in-mass lattices [16, 9], as
well as in a variety of models for the evolution of capillary-gravity waves
[4, 5, 33], though in this latter case the emission runs ahead of the solitary
wave. Simulations of water waves over variable bottom topography evoke a
similar dynamic [13, 25, 22, 26]. These sorts of waves, often referred to as
radiating solitary waves, are examples of metastable structures in nonlinear
dispersive systems.

The observed attenuation is so slow that one of the major tools for analyzing
the dynamics of solitary wave-like solutions, namely approximations of the
problem with the Korteweg-de Vries (KdV) (or similar) equations [12, 30, 15,
6, 29], is incapable of capturing the phenomena. Such approximations are valid
over very long but nevertheless finite time intervals; the erosion is so
subtle during the period of good approximation that it falls within the
natural error bounds. Moreover, radiating solitary waves very often occur in
problems where the construction of genuinely localized solitary waves fails
and what is found instead are generalized solitary waves (also known as a
nanopterons) [2, 1, 31, 10, 7, 17, 8, 23, 19, 20, 24]. These traveling wave
solutions are asymptotic at spatial infinity to very small amplitude co-
propagating periodic waves and are consequently of infinite energy, further
evidence that the dynamics of finite energy solitary wave-like solutions is
subtle.

Unlike their steady counterparts the nanopterons, there is at this time no
fully rigorous mathematical explanation or description of these radiating
solitary waves in any of the many problems in which they arise. There are some
careful non-rigorous treatments (especially [4]), but in the main
investigations are heuristic or numerical. The phenomenon is usually
attributed to an excitation of some sort of high-frequency oscillation by the
solitary wave, due, for instance, to band structure considerations [14], Bragg
scattering [22], internal resonances [9] or non-monotone dispersion relations.
Notably, the amplitude decays so slowly that in most cases it is very hard to
determine the rate with any level of precision and only a few articles hazard
a prediction, not all of which agree. For instance, the formal analysis in [4]
predicts the decay is exponential with a very small decay rate whereas [14]
(co-authored by the author of this paper) conjectures an algebraic rate of
decay based on numerical evidence.

To better understand radiating solitary waves, in this article we formulate a
rather bare-bones model which captures many essential features at play and
which we can rigorously analyze. For our system we are able to exactly pin
down the rate at which the solitary wave decays as well as a rather complete
description of the radiating tail. Note that our model is decidedly not meant
to quantitatively portray any of the specific systems mentioned but instead
provide a partial skeleton for the rigorous analysis of such systems down the
line.

Here are the main ingredients of our model:

  * •

A spatial variable $x\in{\bf{R}}$ and time variable $t\in{\bf{R}}$.

  * •

A solitary wave of fixed speed and profile but variable amplitude:
$a(t)q(x-t)$. Here $q(x)$ is the profile and $a(t)\in{\bf{R}}$ is the
amplitude. We assume that $q(x)$ is real-valued and satisfies some
decay/regularity conditions we make precise below.

  * •

A high-frequency simple harmonic oscillator located at each point
$x\in{\bf{R}}$. We represent each oscillator by its complex amplitude
$\psi(x,t)\in{\bf{C}}$. This field of oscillators is driven by the solitary
wave in a naive way:

(1.1) $\psi_{t}=i\omega\psi+\omega a(t)q(x-t).$

In the above, the oscillators’ natural frequency is $\omega\gg 1.$

  * •

A conserved energy

(1.2) $E:={1\over 2}a(t)^{2}+{1\over
2\omega^{2}}\int_{\bf{R}}\left|\psi(x,t)\right|^{2}dx.$

The first term is proportional to the square of the $L^{2}$-norm (i.e. energy)
of the solitary wave and the second is the energy of the oscillator field.

Together (1.1) and (1.2) form a closed system for the variables $a(t)$ and
$\psi(x,t)$. Our main result states that $a(t)\approx
a_{0}e^{-\theta_{\omega}t}$ where $\theta_{\omega}$ is all but equal to111The
normalization/notation for the Fourier transform of a function
$f:{\bf{R}}\to{\bf{C}}$ we use is:
$\displaystyle{\mathfrak{F}}[f](k):=\widehat{f}(k):={1\over
2\pi}\int_{\bf{R}}e^{-ikx}f(x)dx$. $2\pi^{2}|\widehat{q}(\omega)|^{2}$. Weak
assumptions on $q(x)$ tell us that $\widehat{q}(\omega)\to 0$ as
$\omega\to\infty$ and as such the rate of decay is very slow when $\omega$ is
big. In particular, if $q(x)$ analytic then one knows that
$|\widehat{q}(\omega)|\leq c_{1}e^{-c_{2}|\omega|}$ for some positive
constants $c_{1}$ and $c_{2}$. As such the decay rate of $a(t)$ is incredibly
small, beyond all algebraic orders of $\omega$. For instance, we find that if
$\omega=10$ and $q(x)=\operatorname{sech}(x)$ then the time it will take for
the amplitude to deteriorate to $95\%$ of its original amplitude is on the
order of $10^{11}$.

We also prove several results about the asymptotics of the oscillator field,
the most salient of which is that $e^{\theta_{\omega}t}\psi(x,t)$ converges to
a scalar multiple of a specific profile function $\sigma_{\omega}(x-t)$. This
convergence is pointwise in the moving reference frame $x-t$. We have an
explicit formula for this asymptotic profile which ultimately leads to the
description

$\sigma_{\omega}(x-t)\approx iq(x-t)/\omega+\text{trailing periodic wave of
extremely small amplitude}.$

That is to say, when viewed in a frame moving along with the solitary wave the
oscillator field looks like roughly like a copy of the solitary wave plus a
tiny periodic tail, exactly the qualitative behavior seen in the systems we
are hoping to understand.

In Section 2 we carry out a number of reformulations of (1.1) and (1.2) which
put the analysis and simulation of solutions within grasp. In Section 3 we use
the Laplace transform to describe the dynamics of $a(t)$ for large times, the
specifics of which are contained in Theorem 3.2. In Section 4 we similarly
analyze the oscillatory field $\psi(x,t)$; see Theorem 4.2. Lastly, in Section
5, we describe the results of simulations of solutions of our model and
compare them with our analytical results.

Acknowledgements: The author is grateful to the NSF who funded this work under
grant DMS-2006172.

##  2\. Reformulations and reductions

###  2.1. An ODE on a Banach space

Differentiation of (1.2) with respect to time and the use of (1.1) to
eliminate $\psi_{t}$ get us:

$0=a(t)\dot{a}(t)+{1\over\omega}\int_{\bf{R}}\Re\left[\psi^{*}(x,t)(i\psi(x,t)+a(t)q(x-t))\right]dx.$

The quantity $\psi^{*}i\psi$ is purely imaginary whereas $a$ and $q$ are
purely real so we can simplify this to:

(2.1)
$\dot{a}(t)=-{1\over\omega}\int_{\bf{R}}\Re\left[\psi(x,t)\right]q(x-t)dx.$

Note that now (1.1) and (2.1) are an ODE for the variables $(\psi,a)$ on the
Banach space $L^{2}({\bf{R}};{\bf{C}})\times{\bf{R}}$. Under the rather weak
assumption that $q(x)\in L^{2}({\bf{R}};{\bf{R}})$ it is simple enough to show
that the mapping defined by the right hand sides of these equations is smooth
from $L^{2}({\bf{R}};{\bf{C}})\times{\bf{R}}$ into itself and as such Picard’s
Theorem gives the existence of solutions for short times. The conservation law
(1.2) then allows us to extend the solution for all times. In short, we can
prove:

######  Theorem 2.1.

Fix $q(x)\in L^{2}({\bf{R}};{\bf{R}})$ and $\omega>0$. Then for any
$\psi_{0}(x)\in L^{2}({\bf{R}};{\bf{C}})$ and $a_{0}\in{\bf{R}}$ there exists
unique $(\psi,a)\in C^{1}({\bf{R}};L^{2}({\bf{R}};{\bf{C}})\times{\bf{R}})$
which satisfy (1.1), (2.1), $\psi(x,0)=\psi_{0}(x)$ and $a(0)=a_{0}$. The map
carrying $(\psi_{0},a_{0})$ to $(\psi,a)$ is continuous (for $t$ in compact
sets). Additionally for all $t$

${1\over 2}a(t)^{2}+{1\over
2\omega^{2}}\int_{\bf{R}}\left|\psi(x,t)\right|^{2}dx={1\over
2}a_{0}^{2}+{1\over 2\omega^{2}}\int_{\bf{R}}\left|\psi_{0}(x)\right|^{2}dx.$

Before moving on, note that formulating our system as in (1.1) and (2.1) is a
convenient starting point for performing simulations, which we do in Section
5.

###  2.2. A scalar delay differential equation

We can eliminate the oscillatory field from the governing equations entirely.
We solve (1.1) by Duhamel’s formula:

(2.2) $\psi(x,t)=e^{i\omega
t}\psi_{0}(x)+\omega\int_{0}^{t}e^{i\omega(t-s)}a(s)q(x-s)ds.$

Putting this into (2.1):

$\dot{a}(t)=\begin{array}[t]{@{}c@{}}\underbrace{-{1\over\omega}\int_{\bf{R}}\Re\left[e^{i\omega
t}\psi_{0}(x)\right]q(x-t)dx}\\\
j_{\omega}(t)\end{array}-\int_{\bf{R}}\Re\left[\int_{0}^{t}e^{i\omega(t-s)}a(s)q(x-s)ds\right]q(x-t)dx.$

Rearranging the order of integration in the second term plus some algebra
yields

(2.3) $\dot{a}(t)=j_{\omega}(t)-\int_{0}^{t}\cos(\omega(t-s))q\star
q(t-s)a(s)ds$

where

$\displaystyle q\star q(t):=\int_{\bf{R}}q(x+t)q(x)dx$

is the autocorrelation of $q$. It is worth pointing out here that $q\star
q(t)$ is an even function of $t$.

The scalar continuous delay differential equation (2.3) is equivalent to our
original system and will be the formulation on which we do our analysis.
Nevertheless there is one more change we make that permits a short formal
analysis.

###  2.3. A scalar renewal equation

Integrating (2.3) from $0$ to $t$ gives

$a(t)=a_{0}+\begin{array}[t]{@{}c@{}}\underbrace{\int_{0}^{t}j_{\omega}(\tau)d\tau}\\\
f_{\omega}(t)\end{array}-\int_{0}^{t}\int_{0}^{\tau}\cos(\omega(\tau-s))q\star
q(\tau-s)a(s)dsd\tau.$

Exchanging the order of integration in the second term leads to

(2.4) $a(t)=a_{0}+f_{\omega}(t)+\int_{0}^{t}\phi_{\omega}(t-s)a(s)ds$

where

$\phi_{\omega}(t):=-\int_{0}^{t}\cos(\omega s)q\star q(s)ds.$

Equation (2.4) is a renewal equation of the sort studied in [3], though in
most applications of renewal equations the function $\phi_{\omega}(t)$ is
positive with finite first moment (neither of which is the case for us).

We do a quick non-rigorous analysis of (2.4) now, in the case where
$f_{\omega}(t)=0$. Since $q\star q(t)$ is even, the Fourier convolution
theorem implies $\phi_{\omega}(t)\to-2\pi^{2}|\widehat{q}(\omega)|^{2}$ as
$t\to\infty$. If we simply replace $\phi_{\omega}$ with this asymptotic state
in (2.4) we have $a(t)=\displaystyle
a_{0}-2\pi^{2}|\widehat{q}(\omega)|^{2}\int_{0}^{t}a(s)ds.$ Differentiation of
this gives $\dot{a}=-2\pi^{2}\pi^{2}|\widehat{q}(\omega)|^{2}a$ and so the
solution of this simplified equation is
$a_{0}e^{-2\pi^{2}|\widehat{q}(\omega)|^{2}t}$. This gives us our first
glimpse as to what the rate of decay for $a(t)$ is going to be. Of course this
replacement of $\phi_{\omega}$ is completely unjustified which is why we now
move on to rigorous analysis.

##  3\. The fate of $a(t)$

In this section we determine the long time asymptotics of $a(t)$ using (2.3)
as the starting point. The integral on the right hand side of that equation is
a time-domain convolution of $a(t)$ and

$k_{\omega}(t):=\cos(\omega t)q\star q(t).$

So if we apply the Laplace transform222Here, we use the notation that the
Laplace transform of a function $f:{\bf{R}}^{+}\to{\bf{C}}$ is
${\mathcal{L}}[f](z):=F(z):=\int_{0}^{\infty}e^{-zt}f(t)dt.$ to (2.3) we get

(3.1) $zA(z)-a_{0}=J_{\omega}(z)-K_{\omega}(z)A(z).$

We have made use of a number of well-known properties of the Laplace transform
here. We adhere to the usual convention that the Laplace transform of a
variable whose name is in lower case is represented by the corresponding
capital.

We isolate $A(z)$ in (3.1):

(3.2) $A(z)={a_{0}\over z+K_{\omega}(z)}+{J_{\omega}(z)\over
z+K_{\omega}(z)}.$

From this we ascertain that the key to understanding the evolution of $a(t)$
is the function $v_{\omega}(t)$ whose Laplace transform is

(3.3) $V_{\omega}(z):={1\over z+K_{\omega}(z)}.$

The following result tells us everything we need to know about
$V_{\omega}(z)$:

######  Lemma 3.1.

Suppose that, for some $\rho>0$, $e^{\rho t}q\star q(t)\in
W^{1,\infty}({\bf{R}}^{+};{\bf{R}})$. Then, for all $\omega$ sufficiently
large, $V_{\omega}(z)$ has exactly one pole in the set $\Re(z)\geq-\rho/2$.
This pole is simple and located at $z=-\theta_{\omega}\in{\bf{R}}$ where

(3.4)
$\theta_{\omega}=2\pi^{2}|\widehat{q}(\omega)|^{2}(1+{\mathcal{O}}(1/\omega)).$

###### Proof.

Preliminary estimates: For functions $f:{\bf{R}}^{+}\to{\bf{C}}$ we put:

$\|f\|_{\rho}:=\sup_{t\geq 0}e^{\rho t}|f(t)|.$

For functions with $\|f\|_{\rho}<\infty$ we have the following elementary
estimate, which holds for $\Re(z)>-\rho$:

(3.5)
$|{\mathcal{L}}[f](z)|=\left|\int_{0}^{\infty}e^{-zt}f(t)dt\right|\leq\int_{0}^{\infty}e^{-\Re(z)t}\|f\|_{\rho}e^{-\rho
t}dt\leq{1\over\Re(z)+\rho}\|f\|_{\rho}.$

If we have $\|\dot{f}(t)\|_{\rho}<\infty$ then an integration by parts in the
above gets:

(3.6) $\begin{split}|{\mathcal{L}}[f](z)|=&\left|{1\over z}f(0)+{1\over
z}\int_{0}^{\infty}e^{-zt}\dot{f}(t)dt\right|\leq{1\over|z|}\left(|f(0)|+{1\over\Re(z)+\rho}\|\dot{f}\|_{\rho}\right).\end{split}$

This again holds for $\Re(z)>-\rho$. Lastly, the famous identity
$\displaystyle{\mathcal{L}}[tf](z)=-{d\over dz}{\mathcal{L}}[f](z),$ combined
with (3.6) implies

(3.7) $\left|{d\over
dz}{\mathcal{L}}[f](z)\right|\leq{1\over|z|}{1\over\Re(z)+\rho}\left\|{d\over
dt}\left(tf\right)\right\|_{\rho}.$

There are no poles far away from the origin: Note that the poles of
$V_{\omega}(z)$ are precisely the zeros of $z+K_{\omega}(z)$ and consequently
our analysis will focus on $K_{\omega}(z)$. The assumptions on $q\star q$ and
the definition of $k_{\omega}$ imply that $K_{\omega}(z)$ is analytic in the
set $\Re(z)>-\rho$. Moreover we have $\|k_{\omega}\|_{\rho}\leq\|q\star
q\|_{\rho}<\infty$. Importantly, $\|k_{\omega}\|_{\rho}$ can be controlled
independently of $\omega$. Using (3.5) we have $|K_{\omega}(z)|\leq\|q\star
q\|_{\rho}/(\Re(z)+\rho).$ The restriction that $\Re(z)\geq-\rho/2$ tells us
that $|K_{\omega}(z)|\leq 2\|q\star q\|_{\rho}/\rho$. Therefore $z+K(z)\neq 0$
when $|z|>2\|q\star q\|_{\rho}/\rho$ and the only place to look for the zeros
is in inside the set

$U_{\rho}:=\left\\{|z|\leq 2\|q\star
q\|_{\rho}/\rho\right\\}\cap\left\\{\Re(z)\geq-\rho/2\right\\}.$

There is just one pole near the origin: We use the fact that $\cos(\omega t)$
is of high-frequency; the frequency shifting formula for the Laplace transform
tells us

(3.8) $K_{\omega}(z)={1\over 2}{\mathcal{L}}[q\star q](z+i\omega)+{1\over
2}{\mathcal{L}}[q\star q](z-i\omega).$

Using (3.6) on the right hand side we have:

(3.9) $|K_{\omega}(z)|\leq\left({1\over 2|z+i\omega|}+{1\over
2|z+i\omega|}\right)\left(q\star q(0)+{1\over\Re(z)+\rho}\left\|{d\over
dt}q\star q\right\|_{\rho}\right).$

We have assumed that $\left\|({d/dt})q\star q\right\|_{\rho}<\infty$ and this
quantity is clearly independent of $\omega$, as is $q\star q(0)$. In
$U_{\rho}$ we have $\Re(z)+\rho\geq\rho/2$ and $|z\pm
i\omega|\geq\omega-|z|\geq\omega-2\|q\star q\|_{\rho}/\rho$. All these
together tell us that there are constants $\omega_{0}>0$ and $C_{0}>0$ for
which

(3.10) $\omega\geq\omega_{0}{\text{ and }}z\in
U_{\rho}\implies|K_{\omega}(z)|\leq{C_{0}/\omega}.$

Now let $r:=\min\left\\{|z|:z\in\delta U_{\rho}\right\\}>0$, which is
independent of $\omega$. Thus the preceding estimate allows us to find
$\omega_{1}>0$ so that

$\omega>\omega_{1}{\text{ and }}z\in U_{\rho}\implies|K_{\omega}(z)|\leq r/2.$

Thus we have $|K_{\omega}(z)|<|z|$ on the boundary of $U_{\rho}$. Since $z$
and $K_{\omega}(z)$ are analytic in $U_{\rho}$ we conclude, by Rouche’s
Theorem, that $z$ and $z+K_{\omega}(z)$ have the same number of zeros (counted
with multiplicity) inside $U_{\rho}$. Which is to say $z+K_{\omega}(z)$ has
one simple zero in that set.

And so, as foretold in Lemma 3.1, we have shown that
$V_{\omega}(z)=1/(z+K_{\omega}(z))$ has exactly one simple pole (whose
location is denoted by $-\theta_{\mu}$) in the set $\Re(z)\geq-\rho/2$. What
remains is to further pin down this pole as described in (3.4).

Characterizing the pole: Since $q\star q(t)$ is even in $t$ we have

(3.11) $K_{\omega}(0)=\int_{0}^{\infty}\cos(\omega t)q\star
q(t)dt=\pi{\mathfrak{F}}[q\star q](\omega)=2\pi^{2}|\widehat{q}(\omega)|^{2}.$

The final equality is due to the convolution theorem. The regularity/decay
condition on $q\star q(t)$ implies ${\mathfrak{F}}[q\star q](\omega)\to 0$ as
$\omega\to\infty$, and thus so does $K_{\omega}(0)$.

Differentiating (3.8) we have:

$K^{\prime}_{\omega}(z)={1\over 2}{d\over dz}{\mathcal{L}}[q\star
q](z+i\omega)+{1\over 2}{d\over dz}{\mathcal{L}}[q\star q](z-i\omega).$

Then we use (3.7) to get

$|K^{\prime}_{\omega}(z)|\leq{1\over
2}\left({1\over|z+i\omega|}+{1\over|z-i\omega|}\right){1\over\Re(z)+.9\rho}\left\|{d\over
dt}\left(tq\star q\right)\right\|_{.9\rho}.$

Note that we have used the fact that $\|f\|_{\rho}<\infty$ implies
$\|tf\|_{\rho^{\prime}}<\infty$ for $\rho^{\prime}<\rho$. Much as in the run
up to (3.10), the above estimate implies the existence of $\omega_{2}>0$ and
$C_{2}>0$ such that

(3.12) $\omega\geq\omega_{2}{\text{ and }}z\in
U_{\rho}\implies|K^{\prime}_{\omega}(z)|\leq{C_{2}/\omega}.$

The end approaches. Note that $K(z)$ is real-valued if $z=s\in{\bf{R}}$. The
fundamental theorem of calculus implies that
$|K_{\omega}(s)-K_{\omega}(0)|=\left|\int_{0}^{s}K^{\prime}(s^{\prime})ds^{\prime}\right|$
and so (3.12) leads to:

$K_{\omega}(0)+(1-C_{2}/\omega)s\leq s+K_{\omega}(s)\leq
K_{\omega}(0)+(1+C_{2}/\omega)s.$

The function on the right is zero at $s=-K_{\omega}(0)/(1+C_{2}/\omega)$ and
the function on the left at $s=-K_{\omega}(0)/(1-C_{2}/\omega)$. The
intermediate value theorem then tells us that $s+K_{\omega}(s)$ has a zero
somewhere between these two points, which, because $K_{\omega}(0)\to 0$ as
$\omega\to\infty$, are both in $U_{\rho}$. Consequently that zero is
$-\theta_{\omega}$, the unique zero of $z+K_{\omega}(z)$ in $U_{\rho}$ we
found earlier. All this, and (3.11), give:

${2\pi^{2}|\widehat{q}(\omega)|^{2}\over
1+C_{2}/\omega}\leq\theta_{\omega}\leq{2\pi^{2}|\widehat{q}(\omega)|^{2}\over
1-C_{2}/\omega}.$

This is equivalent to (3.4) and the proof is complete.

∎

Lemma 3.1 tell us that the rightmost pole of $V_{\omega}(z)$ is located at
$-\theta_{\omega}$. The common wisdom that pole placement determines decay
rate tells us to expect $v_{\omega}(t)={\mathcal{L}}^{-1}[V_{\omega}](t)$ will
behave like $e^{-\theta_{\omega}t}$ as $t\to\infty$. Nevertheless inverting
the Laplace transform is a sometimes subtle business and so we make a precise
statement and proof.

######  Theorem 3.2.

Suppose that, for some $\rho>0$, $e^{\rho|x|}q(x)\in
W^{1,\infty}({\bf{R}};{\bf{R}})$. Then the the solution of (1.1) and (2.1)
with initial data $\psi(x,0)=0$ and $a(0)=a_{0}$ satisfies (for $\omega$
sufficiently large)

(3.13) $a(t)=a_{0}r_{\omega}e^{-\theta_{\omega}t}+b_{\omega}(t)$

where
$\theta_{\omega}=2\pi^{2}|\widehat{q}(\omega)|^{2}(1+{\mathcal{O}}(1/\omega))$,
$r_{\omega}=1+{\mathcal{O}}(1/\omega)$ and
$\|b_{\omega}\|_{\rho^{\prime}}<\infty$ for any $\rho^{\prime}<\rho/2$.

Additionally if $e^{\rho|x|}\psi_{0}(x)\in L^{\infty}({\bf{R}};{\bf{C}})$ then
the solution of (1.1) and (2.1) with initial data $\psi(x,0)=\psi_{0}(x)$ and
$a(0)=a_{0}$ satisfies (for $\omega$ sufficiently large)

(3.14) $a(t)=\alpha e^{-\theta_{\omega}t}+{b}_{\omega}(t)$

for some constant $\alpha=\alpha(\psi_{0},a_{0})\in{\bf{R}}$; $b_{\omega}(t)$
satisfies the same estimate as above.

###### Proof.

From Lemma 3.1 we know that the pole of $V_{\omega}(z)$ at $z=-\theta_{\mu}$
is simple. We now compute the residue at the pole in the usual way:

$r_{\omega}:=\textrm{Res}(V_{\omega},-\theta_{\omega})=\lim_{z\to-\theta_{\omega}}(z+\theta_{\omega})V_{\omega}=\lim_{z\to-\theta_{\omega}}{z+\theta_{\omega}\over
z+K_{\omega}(z)}={1\over 1+K^{\prime}(-\theta_{\omega})}.$

Estimate (3.12) tells us that $r_{\omega}=1+{\mathcal{O}}(1/\omega)$.

With this we now know that

$W_{\omega}(z):=V_{\omega}(z)-{r_{\omega}\over z+\theta_{\omega}}$

is analytic for $\Re(z)>-\rho/2$. Because of (3.9) we have

$\lim_{z\to\infty}zW_{\omega}(z)={z\over z+K_{\omega}(z)}-{r_{\omega}z\over
z+\theta_{\omega}}=1-r_{\omega}$

provided $\Re(z)>-\rho/2$. Thus we have $|W_{\omega}(z)|\leq C/|z|$ for $|z|$
big enough. This is not a rapid of enough decay to for us to deploy standard
inversion results from (for instance) [21].

To get around this we let

$B_{\omega}(z):=V_{\omega}(z)-{r_{\omega}\over
z+\theta_{\omega}}-{1-r_{\omega}\over z+2\rho}.$

The final term there is a sort of “fudge factor.” Note that $B_{\omega}(z)$ is
analytic for $\Re(z)>-\rho/2$. It is easy to see that
$\lim_{z\to\infty}zB_{\omega}(z)=0$ provided $\Re(z)>-\rho/2$. Moreover a
routine computation gives:

(3.15)
$\begin{split}&\lim_{z\to\infty}z^{2}B_{\omega}(z)=(-2r_{\omega}+2)\rho+r_{\omega}\theta_{\omega}.\end{split}$

Thus $|B_{\omega}(z)|<C/|z|^{2}$ for $|z|$ big enough. This rate is fast
enough to use Lemma 76.4 in [21] and conclude that
$e^{\rho^{\prime}t}{\mathcal{L}}^{-1}[B_{\omega}](t)\to 0$ as $t\to\infty$ for
any $\rho^{\prime}<\rho/2$.

And so all together we find that

$\begin{split}v_{\omega}(t)={\mathcal{L}}^{-1}\left[{r_{\omega}\over
z+\theta_{\omega}}+{1-r_{\omega}\over
z+2\rho}+B_{\omega}\right](t)=r_{\omega}e^{-\theta_{\omega}t}+(1-r_{\omega})e^{-2\rho
t}+{\mathcal{L}}^{-1}[B_{\omega}](t).\end{split}$

If $\psi_{0}(x,0)=0$ then from (3.2) we have $a(t)=a_{0}v_{\omega}(t)$ and the
above leads directly to (3.13).

Now for (3.14). If $\psi_{0}(x)\neq 0$ the decay conditions placed upon it
imply that $\|j_{\omega}\|_{\rho}\leq C/\omega$; here is the calculation:

$|e^{\rho t}j_{\omega}(t)|\leq{1\over\omega}e^{\rho
t}\int_{\bf{R}}|\psi_{0}(x)||q(x-t)|dx\leq{C\over\omega}e^{\rho
t}\int_{\bf{R}}e^{-\rho|x|}e^{-\rho|x-t|}dx\leq{C\over\omega}.$

This tells us that $J_{\omega}(z)$ is analytic for $\Re(z)>-\rho$. Thus the
term $J_{\omega}(z)V_{\omega}(z)$ in (3.2) is analytic in $\Re(z)>-\rho/2$
except for the pole of $V_{\omega}(z)$ at $-\theta_{\omega}$. The same sorts
of steps as above can be repeated to show that
$\displaystyle\lim_{t\to\infty}e^{\theta_{\omega}t}{\mathcal{L}}^{-1}[J_{\omega}V_{\omega}](t)$
exists, from which (3.14) follows. ∎

##  4\. The fate of $\psi(x,t)$

Now that we have determined the dynamics of $a(t)$ for large values of $t$, we
do the same for $\psi(x,t)$. We begin by observing that because $a(t)\to 0$
and we have the conservation of the energy (1.2) we know that $\psi(x,t)$ does
not converge to zero in the $L^{2}({\bf{R}};{\bf{C}})$ norm. This is not much
of a statement, but it does indicate that the ultimate behavior of $\psi(x,t)$
is not disintegration.

For a more refined analysis, our starting point is (2.2) which expresses
$\psi(x,t)$ explicitly in terms of $a(t)$. Here is our first result:

######  Corollary 4.1.

If $e^{\rho|x|}q(x)\in W^{1,\infty}({\bf{R}};{\bf{R}})$ and
$e^{\rho|x|}\psi_{0}(x)\in L^{\infty}({\bf{R}};{\bf{C}})$ then the solution of
(1.1) and (2.1) with initial data $\psi(x,0)=\psi_{0}(x)$ and $a(0)=a_{0}$
satisfies (for $\omega$ sufficiently large)

$\lim_{t\to\infty}e^{-i\omega
t}\psi(x,t)=\psi_{0}(x)+\omega\int_{0}^{\infty}e^{-i\omega s}a(s)q(x-s)ds$

pointwise in $x$.

###### Proof.

Divide both sides of (2.2) by $e^{i\omega t}$ and take the limit. The integral
converges due to the restrictions placed upon $\psi_{0}$ and $q$. ∎

The most important takeaway from this result is that it indicates that
$\psi(x,t)$ does not decay to zero (in the supremum norm) in the large time
limit. This can be make rigorous by drilling down into the integral term above
to get a more refined picture, but it turns out it is more interesting to view
$\psi(x,t)$ in a frame that moves along with the solitary wave $q(x-t)$. And
so we put $x-t=l$ and $\varphi(l,t)=\psi(x,t).$ This converts (2.2) to

(4.1) $\varphi(l,t)=e^{i\omega
t}\psi_{0}(t+l)+\omega\int_{0}^{t}e^{i(t-s)}a(s)q(t+l-s)ds.$

We begin by considering the situation where $\psi_{0}(x)=0$ and $a_{0}=1$, in
which case $\varphi(l,t)=\gamma_{\omega}(l,t)$ with

$\gamma_{\omega}(l,t):=\omega\int_{0}^{t}e^{i\omega(t-s)}v_{\omega}(s)q(t+l-s)ds.$

If we fix $l$ and take the Laplace transform of the above with respect to $t$,
the convolution and frequency shifting identities get us

$\Gamma_{\omega}(l,z)=\omega Q_{l}(z-i\omega)V_{\omega}(z)$

where

$Q_{l}(z):={\mathcal{L}}[q(\cdot+l)](z).$

With our standard assumption that $e^{\rho|x|}q(x)\in
W^{1,\infty}({\bf{R}};{\bf{R}})$ we have $\|q(\cdot+l)\|_{\rho}<\infty$ for
any $l\in{\bf{R}}$. Thus $Q_{l}(z)$ will be analytic when $\Re(z)>-\rho$. In
turn this implies that $\Gamma_{\omega}(l,z)$ will inherit the simple pole at
$z=-\theta_{\omega}$ from $V_{\omega}(z)$ and that this is the only
singularity when $\Re(z)\geq-\rho/2$. It is easy to compute that

(4.2) $\textrm{Res}(\Gamma_{\omega}(l,z),-\theta_{\omega})=\omega
r_{\omega}\begin{array}[t]{@{}c@{}}\underbrace{Q_{l}(-\theta_{\omega}-i\omega)}\\\
\sigma_{\omega}(l)\end{array}.$

Thus the common wisdom implies that $\gamma_{\omega}(l,t)$ looks like
$\sigma_{\omega}(l)e^{-\theta_{\omega}t}$ as $t\to\infty$. We codify this in
the following result, whose proof is so close to that of Theorem 3.2 we omit
it:

######  Theorem 4.2.

Suppose that, for some $\rho>0$, $e^{\rho|x|}q(x)\in
W^{1,\infty}({\bf{R}};{\bf{R}})$. Then the the solution of (1.1) and (2.1)
with initial data $\psi(x,0)=0$ and $a(0)=a_{0}$ satisfies (for $\omega$
sufficiently large)

(4.3) $\psi(t+l,t)=a_{0}\omega
r_{\omega}\sigma_{\omega}(l)e^{-\theta_{\omega}t}+\eta(l,t).$

Here $\sigma_{\omega}(l)$ is given in (4.2),
$r_{\omega}=1+{\mathcal{O}}(1/\omega)$ and
$\|\eta(l,\cdot)\|_{\rho^{\prime}}<\infty$ for all $\rho^{\prime}<\rho/2$ and
every $l\in{\bf{R}}$.

Additionally if $e^{\rho|x|}\psi_{0}(x)\in L^{\infty}({\bf{R}};{\bf{C}})$ then
the solution of (1.1) and (2.1) with initial data $\psi(x,0)=\psi_{0}(x)$ and
$a(0)=a_{0}$ satisfies (for $\omega$ sufficiently large)

(4.4) $\psi(t+l,t)=\beta\omega
r_{\omega}\sigma_{\omega}(l)e^{-\theta_{\omega}t}+\eta(l,t)$

for some finite constant $\beta=\beta(\psi_{0},a_{0})\in{\bf{R}}$;
$r_{\omega}$ and $\eta(l,t)$ satisfy the same estimates as above.

######  Remark 1.

The result above is pointwise in $l$. It is quite possible that a stronger
mode of convergence holds here, though the technical difficulty in
establishing this is motivation enough to leave that for another article.

To close out this section, we now describe the “asymptotic profile”
$\sigma_{\omega}(l)$ in greater detail. The analysis here is formal though it
could be made rigorous if we stack enough hypotheses on $q(x)$. Recall that
$\sigma_{\omega}(l)=Q_{l}(-\theta_{\omega}-i\omega).$ Since
$\theta_{\omega}\to 0$ as $\omega\to\infty$ and $Q_{l}(z)$ is analytic in $z$
we have
$Q_{l}(\theta_{\omega}-i\omega)=Q_{l}(-i\omega)+{\mathcal{O}}(\theta_{\omega}).$
Then by definition we have

$Q_{l}(-i\omega)=\int_{0}^{\infty}e^{i\omega t}q(t+l)dt=e^{-i\omega
l}\int_{l}^{\infty}e^{i\omega x}q(x)dx.$

And so we see that $e^{i\omega l}Q_{l}(-i\omega)\to 2\pi\widehat{q}(-\omega)$
as $l\to-\infty$ and $e^{i\omega l}Q_{l}(-i\omega)\to 0$ as $l\to\infty$. This
latter convergence will be exponentially fast. Lastly the same sort of
calculation that led to (3.6) gives
$Q_{l}(-i\omega)=iq(l)/\omega+{\mathcal{O}}(1/\omega^{2}).$

Putting everything together we have

$\sigma_{\omega}(l)\sim\begin{cases}2\pi\widehat{q}(-\omega)e^{-i\omega
l}&l\ll 0\\\ {i}q(l)/\omega&l\sim 0\\\ \text{exponential decay}&l\gg
0.\end{cases}$

In short, $\sigma_{\omega}(l)$ looks like a scalar multiple of the solitary
wave profile plus a trailing periodic wave of small amplitude and frequency
$\omega$.

##  5\. Simulations

We have simulated solutions of our model with a variety choices of the profile
$q(x)$ and frequency $\omega$. We always take $\psi_{0}(x)=0$ since our
analytic results indicate the effects of this part of the initial data are
transient and do not alter the long time behavior (or at least the rate of
decay of $a(t)$). Likewise $a_{0}=1$ in all cases.

Our method is straightforward: we treat the system as an ODE for $(\psi,a)$ as
in (1.1) and (2.1) and simulate using an RK4 algorithm. We implement the
integrals in (2.1) via Simpson’s rule. We also compute $\sigma_{\omega}(l)$
numerically, as for most functions closed expressions are hard to obtain. We
again use Simpson’s rule for the computation. All simulations were done in
MATLAB.

Before we get into the results we note that with our method high accuracy/long
time simulations of the problem are challenging to obtain even for modestly
large values of $\omega$. Roughly speaking, to accurately resolve a decay rate
like $e^{-\theta_{\omega}t}$ we would need to simulate out to times of
${\mathcal{O}}(1/\theta_{\omega})$ and (for the RK4 method we use) have a
temporal step size which is ${\mathcal{O}}(\theta_{\omega})$. The long time of
integration in turn implies a large spatial domain (also
${\mathcal{O}}(1/\theta_{\omega})$) is needed, since $q(x-t)$ propagates in
space. And we need to resolve that spatial domain at the same step size as the
temporal one. All these considerations tell us that we need
${\mathcal{O}}(1/\theta_{\omega}^{4})$ operations, at a bare minimum, to
experimentally determine $\theta_{\omega}$ in a quantitatively reliable
fashion. And in the most interesting cases, $\theta_{\omega}$ is exponentially
small in $\omega$, meaning that we quickly reach a computational bottleneck.

In light of these inherent difficulties we simply integrate out to $t=1000$
and choose our step size to be $\pi(10\omega)^{-1}$. For the largest values of
$\omega$, our simulations have been pushed past the point at which we can be
convinced of their quantitative reliability and instead we view them as a
qualitative illustration of our results (which are, after all, fully
justified) and the phenomena they describe.

We show our results in the figures which follow. Each figure contains the same
four sorts of graphs.

  * •

Upper left panel: A semilog plot of the numerical solution $a(t)$ vs $t$ for
the entirety of the run. In the same panel is shown the graph of $\alpha
e^{-2\pi|\widehat{q}(\omega)|^{2}t}$ where $\alpha$ is just a scaling factor
used to make the graph readable. Theorem 3.2 tells us that these should be
nearly parallel when $\omega$ is large and in fact we do see this.

  * •

Upper right panel: A plot of $a(t)$ vs $t$ during the very beginning of the
simulation. In each case we see that the $a(t)$ oscillates a few times and
quickly “settles down” into the slow decay. The time it takes to settle down
does not seem to depend on $\omega$, as predicted by Theorem 3.2.

  * •

Lower left panel: The graphs of the real and imaginary parts of $\psi$ vs $x$,
at $t=1000$, the end of the run. The figure is zoomed in on the leading edge
of the solution, located near $x=1000$.

  * •

Lower right panel: The numerically computed asymptotic profile function
$\sigma_{\omega}(l)$ vs $l$. The horizontal scale is arranged to match that of
the previous panel. Theorem 4.2 tells us that $\psi(x,t)$ should look much
like (a scalar multiple of) $\sigma_{\omega}(l)$ as $t\to\infty$ and indeed we
see exactly this.

###  5.1. $q(x)=\operatorname{sech}(x)$

We lead off with this choice for $q(x)$ because of the ubiquity of hyperbolic
secant in profiles for solitary waves. One has
$\widehat{q}(\omega)=\operatorname{sech}(\pi\omega/2)/2$ and consequently we
have

$\theta_{\omega}={\pi^{2}\over 2}\operatorname{sech}^{2}\left({\pi\omega\over
2}\right)(1+{\mathcal{O}}(1/\omega)).$

We take $\omega=2,\ 4$ and $8$. Results are shown in Figures 2-4. When
$\omega=8$ both the decay of $a(t)$ and the trailing oscillatory tails of
$\psi(x,1000)$ and $\sigma_{\omega}(l)$ have vanished to the naked eye.

###  5.2. $q(x)=e^{-x^{2}}$

We selected this because the ultra-rapid decay of its Fourier transform
renders the decay rate incredibly small even at modest values of $\omega$. We
have

$\theta_{\omega}={\pi\over 2}e^{-\omega^{2}/2}(1+{\mathcal{O}}(1/\omega)).$

We show results for $\omega=4$ and $\omega=6$, in Figures 5 and 6. Note how
there is essentially no decay even at $\omega=6$. Likewise at this value the
trailing oscillatory waves in $\psi(x,1000)$ and $\sigma_{\omega}(l)$ are
invisible.

###  5.3. $q(x)=e^{-|x|}$.

This “peakon” profile is not analytic and consequently the decay rate merely
goes to zero algebraically fast and as such the decay is more obvious at
larger values of $\omega$. To wit

$\theta_{\omega}={2\over(1+\omega^{2})^{2}}(1+{\mathcal{O}}(1/\omega)).$

We show results for $\omega=4,\ 8$ and $16$ in Figures 7-9. Both the decay and
the tails are visible at $\omega=16$.

###  5.4. $q(x)=(1-|x|)_{+}$.

The subscript “$+$” mean to take the positive part, which is to say that
$q(x)$ is the “tent” map. This is also non-analytic and its Fourier transform
is sometimes zero (unlike the others) and as such if we fine tune $\omega$ we
can get solutions which do not decay at all. Specifically we have

$\theta_{\omega}={1\over 2}\textrm{sinc}(\omega/2)^{4}$

so if we take $\omega$ to be an even multiple of $\pi$ we should see no
decay/no tails. We show results here (Figures 10-13) for $\omega=2\pi,\ 3\pi,\
4\pi$ and $5\pi$ and we see exactly this behavior.

##  6\. Conclusions, remarks and future directions

A key takeaway of this article is that in this simple model of radiating
solitary waves the rate of attenuation in the amplitude is very slow, but
nevertheless exponential. This is in line with the results of [4] and not the
algebraic decay rate the author predicted in [14]. The extremely slow rate of
decay predicted by the main results provides further evidence that radiating
solitary waves present challenging complications both on the numerical and
analytic sides. Indeed, even in this simple model really capturing the rate of
decay numerically for large values of $\omega$ would require a much more
sophisticated approach than is used here.

One can reasonably ask, however, whether or not our model here will truly be
reflective of the systems which possess radiating solitary waves. After all,
the derivation of (1.1)-(1.2) is ad hoc and (as is always the case in such
models) important features/considerations have been omitted. Indeed, in our
model only the amplitude of the solitary wave is variable and it has fixed
speed. In nearly all systems with solitary waves, the speed, amplitude and
wavelength of the solitary wave are linked. This is particularly important in
the proofs of stability of solitary waves in KdV [28] and FPUT [11] and in the
derivation of effective equations for solitary waves in potentials [18]. By
restricting to a fixed speed and width, our model is ultimately linear in its
unknowns and this in turn makes the analysis by Laplace transform possible.
Preliminary attempts at incorporating variable speed/wavelength result in
nonlinear systems and, consequently, there are substantial technical
challenges that make the application of the methods used here non-obvious. One
avenue is to adapt the refined asymptotics for stationary problems in [1] to
this time-dependent setting. Success in that venture will lead to the next,
and most exciting challenge: connecting such results rigorously to a full
system with radiating solitary waves. Work is underway.

Figure 2.  Figure 3.  Figure 4.  Figure 5.  Figure 6.  Figure 7.  Figure 8.
Figure 9.  Figure 10.  Figure 11.  Figure 12.  Figure 13.

## References

  * [1] T. R. Akylas and T.-S. Yang, On short-scale oscillatory tails of long-wave disturbances, Stud. Appl. Math., 94 (1995), pp. 1–20.
  * [2] J. T. Beale, Exact solitary water waves with capillary ripples at infinity, Comm. Pure Appl. Math., 44 (1991), pp. 211–257.
  * [3] R. Bellman and K. L. Cooke, Differential-difference equations, Academic Press, New York-London, 1963.
  * [4] E. S. Benilov, R. Grimshaw, and E. P. Kuznetsova, The generation of radiating waves in a singularly-perturbed Korteweg-de Vries equation, Phys. D, 69 (1993), pp. 270–278.
  * [5] J. L. Bona, V. A. Dougalis, and D. E. Mitsotakis, Numerical solution of Boussinesq systems of KdV-KdV type. II. Evolution of radiating solitary waves, Nonlinearity, 21 (2008), pp. 2825–2848.
  * [6] M. Chirilus-Bruckner, C. Chong, O. Prill, and G. Schneider, Rigorous description of macroscopic wave packets in infinite periodic chains of coupled oscillators by modulation equations, Discrete Contin. Dyn. Syst. Ser. S, 5 (2012), pp. 879–901.
  * [7] T. E. Faver, Nanopteron-stegoton traveling waves in spring dimer Fermi-Pasta-Ulam-Tsingou lattices, Quart. Appl. Math., 78 (2020), pp. 363–429.
  * [8] , Small mass nanopteron traveling waves in mass-in-mass lattices with cubic FPUT potential, J. Dynam. Differential Equations, 33 (2021), pp. 1711–1752.
  * [9] T. E. Faver, R. H. Goodman, and J. D. Wright, Solitary waves in mass-in-mass lattices, Z. Angew. Math. Phys., 71 (2020), pp. Paper No. 197, 20\.
  * [10] T. E. Faver and J. D. Wright, Exact diatomic Fermi-Pasta-Ulam-Tsingou solitary waves with optical band ripples at infinity, SIAM J. Math. Anal., 50 (2018), pp. 182–250.
  * [11] G. Friesecke and R. L. Pego, Solitary waves on Fermi-Pasta-Ulam lattices. IV. Proof of stability at low energy, Nonlinearity, 17 (2004), pp. 229–251.
  * [12] J. Gaison, S. Moskow, J. D. Wright, and Q. Zhang, Approximation of polyatomic FPU lattices by KdV equations, Multiscale Model. Simul., 12 (2014), pp. 953–995.
  * [13] J. Garnier, J. C. Muñoz Grajales, and A. Nachbin, Effective behavior of solitary waves over random topography, Multiscale Model. Simul., 6 (2007), pp. 995–1025.
  * [14] N. Giardetti, A. Shapiro, S. Windle, and J. D. Wright, Metastability of solitary waves in diatomic FPUT lattices, Math. Eng., 1 (2019), pp. 419–433.
  * [15] F. Hadadifard and J. D. Wright, Mass-in-mass lattices with small internal resonators, Stud. Appl. Math., 146 (2021), pp. 81–98.
  * [16] S. Hauver, X. He, D. Mei, E. G. Charalampidis, P. G. Kevrekidis, E. Kim, J. Yang, and A. Vainchtein, Lattices with internal resonator defects, Phys. Rev. E, 98 (2018), p. 032902.
  * [17] A. Hoffman and J. D. Wright, Nanopteron solutions of diatomic Fermi-Pasta-Ulam-Tsingou lattices with small mass-ratio, Phys. D, 358 (2017), pp. 33–59.
  * [18] J. Holmer, Dynamics of KdV solitons in the presence of a slowly varying potential, Int. Math. Res. Not. IMRN, (2011), pp. 5367–5397.
  * [19] M. A. Johnson and J. D. Wright, Generalized solitary waves in the gravity-capillary Whitham equation, Stud. Appl. Math., 144 (2020), pp. 102–130.
  * [20] N. Joshi and C. J. Lustri, Generalized solitary waves in a finite-difference Korteweg–de Vries equation, Stud. Appl. Math., 142 (2019), pp. 359–384.
  * [21] T. W. Körner, Fourier analysis, Cambridge University Press, Cambridge, 1988.
  * [22] Y. Liu and D. K. P. Yue, On generalized Bragg scattering of surface waves by bottom ripples, J. Fluid Mech., 356 (1998), pp. 297–326.
  * [23] E. Lombardi, Oscillatory integrals and phenomena beyond all algebraic orders, vol. 1741 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, 2000.  With applications to homoclinic orbits in reversible systems.
  * [24] C. J. Lustri and M. A. Porter, Nanoptera in a period-2 Toda chain, SIAM J. Appl. Dyn. Syst., 17 (2018), pp. 1182–1212.
  * [25] A. Nachbin and G. C. Papanicolaou, Water waves in shallow channels of rapidly varying depth, J. Fluid Mech., 241 (1992), pp. 311–332.
  * [26] O. Nakoulima, N. Zahibo, E. Pelinovsky, T. Talipova, and A. Kurkin, Solitary wave dynamics in shallow water over periodic topography, Chaos, 15 (2005), pp. 037107, 8.
  * [27] Y. Okada, S. Watanabe, and H. Tanaca, Solitary wave in periodic nonlinear lattice, Journal of the Physical Society of Japan, 59 (1990), pp. 2647–2658.
  * [28] R. L. Pego and M. I. Weinstein, Asymptotic stability of solitary waves, Comm. Math. Phys., 164 (1994), pp. 305–349.
  * [29] D. E. Pelinovsky and G. Schneider, The monoatomic FPU system as a limit of a diatomic FPU system, Appl. Math. Lett., 107 (2020), pp. 106387, 8\.
  * [30] G. Schneider and C. E. Wayne, The rigorous approximation of long-wavelength capillary-gravity waves, Arch. Ration. Mech. Anal., 162 (2002), pp. 247–285.
  * [31] S. M. Sun, Existence of a generalized solitary wave solution for water with positive Bond number less than $1/3$, J. Math. Anal. Appl., 156 (1991), pp. 471–504.
  * [32] Y. Tabata, Stable solitary wave in diatomic toda lattice, Journal of the Physical Society of Japan, 65 (1996), pp. 3689–3691.
  * [33] Y. Tan, J. Yang, and D. E. Pelinovsky, Semi-stability of embedded solitons in the general fifth-order KdV equation, Wave Motion, 36 (2002), pp. 241–255.