<figure><img src="image_1.png" /></figure>

HAL
open science

# Higher order corrections to "beyond all order" effects in a fifth order Korteweg-de Vries equation

Gyula Fodor, Péter Forgács, Muneeb Mushtaq

► To cite this version:

Gyula Fodor, Péter Forgács, Muneeb Mushtaq. Higher order corrections to "beyond all order" effects in a fifth order Korteweg-de Vries equation. Physical Review D, 2023, 107 (10), pp.105002. 10.1103/PhysRevD.107.105002. hal-04038996

HAL Id: hal-04038996

https://hal.science/hal-04038996v1

Submitted on 13 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Higher order corrections to beyond-all-order effects in a fifth order
Korteweg–de Vries equation

Gyula Fodor,¹ Péter Forgács,¹,² and Muneeb Mushtaq¹,³

¹Wigner Research Centre for Physics, 1525 Budapest 114, P.O. Box 49, Hungary

²Institut Denis-Poisson CNRS/UMR 7013, Université de Tours, Parc de Grandmont, 37200 Tours, France

³Institute for Theoretical Physics, Eötvös University,
Pázmány Péter sétány 1/A, H-1117 Budapest, Hungary

<figure><img src="image_3.png" /></figure>

(Received 3 February 2023; accepted 19 April 2023; published 8 May 2023)

A perturbative scheme is applied to calculate corrections to the leading, exponentially small (beyond-all-orders) amplitude of the “trailing” wave asymptotics of weakly localized solitons. The model considered is a Korteweg–de Vries equation modified by a fifth order derivative term, $\epsilon^2 \partial_x^5$ with $\epsilon \ll 1$ (fKdV). The leading order corrections to the tail amplitude are calculated up to $O(\epsilon^5)$. An arbitrary precision numerical code is implemented to solve the fKdV equation and to check the perturbative results. Excellent agreement is found between the numerical and analytical results. Our work also clarifies the origin of a long-standing disagreement between the $O(\epsilon^2)$ perturbative result of Grimshaw and Joshi [SIAM J. Appl. Math. **55**, 124 (1995)] and the numerical results of Boyd [Comput. Phys. **9**, 324 (1995)].

DOI: 10.1103/PhysRevD.107.105002

I. INTRODUCTION

Oscillons—slowly radiating lumps in theories contain-
ing scalar fields with a well-defined core and living for very
long times—have generated quite some interest in view of
their numerous physical applications [1,2]; for recent
reviews see [3–5] and the references cited therein. An
important problem is to determine the radiation rate and
lifetimes of oscillons; see Refs. [6,7] for recent results
and also for further references. In Refs. [4,8,9] the
perturbative scheme of Ref. [10] has been generalized to
a large class of theories, yielding the leading order estimate
for oscillon lifetimes. The main idea has been to approxi-
mate slowly radiating oscillons through the adiabatic
evolution of appropriate stationary configurations, called
“quasibreathers”—weakly localized lumps with asymptotic
standing wave tails [11,12]. Quasibreathers are time peri-
odic and can be thought of as oscillons made stationary by
incoming radiation from infinity. Determination of the
amplitude of the quasibreather wave tails is necessary to
deduce the radiation rate of time-dependent oscillons. A
rather nontrivial aspect of the perturbative computation of
the standing wave tail amplitude for quasibreathers corre-
sponding to long-lived oscillons is that it is beyond all
orders in perturbation theory.

It still remains a challenging problem to calculate higher
order corrections to the oscillating tails of quasibreathers.
To prepare the ground for such higher order computations
we have taken up a much simpler problem—the compu-
tation of higher order corrections to asymptotic wave tails
in the familiar Korteweg–de Vries (KdV) equation modi-
fied by a fifth order derivative term (fKdV), also called the
Kawahara equation [13,14]. The fKdV equation plays an
important role in many applications in plasma physics and
in hydrodynamics. For a detailed derivation of the fKdV
equation in a hydrodynamical context see [15]. A crucial
point of interest of the fKdV equation from our point of
view is that the familiar solitary wave solutions of the KdV
equation are deformed into oscillon-type objects, losing
continuously some of their mass by radiating small
amplitude waves in the direction of propagation [16]. It
has been proven that the spatially localized solitary
traveling wave solution of the KdV equation ceases to
exist when a fifth order dispersion term proportional to $\epsilon^2$ is
added [17,18]. What happens is that the extra dispersion
term causes the KdV solitary wave to develop a radiating
tail whereby it loses energy [16]. Bounded, stationary
solutions of the fKdV equation are weakly localized, in
that asymptotically they tend to a standing wave “tail.”
This is to be contrasted to the exponential falloff of the
well localized KdV solitons. In this context such weakly
localized, stationary solutions, baptized “nanopterons”
have been studied systematically in Refs. [19,20].

For both the quasibreathers and the nonlocal solitons of
the fKdV equation, the tail amplitude is exponentially small
in terms of a parameter ($\epsilon$) which characterizes the specific

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI. Funded by SCOAP³.
perturbation or describes the core amplitude (for quasi-breathers), and to leading order it can be written as

$$ \alpha_m = \frac{\lambda}{\epsilon^\nu} \exp\left(-\frac{\sigma}{\epsilon}\right), \quad (1) $$

where $\lambda$, $\nu$, and $\sigma$ are constants. For the fKdV case the exponent in the denominator is $\nu = 2$, $\sigma = \pi/2$ (in suitable units) while for spherically symmetric oscillons $\nu = (d-1)/2$ where $d$ is the number of spatial dimensions. The proportionality constant $\lambda$ is rather nontrivial to compute. Following the technique pioneered in Ref. [10], the value of $\lambda$ can be obtained by going through a complicated asymptotic matching calculation in the complex plane. For the fKdV equation the leading order result has been obtained in Ref. [21] with the result $\lambda \approx 19.969\pi$.

It is desirable to compute higher order corrections to the leading order result (1) from both a theoretical and a practical point of view. One would like to establish a systematic framework for such calculations, quantify the contributions of higher order terms, but it is also necessary to know higher order corrections in order to compare the results to numerical simulations for small, but obviously finite values of $\epsilon$. For very small $\epsilon$ values where the leading order term is supposed to dominate, the tail amplitude becomes so small that it is extremely demanding and difficult to calculate by numerical methods. This emphasizes the necessity to investigate corrections to Eq. (1). For the fKdV problem and also for quasibreathers, it is to be expected that the leading order result (1) gets several types of corrections, with the dominant ones involving only $\epsilon$. In the simplest case such corrections can be written as a power series of $\epsilon$ of the form

$$ \alpha_m = \frac{\lambda}{\epsilon^\nu} \exp\left(-\frac{\sigma}{\epsilon}\right) (1 + \zeta_1 \epsilon + \zeta_2 \epsilon^2 + \zeta_3 \epsilon^3 + \dots), \quad (2) $$

where $\zeta_j$ are some yet unknown constants. (For $d \ge 2$ dimensional oscillons there are also $\epsilon \ln \epsilon$ and similar logarithmic terms [9].) For most cases including the scalar field problem the tail amplitude has been calculated only to leading order; i.e., even $\zeta_1$ is unknown yet. On the other hand, for the fKdV equation detailed properties of weakly nonlocal solitons have been under intense scrutiny, and even the computation of $\zeta_1$ and $\zeta_2$ was undertaken by Grimshaw and Joshi in Ref. [22]. Since we expect that similar methods to that of Ref. [22] will be applicable for the more complicated scalar field problem as well, a thorough understanding of the higher order contributions in the fKdV problem appears to be necessary to us, which is the subject of the present paper.

According to the results of Ref. [22], the first two corrections in Eq. (2) are given as $\zeta_1 = -\pi$ and $\zeta_2 = \pi^2/2 \approx 4.935$ in suitable units. The value of $\zeta_1$ is in agreement with the numerical computations by

Boyd [20,23,24], using a multiple precision pseudospectral code. According to the numerical estimates of Boyd, however, $-0.1 < \zeta_2 < 0$, which is clearly in disagreement with the result of Ref. [22]. As far as we know, the reason for this discrepancy has remained unknown up to now.

We have succeeded to compute higher order corrections in Eq. (2) up to order $\epsilon^5$, and while confirming the result for the value of $\zeta_1 = -\pi$ of Ref. [22], we have found $\zeta_2 = \pi^2/2 - 5 \approx -0.065198$, a value which is quite consistent with the numerical estimates of Ref. [24].

To check our analytical results we have also developed a multiple precision pseudospectral numerical code to solve the fKdV equation. Since the tail amplitude for the intended $\epsilon$ parameter range can be many orders of magnitude smaller than $10^{-16}$, the use of arbitrary precision arithmetic is indispensable. The speed of present day personal computers and the efficiency of the available multiple precision arithmetic libraries allows us to reach significantly higher precision than what was available earlier. The running time of our code also allows us to apply a numerical minimization algorithm to find the phase of the tail where the tail amplitude is really minimal. This was not included in earlier publications. Comparison of the coefficients of the higher order corrections, $\zeta_1 \cdots \zeta_5$, with the results of our numerical simulations shows remarkably good agreement. In particular, the numerical value of $\zeta_2$ agrees with our perturbative result to five digits of precision, leaving little doubt as to its correctness. We have also obtained a very high order ($\sim 100$) asymptotic expansion in $\epsilon$ of the phase of the minimal amplitude wave tail in the Wentzel-Kramers-Brillouin (WKB) approximation. We could compare up to order 15 the expansion of the phase with the numerical results and an excellent agreement has been found.

The plan of the paper is the following: In Sec. II the fifth order KdV equation is introduced and its basic properties are discussed. Next, in Sec. III our implementation of the spectral method is described in some detail. In Sec. IV the small $\epsilon$ expansion for the “core” part of the solution is carried out. In Sec. V in the framework of the WKB approximation, the linearized solution around the core is determined to arbitrary order. The amplitude of the wave tail is given up to order 6, and the phase is given explicitly up to order 11. Sections VI and VII contain the most important calculations, namely carrying out the asymptotic matching to fourth order in $\epsilon$ of the complex extension of the “inner” and of the “outer” part of the solution, near the first singularity in the complex plane. Section VIII contains our conclusions.

## II. FIFTH ORDER KdV EQUATION

The Korteweg-de Vries equation modified by a small fifth derivative term can be written as [15,21]

$$ \epsilon^2 u_{yyyyy} + u_{yyy} + 6uu_y + u_t = 0, \quad (3) $$
where $u$ is a function of the time $t$ and spatial coordinate $y$, and $\epsilon$ is a small non-negative parameter. The indices denote derivatives with respect to $t$ and $y$. We only consider stationary solutions traveling with speed $c$ to the right, so that $u$ is time independent when $x = y - ct$ is used as a comoving spatial coordinate,

$$ \epsilon^2 u_{xxxx} + u_{xxx} + (6u - c)u_x = 0. \quad (4) $$

This form clearly shows that if $u$ is a solution and $u_c$ is a constant, then $u + u_c$ is also a solution moving with speed $c + 6u_c$. The equation can be integrated once, yielding

$$ \epsilon^2 u_{xxxx} + u_{xx} + 3u^2 - cu = M, \quad (5) $$

where $M$ is a constant that can be interpreted as the mass flux [22]. If $u$ and its derivatives tend to zero in either the positive or negative $x$ direction, then necessarily $M = 0$. For $\epsilon > 0$ bounded solutions necessarily have an asymptotic oscillating tail in both directions. In this case $M$ becomes determined by the boundary conditions. Two different boundary conditions have been used in the literature. One possibility, used by Boyd [20,23,24], is to impose no mass flux by requiring $M = 0$. Another choice, applied by Grimshaw and Joshi [22], is to set the integral of $u$ zero for an interval of $x$ corresponding to the wavelength of one oscillation, thereby keeping the fluid volume constant with respect to the $u = 0$ solution. Even in this latter case $M$ is extremely small, of order $\alpha^2$, where $\alpha$ is the tail amplitude, which turns out to be exponentially small in $\epsilon$. An $\alpha^2$ order small shift in $u$ and $c$ can be used to set $M$ zero. In the present paper we are interested in solutions for relatively small values of $\epsilon$. From now on we set $M = 0$, and our aim is to solve the equation

$$ \epsilon^2 u_{xxxx} + u_{xx} + 3u^2 - cu = 0, \quad (6) $$

for positive values of the two parameters $\epsilon$ and $c$. Equation (6) remains invariant under the rescalings

$$ u = \xi^2 \bar{u}, \quad x = \frac{1}{\xi} \bar{x}, \quad c = \xi^2 \bar{c}, \quad \epsilon = \frac{1}{\xi} \bar{c}, \quad (7) $$

for any $\xi > 0$ constant. Note that $\epsilon^2 c$ remains invariant. It would be possible to use this freedom to scale either of the constants $\epsilon$ or $c$ to some given value, and consider (6) as an equation containing only the other parameter. Following Grimshaw and Joshi [22], we use this rescaling freedom to set the speed parameter $c$ to a specific known function of $\epsilon$. The actual form of the function $c(\epsilon)$ will be fixed a bit later by the requirement that the spatial decay rate of the core of the solution should be parameter independent. Still, with this choice we reduce the number of parameters in the problem from two to only one, keeping only $\epsilon$.

We denote the value of $c \equiv c(\epsilon)$ in the $\epsilon \to 0$ limit by $c_0$. For $\epsilon = 0$ and $c_0 > 0$ Eq. (6) has the well-known KdV solitary wave solution

$$ u_0 = 2\gamma^2 \operatorname{sech}^2(\gamma x), \quad (8) $$

where the positive $\epsilon$ independent constant $\gamma$ is defined by

$$ c_0 = 4\gamma^2. \quad (9) $$

Localized asymptotically decaying solutions of the KdV equation only exist for $c_0 > 0$. Since our aim is to look for solutions that are as similar to solitary waves as possible, we assume that both $c_0$ and $c(\epsilon)$ are positive.

We are looking for solutions for which asymptotically $u$ becomes small. The linearization of Eq. (6) around $u = 0$ has four independent solutions. The exponentially decaying or growing solutions have the form $u = \exp(2\tilde{\gamma}x)$, where

$$ \epsilon^2(4\tilde{\gamma}^2)^2 + 4\tilde{\gamma}^2 - c = 0. \quad (10) $$

Since we assume $c > 0$, this has two real valued solutions for $\tilde{\gamma}$, determined by

$$ 4\tilde{\gamma}^2 = \frac{-1 + \sqrt{1 + 4\epsilon^2 c}}{2\epsilon^2}. \quad (11) $$

The remaining two oscillating solutions can be written in the form $u = \exp(ikx/\epsilon)$. This is a natural parametrization, since the spatial frequency tends to infinity when $\epsilon$ goes to zero. The linearization of (6) gives two real solutions for $k$, given by

$$ \frac{k^2}{\epsilon^2} = \frac{1 + \sqrt{1 + 4\epsilon^2 c}}{2\epsilon^2}. \quad (12) $$

Subtracting (11) from (12) and rearranging, we obtain the simple relation between $k$ and $\tilde{\gamma}$,

$$ k^2 = 1 + 4\epsilon^2 \tilde{\gamma}^2. \quad (13) $$

Taking the $\epsilon \to 0$ limit of Eq. (10) and comparing with (9), it follows that $\tilde{\gamma}$ tend to $\gamma$ when $\epsilon$ goes to zero. Grimshaw and Joshi [22] made the choice to impose that $\tilde{\gamma}$ is $\epsilon$ independent, in which case obviously $\tilde{\gamma} = \gamma$. This is a natural option, since it makes the decay rate of the core independent of the parameter $\epsilon$, and the expansion formalism becomes considerably simpler. We will also assume $\tilde{\gamma} = \gamma$ from now on. As a consequence of this choice, the parameter $c$ necessarily becomes $\epsilon$ dependent, and from (10) we obtain

$$ c = 4\gamma^2 + 16\gamma^4\epsilon^2. \quad (14) $$
This is the specific $c \equiv c(\epsilon)$ function that we ensure by the appropriate use of the rescaling in (7).

Since we are looking for real solutions, instead of the exponential form $u = \exp(ikx/\epsilon)$ for the tail, we write

$$u = \alpha_{\pm} \sin \left( \frac{k|x|}{\epsilon} - \delta_{\pm} \right). \qquad (15)$$

The amplitude $\alpha_+$ and phase $\delta_+$ in the positive direction can be different from that of $\alpha_-$ and $\delta_-$ in the negative direction. Using the trigonometric form it is enough to use the positive root in (13); hence in the following we will set

$$k = \sqrt{1 + 4\gamma^2\epsilon^2}. \qquad (16)$$

A conserved quantity defined as

$$F = -\frac{1}{2} cu^2 + 2u^3 + uu_{xx} - \frac{1}{2} u_x^2 + \epsilon^2 \left( uu_{xxxx} - u_x u_{xxx} + \frac{1}{2} u_{xx}^2 \right) \qquad (17)$$

can be interpreted as an energy flux [16,22]. Using (4) it is easy to check that $F_x = 0$. Substituting the form (15) of the tail, to leading order the energy flux turns out to be $F = \alpha_{\pm}^2/(2\epsilon^2)$. Since $F$ is constant, it follows that the amplitude of the tail at the two directions must necessarily agree, at least for small amplitudes. The phases can also be made to agree by a small shift in $x$, but the argument is not giving any information about the symmetry of the core. Numerical simulations also support the conjecture that solutions which have a single large core and small tails in both the positive and the negative directions are necessarily reflection symmetric. The important consequence is that there are no solutions for which there is a small tail in one direction and exponential decay without any tail in the other direction. Our numerical simulations show that for each small $\epsilon$ there is a solution which decays exponentially to zero for $x > 0$, and it has a core region similar to the KdV solitary wave, but continuing further to the negative direction the solution blows up at some finite $x < 0$ before a standing wave tail could appear. Hence in the following we consider only solutions which are symmetric with respect to $x = 0$. For large $x > 0$ their tail is characterized by the parameters $\alpha$ and $\delta$,

$$u = \alpha \sin \left( \frac{kx}{\epsilon} - \delta \right). \qquad (18)$$

Even after using the rescaling (7) to ensure the intended $\epsilon$ dependence of $c$ according to (14), there is still a remaining freedom to set $c_0$ to any desired positive value. It is sufficient to solve Eq. (6) numerically or analytically for only one special choice of $c_0$, since other solutions can be obtained by

the above rescalings. The choice made by Boyd in [20,24] is $c_0 = 4$, which is a natural option since it corresponds to $\gamma = 1$. For easier comparison with Boyd's results we will also use $c_0 = 4$ in our numerical calculations.

In this paper we mostly use the notations introduced by Grimshaw and Joshi in [22], but we discuss the connection with the variables and equations used by Boyd [20,23,24]. Applying the rescaling (7) with $\xi = 1/\epsilon$ we get $\bar{\epsilon} = 1$, and hence the $\epsilon^2$ factor disappears in front of the fourth derivative term. Introducing a further rescaled function by $\bar{u} = \bar{v}/6$, the coefficient of the quadratic term can change, and we obtain the form of the fKdV equation used by Boyd,

$$\bar{v}_{\bar{x}\bar{x}\bar{x}\bar{x}} + \bar{v}_{\bar{x}\bar{x}} + \frac{1}{2}\bar{v}^2 - \bar{c}\bar{v} = 0. \qquad (19)$$

In this case the only parameter is $\bar{c} = \epsilon^2 c = 4\gamma^2\epsilon^2 + 16\gamma^4\epsilon^4$.

### III. NUMERICAL METHOD

We apply a pseudospectral numerical method (see, e.g., [25]) to solve Eq. (6), looking for solutions $u$ that are symmetric at $x = 0$. Since the tail has infinitely many oscillations, spatial compactification with standard Chebyshev expansion cannot be used efficiently in this case. Boyd [20,23,24] used an additional basis function to represent the oscillating tail. In our numerical simulations we match the solution to the tail given in (18) at the outer numerical boundary $x = L$. By appropriate rescalings it is always possible to arrange that $c$ is given by (14) with $\gamma = 1$. For any chosen parameter $\epsilon$ and phase $\delta$ the numerical problem can be solved to obtain a corresponding tail amplitude $\alpha$. We use two boundary conditions at $x = L$,

$$\epsilon^2 u_{xx} + k^2 u = 0, \qquad (20)$$

$$\epsilon \sin \left( \frac{kx}{\epsilon} - \delta \right) u_x - k \cos \left( \frac{kx}{\epsilon} - \delta \right) u = 0. \qquad (21)$$

Together with the symmetry assumption these conditions make the solution unique.

We introduce an alternative independent variable $\theta$ by $x = L \cos \theta$. The center $x = 0$ corresponds to $\theta = \pi/2$, and the outer boundary to $\theta = 0$. We fix some order $N$, and represent the function $u$ by $N$ Fourier components $U_n$,

$$u = \sum_{n=0}^{N-1} \frac{1}{\rho_n} U_n \cos(2n\theta), \qquad (22)$$

where

$$\rho_n = \begin{cases} 1 & \text{if } 1 \le n \le N-2, \\ 2 & \text{if } n = 0 \text{ or } n = N-1. \end{cases} \qquad (23)$$
We only include even Fourier components because of the reflection symmetry at the center. Actually, this corresponds to expansion in even indexed Chebyshev polynomials $T_{2n}(x/L)$, since $\cos(n\theta) = T_n(\cos\theta)$. The solution can be alternatively represented by its values at N collocation points $x_n = L \cos\theta_n$, where $\theta_n = \pi n / (2(N-1))$ for $0 \le n \le N-1$. For the function values we introduce the notation $\tilde{u}_n = u(x_n)$. The value of u at the center is $\tilde{u}_{N-1}$, while at the outer boundary it is $\tilde{u}_0$. Using (22), the collocation values can be obtained by matrix multiplication,

$$ \tilde{u}_n = \sum_{j=0}^{N-1} C_{nj} U_j, \quad C_{nj} = \frac{1}{\rho_j} \cos \frac{nj\pi}{N-1}. \quad (24) $$

The inverse transformation is

$$ U_n = \sum_{j=0}^{N-1} F_{nj} \tilde{u}_j, \quad F_{nj} = \frac{2}{N-1} C_{nj}. \quad (25) $$

This corresponds to calculating the Fourier coefficients by numerical integration based on the collocation points. In both directions, the transformation corresponds to type I discrete cosine transform, DCT-I.

Multiplication of functions can easily be calculated using the collocation values, while derivatives can be naturally obtained using Fourier coefficients. The Fourier coefficients $(U_{xx})_n$ of the second derivative of the function u represented by $U_n$ can be calculated as

$$ (U_{xx})_n = \sum_{j=0}^{N-1} D_{nj}^{(2)} U_j, \quad (26) $$

where

$$ D_{nj}^{(2)} = \begin{cases} \frac{8}{\rho_j L^2} j(j^2 - n^2) & \text{if } j \ge n + 1, \\ 0 & \text{if } j < n + 1. \end{cases} \quad (27) $$

The second derivative matrix in the collocation picture can be obtained as $\tilde{D}_{nj}^{(2)} = C_{nl} D_{lp}^{(2)} F_{pj}$, where the summation for repeated indices is understood. The fourth derivative matrix can be obtained most easily by multiplying the second derivative matrix by itself, $D_{nj}^{(4)} = D_{nl}^{(2)} D_{lj}^{(2)}$. For the boundary condition we also need the first derivative, but only at the boundary. This can be calculated by the scalar product $V_n \tilde{u}_n$, where

$$ V_n = \sum_{j=0}^{N-1} \frac{4j^2}{L} F_{jn}. \quad (28) $$

Since Eq. (6) is nonlinear, we use an iterative procedure, called the Newton-Kantorovich method to solve it (see, e.g., Appendix C of [25]). Suppose that at the nth step we

have an approximate solution $u^{(n)}$. The next approximation shall be $u^{(n+1)} = u^{(n)} + \Delta$. Substituting into (6) and linearizing for $\Delta$ we obtain the equation

$$ \epsilon^2 \Delta_{xxxx} + \Delta_{xx} + 6u^{(n)}\Delta - c\Delta = R, \quad (29) $$

where the residual is

$$ R = -\epsilon^2 u_{xxxx}^{(n)} - u_{xx}^{(n)} - 3(u^{(n)})^2 + c u^{(n)}. \quad (30) $$

According to our experience, less than ten steps of iteration is enough to get extremely high precision solutions for this problem. The iteration can be started from the KdV solitary wave solution (8). We consider (29) as N algebraic equations at the collocation points. The left-hand side can be considered as a linear matrix operator $L_{nj}$ multiplying the collocation values $\Delta_j$ of the unknown function $\Delta$. The collocation values $R_n$ of residuals can be calculated at each step from the previous approximation. We replace two lines of $L_{nj}$ by values enforcing the boundary conditions, while replacing the corresponding elements in $R_n$ by the previous error in the boundary conditions. The result turns out to be quite insensitive of which lines we replace; for example, we can use line 0 and $N-2$. After this, the matrix equation can be solved for $\Delta_j$. Updating $u^{(n)}$ by adding the calculated $\Delta$ yields the next approximation.

We have written equivalent C and C++ codes to solve the numerical problem. Since we are interested in comparing the numerical results to the analytical ones, we intend to calculate the tail amplitude $\alpha$ for relatively small $\epsilon$ values. It can be seen easily that the usual 16 or 19 digits arithmetic is not enough for our aims. The core amplitude is always close to 2. If the tail amplitude $\alpha$ is of the order $10^{-a}$, and we intend to calculate it to b digits of precision, then the whole numerical procedure should be carried out with at least $a+b$ digits of precision. There are freely available numerical packages for calculations with arbitrarily many digits of precision that are fast enough for our purposes. For our C++ codes we have used the Class Library for Numbers (CLN) [26]. Our codes using the C library for arbitrary-precision ball arithmetic (ARB) [27,28] turn out to be about 20 times faster, due to more advanced matrix manipulation methods. However, writing programs using ARB is more difficult, since it requires separate lines of codes for each elementary algebraic manipulation, such as addition or multiplication of numbers.

To illustrate the structure and precision of the obtained solutions we present some results for $\epsilon = 0.05$. We use $c_0 = 4$ (and consequently $\gamma = 1$) in all our numerical work. For the phase we choose $\delta = 6\gamma\epsilon = 0.3$ for this example, which corresponds to the linear approximation of the phase belonging to the minimal tail amplitude [see (65) and (57) for details on that]. To include a large enough portion of the tail we set the outer boundary at $L = 30$. In Fig. 1 we plot
<figure><img src="image_4.png" /><figcaption>FIG. 1. The red curve shows $u$ as a function of $x$. Subtracting the function $u_\alpha = \alpha \sin(kx/\epsilon - \delta)$ corresponding to the matched tail we obtain the blue curve.</figcaption></figure>

logarithmically the $x$ dependence of the function $u$. The downward spikes correspond to zero crossings in the tail. To illustrate the precision of the numerical solution and the matching, we also plot the difference of $u$ and the inward continuation of the matched tail (18), using the tail amplitude $\alpha = 4.811363414 \times 10^{-10}$ provided by the numerical code. The extremely good agreement of $u$ and the tail for $x > 23$ indicates that the obtained $\alpha$ is precise to 10 digits. To reach this precision we need at least $N = 500$ collocation points and 26 digits of precision during the whole numerical calculation. At this resolution the running time on our desktop computer is about 90 s for the CLN code, and only 2 s for the ARB code. Even if we increase the resolution further, the difference from the matched tail will not go below $10^{-20}$ for the present $\epsilon$ value. The reason for this is that we represent the matched tail with the linearized solution (18). If the tail amplitude is $\alpha$, then we make an error of order $\alpha^2$ by this choice.

Numerical calculation of symmetric solutions with very small tails have been first reported by Boyd in [23], applying spectral methods. In that paper the solution is matched to a higher order nonlinear representation of the tail. Boyd called the method cnoidal matching, since in the $\epsilon = 0$ case the spatially periodic higher amplitude KdV solutions are given by the elliptic cosine function cn. In our present paper we are interested in solutions with such tiny tails that the linear tail approximation is adequately precise. In a subsequent paper [24] Boyd has presented high precision results for the case when the asymptotic phase of the tail is $\delta = 0$. The use of multiple precision arithmetic allowed the calculation of extremely small tail amplitudes for $\epsilon < 0.05$. We have checked that our code reproduces the tail amplitudes listed in Table II of [24]. Note that since Boyd uses the variables discussed in Eq. (19), the amplitude $\alpha$ given in that paper is equal to $6\epsilon^2\alpha$ using our notations.

The speed of our code allows us to search numerically for the phase $\delta_m$ for which the tail amplitude is minimal,

<table><caption>TABLE I. Numerically calculated values of the minimal tail amplitude $\alpha_m$, the corresponding phase $\delta_m$, and the core radius $x_{\text{core}}$. The order of optimal truncation for the expansion describing the core is listed in the last column (see Sec. IV).</caption><thead><tr><th>ε</th><th>α<sub>m</sub></th><th>δ<sub>m</sub></th><th>x<sub>core</sub></th><th>N<sub>opt</sub></th></tr></thead><tbody><tr><td>0.15</td><td>4.1 × 10<sup>-2</sup></td><td>0.958</td><td>2.56</td><td>3</td></tr><tr><td>0.1</td><td>6.572 × 10<sup>-4</sup></td><td>0.60552</td><td>4.77</td><td>6</td></tr><tr><td>0.07</td><td>1.802403 × 10<sup>-6</sup></td><td>0.4207362</td><td>7.72</td><td>10</td></tr><tr><td>0.05</td><td>4.811363375 × 10<sup>-10</sup></td><td>0.3001268310</td><td>12.1</td><td>14</td></tr><tr><td>0.035</td><td>1.472008979 × 10<sup>-15</sup></td><td>0.2100205651</td><td>18.2</td><td>21</td></tr><tr><td>0.025</td><td>4.771438977 × 10<sup>-23</sup></td><td>0.1500037632</td><td>26.8</td><td>30</td></tr><tr><td>0.017</td><td>1.527829748 × 10<sup>-35</sup></td><td>0.1020005424</td><td>41.2</td><td>45</td></tr><tr><td>0.012</td><td>5.935328843 × 10<sup>-52</sup></td><td>0.0720000947</td><td>60.1</td><td>64</td></tr></tbody></table>

$\alpha = \alpha_m$. We are not aware of such a study in the literature. Using Brent's minimization method [29], usually about 20 iterations are enough to get the necessary precision. In Table I we list the phase $\delta_m$ and minimal amplitude $\alpha_m$ for several choices of the parameter $\epsilon$. In the table we also include the radius of the core, which we estimate by the lowest value of $x$ where there is a zero crossing in the function $u$. For $x > x_{\text{core}}$ the oscillating tail dominates. For higher $\epsilon$ values we can reach fewer digits of precision because the matching to the linear approximation of the tail brings an error of order $\alpha_m^2$. In these cases there is a slight dependence on how the outer boundary $L$ is chosen. For $\epsilon < 0.05$ we present only 10 digits, even if we can reach higher precision with our spectral code. The reason for such precise calculations is that we intend to compare to the higher order analytical results presented in the next sections. The smaller $\epsilon$ is the more computational resources are necessary. To get the values for $\epsilon = 0.012$ we have used 5000 collocation points with 105 digits arithmetic, and the calculation using the ARB library took several hours on a desktop computer.

## IV. EXPANSION PROCEDURE FOR THE CORE

From this section on we concentrate on analytical methods and compare them to our numerical results. We construct an asymptotic expansion, in powers of the small parameter $\epsilon$, which can be used to describe the core region of an almost localized solitary wave solution. Since the amplitude of the oscillating tail in the faraway region is exponentially small in $\epsilon$, this expansion is not able to describe those oscillations. Into Eq. (6) we substitute the expansions

$$u = \sum_{n=0}^{\infty} u_n \epsilon^{2n}, \qquad (31)$$

$$c = \sum_{n=0}^{\infty} c_n \epsilon^{2n}, \qquad (32)$$
where $u_n$ are functions of $x$ and $c_n$ are numbers [20,22]. To obtain a meaningful finite result, the asymptotic series has to be truncated at some positive integer order. The error of the approximation is the smallest when the series is truncated at the optimal order $N_{opt}$, which is expected to increase proportionally with $1/e$. The error of the optimally truncated series is anticipated to be exponentially small in $e$, just as the oscillating tail that this expansion cannot describe. Substituting into Eq. (6), the vanishing of the $\epsilon$ independent part gives the KdV equation

$$u_{0,xx} + 3u_0^2 - c_0u_0 = 0, \quad (33)$$

which has the solution given by (8). For a given $c_0$ this is the unique localized single-core solution when symmetry with respect to $x = 0$ is assumed. For $n > 0$, the vanishing of the coefficient of $\epsilon^{2n}$ yields

$$u_{n-1,xxxx} + u_{n,xx} + \sum_{j=0}^{n} (3u_j - c_j)u_{n-j} = 0. \quad (34)$$

If we assume that the functions are known up to order $n-1$, then this equation can be considered as a linear inhomogeneous differential equation for determining $u_n$,

$$u_{n,xx} + 6u_0u_n - c_0u_n = R_n, \quad (35)$$

where

$$R_n = -u_{n-1,xxxx} - 3 \sum_{j=1}^{n-1} u_j u_{n-j} + \sum_{j=1}^{n} c_j u_{n-j}. \quad (36)$$

General solutions of the homogeneous problem with $R_n = 0$ can be obtained as linear combinations of two solutions. The first solution is the derivative of $u_0$, which is antisymmetric. The other solution blows up exponentially at infinity. It follows that the solution of the inhomogeneous problem (35) which is symmetric with respect to $x = 0$ and localized has to be unique.

Proceeding order by order in $n$, it turns out that the inhomogeneous source term can be written as a finite sum of powers of $\mathrm{sech}^2(\gamma x)$,

$$R_n = \sum_{j=1}^{n+2} R_{n,j} \gamma^{2n+4} \mathrm{sech}^2(\gamma x). \quad (37)$$

The functions $u_n$ can be expanded similarly,

$$u_n = \sum_{j=1}^{n+1} u_{n,j} \gamma^{2n+2} \mathrm{sech}^2(\gamma x). \quad (38)$$

The powers of $\gamma$ are included in order to make $R_{n,j}$ and $u_{n,j}$ rational numbers without $\gamma$ factors. Comparing with (8) follows that $u_{0,1} = 2$.

Using the identities for the derivatives of $\mathrm{sech}^j x$ (see, e.g., Appendix A of [20]), from (36) it follows that for $n \ge 1$,

$$\begin{aligned} R_{n,j} = & -16j^4 u_{n-1,j} + 8(j-1)(2j-1)(2j^2-2j+1)u_{n-1,j-1} \\ & - (2j-4)(2j-3)(2j-2)(2j-1)u_{n-1,j-2} \\ & - 3 \sum_{l=1}^{n-1} \sum_{m=1}^{l+1} u_{l,m} u_{n-l,j-m} + \sum_{l=1}^{n-j+1} \hat{c}_l u_{n-l,j} \end{aligned} \quad (39)$$

where $\hat{c}_l = \gamma^{-2l-2}c_l$. This expression for $R_{n,j}$ is valid only if we substitute zero for every occurrence of $u_{n,j}$ when $j < 1$ or $j > n + 1$. If the coefficients $u_{n,j}$ are known up to order $n-1$ in the first index, then (39) can be used to calculate the source term $R_n$ in (35). For $j \ge 2$ all $\mathrm{sech}^{2j}(\gamma x)$ terms in $R_n$ can be generated from appropriate sech terms in $u_n$. However, when setting the right-hand side of (35) to $\mathrm{sech}^2(\gamma x)$ the symmetric, asymptotically decaying solution can be written as

$$u_n = \frac{1}{2} \mathrm{sech}^2(\gamma x) (1 - \gamma x \tanh(\gamma x)). \quad (40)$$

This solution goes to zero at infinity as $x \exp(-2\gamma x)$ which is a slower decay than the $\exp(-2\gamma x)$ decay of the other terms. Because of this, we must avoid this source term by setting $R_{n,1} = 0$ for all $n$. Since $R_{n,1}$ contains $\hat{c}_n$ linearly, this can be achieved at any order by the appropriate choice of $\hat{c}_n$. Proceeding with the calculation it turns out that $\hat{c}_0 = 4$, $\hat{c}_1 = 16$, and $\hat{c}_j = 0$ for all $j \ge 2$, consistently with (14).

From (35) it follows that for $2 \le j \le n + 1$,

$$4(j^2 - 1)u_{n,j} + [12 - (2j-1)(2j-2)]u_{n,j-1} = R_{n,j}. \quad (41)$$

For $j = n + 2$ we obtain

$$[12 - (2n + 3)(2n + 2)]u_{n,n+1} = R_{n,n+2}. \quad (42)$$

The equation for $j = 1$ is simply $R_{n,1} = 0$. If all $R_{n,j}$ coefficients are already calculated at order $n$, then $u_{n,n+1}$ can be obtained from (42). After this, all $u_{n,j-1}$ can be calculated one by one in decreasing order in $j$ using (41). This algorithm can easily be implemented using any algebraic manipulation software. Results for the coefficients up to order $n = 4$ are given in Table II. A similar algorithm to calculate $u_{n,j}$ has been presented in Table 10.5 of Boyd's book [20]. Note that there is a typo in the pseudocode there: the summation in the phase speed contributions should start from $m = 1$, not from 0.

If the solution $u$ has already been calculated precisely by some numerical method, we can compare it to various orders of the above $\epsilon$ expansion, defining the error of the $N$th order analytic approximation by
<table><caption>TABLE II. First few values of $u_{n,j}$.</caption><thead><tr><th>n↓ j→</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>0</td><td>2</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>1</td><td>-20</td><td>30</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2</td><td>60</td><td>-930</td><td>930</td><td>...</td><td>...</td></tr><tr><td>3</td><td>-2472</td><td>21036</td><td>-66216</td><td>49662</td><td>...</td></tr><tr><td>4</td><td>-240780/7</td><td>-3177030/7</td><td>23319570/7</td><td>-48197250/7</td><td>28918350/7</td></tr></tbody></table>

$$ \Delta u_N = u - \sum_{n=0}^{N} u_n \epsilon^{2n}. \qquad (43) $$

Choosing $\epsilon = 0.05$ and the symmetric solution $u$ with the minimal tail, in Fig. 2 we plot the functions $\Delta u_N$ for three values of $N$ for which the error is the smallest. Clearly, for this $\epsilon$ value the optimal truncation is at $N_{opt} = 14$. Since all $u_n$ decay exponentially, for large $x$ the difference should agree with the oscillating tail. The figure shows that for the optimal truncation this holds even in most of the inner region, since now the core radius is $x_{core} = 12.1$. It turns out that close to the center the error is even smaller. We have obtained similar plots for the other $\epsilon$ values listed in Table I. The values of $N_{opt}$ for various $\epsilon$ values are listed in the last column of the table. It can be checked that $N_{opt}$ increases proportionally to $1/\epsilon$. The contribution of the $n$th term in the expansion (31), i.e., the function $u_n\epsilon^{2n}$, is generally the smallest for $n = N_{opt} + 1$, or in some cases for $n = N_{opt}$, as one can expect it for asymptotic series.

## V. WKB SOLUTION

We intend to linearize Eq. (6) around some solution $u$. This solution may be a numerically obtained symmetric solution with a small tail in both directions, or an asymmetric solution that tends to zero for positive $x$.

<figure><img src="image_5.png" /><figcaption>FIG. 2. Difference of the $N$th order approximation of the $\epsilon$ expansion (31) from the minimal tail symmetric solution $u$ for $\epsilon = 0.05$. The horizontal lines show the tail amplitude $\alpha_m = 4.811363375 \times 10^{-10}$.</figcaption></figure>

The important point is that the core region should be approximated well by the expansion (31). Substituting $u \to u + u_w$ into (6) and linearizing gives

$$ \epsilon^2 u_{wxxxx} + u_{wxx} + 6uu_w - cu_w = 0. \qquad (44) $$

We will apply this formalism for $u_w$ which turns out to be exponentially small in terms of $\epsilon$; hence the linear approximation is well justified. As we have seen, in the asymptotic region there are oscillations with spatial frequency proportional to $1/\epsilon$. Hence we use the WKB method to search for solutions of (44). The first step is to substitute $u_w = \beta_c \exp A$, where $A$ is a function of $x$, and $\beta_c$ is a complex constant. Then we expand $A$ in powers of $\epsilon$, starting with a $1/\epsilon$ term,

$$ A = \sum_{n=-1}^{\infty} A_n \epsilon^n, \qquad (45) $$

and solve the resulting equation order by order in $\epsilon$. Since all $A_k$ appear only in differentiated form, there will be additive complex scalar freedom in all of these functions. All these can be absorbed into the complex valued $\epsilon$ dependent factor $\beta_c$.

To leading $\epsilon^{-2}$ order we obtain $(A_{-1x})^2[(A_{-1x})^2 + 1] = 0$. Since we are looking for high frequency solutions, we are not interested in the $A_{-1x} = 0$ solution. We continue with the choice $A_{-1x} = -i$, since the solution obtained from $A_{-1x} = i$ turns out to be the complex conjugate to all orders in $\epsilon$. Proceeding order by order in $\epsilon$, at each order we obtain a condition determining $A_{nx}$. The first few functions are

$$ A_{-1} = -ix, \qquad (46) $$

$$ A_0 = 0, \qquad (47) $$

$$ A_1 = -2i\gamma^2 x + 6i\gamma \tanh(\gamma x), \qquad (48) $$

$$ A_2 = 15\gamma^2 \operatorname{sech}^2(\gamma x), \qquad (49) $$

$$ A_3 = 2i\gamma^4 x + 111i\gamma^3 \operatorname{sech}^2(\gamma x) \tanh(\gamma x), \qquad (50) $$

$$ A_4 = \frac{525}{2} \gamma^4 \operatorname{sech}^2(\gamma x) [3\operatorname{sech}^2(\gamma x) - 2], \qquad (51) $$

$$ A_5 = -4i\gamma^6 x + \frac{3}{5}i\gamma^5 [12267\operatorname{sech}^4(\gamma x) - 4089\operatorname{sech}^2(\gamma x) + 632] \tanh(\gamma x), \qquad (52) $$

$$ A_6 = \frac{3}{2}\gamma^6 \operatorname{sech}^2(\gamma x) [49317\operatorname{sech}^4(\gamma x) - 49317\operatorname{sech}^2(\gamma x) + 8050]. \qquad (53) $$

It is natural to choose the value of the additive constants in $A_{2n-1}$ to make the functions antisymmetric at $x = 0$. The
form of $A_{2n}$ can be made unique by requiring that the functions tend to zero at infinity.

The terms proportional to $x$ in the odd indexed $A_n$ functions can be absorbed into the $A_{-1}e^{-1} = -ix/\epsilon$ term if we replace it by $-ikx/\epsilon$, where $k$ is the $\epsilon$ dependent constant defined in (16). This can be done by setting

$$ -\frac{ix}{\epsilon} + \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} A_n \epsilon^n = -\frac{ikx}{\epsilon} + \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} i \tilde{A}_n \epsilon^n, \quad (54) $$

where for positive integer $n$ the functions $\tilde{A}_{2n-1}$ are defined by

$$ A_{2n-1} = i\tilde{A}_{2n-1} - ix\gamma^{2n} \frac{(-1)^{n+1}(2n)!}{(2n-1)(n!)^2}. \quad (55) $$

Inserting the $k$ factor into the linear term is natural, since the asymptotic spatial frequency is $k/\epsilon$, which is valid to all orders in $\epsilon$. All these new functions with odd indices have a finite limit at infinity, which we denote by

$$ \tilde{\delta}_{2n-1} = \lim_{x \to +\infty} \tilde{A}_{2n-1}, \quad (56) $$

since they will determine the asymptotic phase shift of the minimal tail configuration. The first few values are

$$ \tilde{\delta}_1 = 6\gamma, \quad (57) $$

$$ \tilde{\delta}_3 = 0, \quad (58) $$

$$ \tilde{\delta}_5 = \frac{1896}{5}\gamma^5, \quad (59) $$

$$ \tilde{\delta}_7 = \frac{67140}{7}\gamma^7, \quad (60) $$

$$ \tilde{\delta}_9 = \frac{2662320}{7}\gamma^9, \quad (61) $$

$$ \tilde{\delta}_{11} = \frac{1301363652}{77}\gamma^{11}. \quad (62) $$

We have calculated $\tilde{\delta}_n$ for $n < 100$ using algebraic manipulation software.

For real $x$ all even indexed $A_n$ are real, and hence they will contribute to the amplitude of the linearized solution. The odd $A_n$ are all purely imaginary, so they will determine the phase. The general solution of the linearized problem that takes real values for real $x$ can be obtained by a linear combination of the solutions belonging to $A_{-1} = -ix$ and $A_{-1} = ix$,

$$ u_w = \beta \exp \left( \sum_{\substack{n=2 \\ \text{even}}}^{\infty} A_n \epsilon^n \right) \sin \left( \frac{kx}{\epsilon} - \delta_w - \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{A}_n \epsilon^n \right). \quad (63) $$

Here $\beta$ and $\delta_w$ are real constants with arbitrary $\epsilon$ dependence. They are related to the magnitude and phase of the complex constant $\beta_c$. The exponential term can be directly expanded in powers of $\epsilon$, providing

$$ u_w = \beta \left( 1 + \sum_{\substack{n=2 \\ \text{even}}}^{\infty} \tilde{A}_n \epsilon^n \right) \sin \left( \frac{kx}{\epsilon} - \delta_w - \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{A}_n \epsilon^n \right), \quad (64) $$

where the even indexed coefficients $\tilde{A}_n$ can easily be obtained from the original $A_n$. The first few functions are $\tilde{A}_2 = A_2$, $\tilde{A}_4 = A_4 + \frac{1}{2}A_2^2$, and $\tilde{A}_6 = A_6 + A_4A_2 + \frac{1}{6}A_2^3$. For $\delta_w = 0$ the function $u_w$ is antisymmetric with respect to $x = 0$, while for $\delta_w = \pi/2$ it is symmetric.

### A. Phase of the tail

Since the $A_{2n}$ functions tend to zero at infinity, $\beta$ gives the asymptotic amplitude of the oscillation represented by $u_w$ in (64). The amplitude in the core region is modified by the factor which is $O(1)$, so the linear correction $u_w$ remains small even in the core region. The asymptotic behavior of the function for positive $x$ is $u_w = \beta \sin(kx/\epsilon - \delta_w - \delta_m)$ with

$$ \delta_m = \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{\delta}_n \epsilon^n, \quad (65) $$

where the constants $\tilde{\delta}_n$ for odd $n$ are defined in (56). The same function near the center behaves as $u_w \sim \sin(kx/\epsilon - \delta_w)$, since $\tilde{A}_{2n-1} = 0$ at $x = 0$ in (64). The asymptotic phase comes from two contributions. The part $\delta_w$ gives the phase near the center, while the $\epsilon$ dependent constant $\delta_m$ gives the additional phase shift between the center and positive infinity.

Another way to consider $u_w$ in (64) is to decompose it to a sine part with $\delta_w = 0$ and a cosine part corresponding to $\delta_w = \pi/2$, both with arbitrary amplitudes,

$$ u_w = \left( 1 + \sum_{\substack{n=2 \\ \text{even}}}^{\infty} \tilde{A}_n \epsilon^n \right) \left[ \beta \sin \left( \frac{kx}{\epsilon} - \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{A}_n \epsilon^n \right) + \beta \cos \cos \left( \frac{kx}{\epsilon} - \sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{A}_n \epsilon^n \right) \right]. \quad (66) $$

The sine part is antisymmetric at the center $x = 0$, while the cosine part is symmetric. The cosine part with arbitrary $\beta_{\cos}$ amplitude can be added to any symmetric solution $u$ of (6), showing that the symmetric solution is not unique. On the other hand, the sine part of the tail of any symmetric solution has a fixed amplitude. Let us suppose that a certain symmetric solution $u$ of (6) has a small-amplitude tail given by (18) at large distances for $x > 0$. This tail can also be written as
$$u = \alpha \cos(\delta - \delta_m) \sin\left(\frac{kx}{\epsilon} - \delta_m\right) - \alpha \sin(\delta - \delta_m) \cos\left(\frac{kx}{\epsilon} - \delta_m\right). \quad (67)$$

Adding $u_w$ with $\beta_{\sin} = 0$ and $\beta_{\cos} = \alpha \sin(\delta - \delta_m)$, the cosine part of the tail becomes completely canceled, and we obtain the minimal tail symmetric solution, $u_m = u + u_w$. Adding $u_w$ with nonzero $\beta_{\sin}$ would destroy the central symmetry of $u$. For given $\epsilon$ the minimal tail solution is unique, and the tail for $x > 0$ has the asymptotic behavior

$$u_m = \alpha_m \sin \left( \frac{kx}{\epsilon} - \delta_m \right), \quad (68)$$

where $\alpha_m = \alpha \cos(\delta - \delta_m)$. Consequently, the minimal tail-amplitude symmetric solution necessarily has the phase $\delta_m$ in the tail, which has already been calculated in terms of the asymptotic series (65).

In Fig. 3 we compare the numerically calculated minimal-amplitude asymptotic phase $\delta_m$ and its various order approximations provided by (65),

$$\delta_m^{(j)} = \sum_{\substack{n=1 \\ \text{odd}}}^j \tilde{\delta}_n \epsilon^n. \quad (69)$$

We plot logarithmically the relative difference $\Delta\delta_j = |(\delta_m - \delta_m^{(j)})/\delta_m|$ for $j \le 15$. Since $\tilde{\delta}_3 = 0$, naturally $\Delta\delta_1 = \Delta\delta_3$. For $\epsilon \lesssim 0.07$ the numerical values of $\Delta\delta_j$ are decreasing as $\epsilon^{j+1}$. The deviations from the straight line of $\tilde{\delta}_{13}$ and $\tilde{\delta}_{15}$ for $\epsilon = 0.012$ show that the numerical result for this $\epsilon$ value is less than 18 digits precise. This is due to the necessary high resolution and long running times for such small values of $\epsilon$. Although we do not show $\Delta\delta_j$

<figure><img src="image_6.png" /><figcaption>FIG. 3. Log-log plot of $\Delta\delta_j = |(\delta_m - \delta_m^{(j)})/\delta_m|$, showing the relative difference of the numerically calculated $\delta_m$ from its various order analytic approximations $\delta_m^{(j)}$.</figcaption></figure>

for higher $j$ in the figure, the precision of this agreement can be further improved for not so small $\epsilon$ values. For example, for $\epsilon = 0.02$ our numerically calculated $\delta_m$ agrees to 34 digits with $\delta_m^{(47)}$. The fact that we obtain numerical results correct to so many digits of precision shows the remarkable power of the exponentially convergent spectral method when it is combined with arbitrary precision arithmetic. On the other hand, for high $\epsilon$ values it is apparent that (65) is indeed an asymptotic series. For example, for $\epsilon = 0.15$ the best approximation of the numerical result is given by $\delta_m^{(9)}$, higher orders giving larger and larger errors. In this case the smallest $\tilde{\delta}_n \epsilon^n$ contribution to the sum in (65) belongs to $n = 11$, which is just a little smaller than the contribution of the $n = 9$ term. In general, an asymptotic series is expected to give the best approximation when the summing is stopped at the term that gives the smallest contribution to the result.

Another important consequence of (67) is that for a given $\epsilon$ the tail amplitude $\alpha$ of any symmetric solution with phase $\delta$ is related to the minimal tail amplitude $\alpha_m$ by

$$\alpha = \frac{\alpha_m}{\cos(\delta - \delta_m)}. \quad (70)$$

Since $\alpha_m$ is exponentially small, this relation is valid to higher polynomial orders in $\epsilon$. To numerically check Eq. (70), for various $\epsilon$ values we have calculated $\delta_m$, $\alpha_m$, and for different $\delta$ phase shifts the amplitudes $\alpha$. According to our results, the error in (70) turns out to be order $\alpha^2$, similar to the error in the numerically calculated $\alpha$ caused by the linear tail approximation. Because of the high precision of this relation it is appropriate to concentrate on the minimal tail amplitude in the following. We should like to point out that in the limit of $\delta \to \pi/2$, Eq. (70) yields

$$\alpha|_{\delta=\frac{\pi}{2}} = \alpha_m \left( \frac{1}{6\gamma\epsilon} + \mathcal{O}(\epsilon) \right), \quad (71)$$

which agrees with the result of Ref. [30]. Taking $u_m$ as the minimal tail symmetric solution with asymptotics in (68) and subtracting $u_w$ with $\beta = \alpha_m$ and $\delta_w = 0$ in (64), the oscillating tail of the asymmetric solution $u_- = u_m - u_w$ becomes totally canceled in the positive $x$ direction. Using the linear approximation, the solution $u_-$ would appear to have double tail amplitude in the negative directions. However, it has been shown in [22] using an energy flux conservation law that no such solution can exist. From numerical simulations we can actually see that although the core domain remains quite similar to that of the symmetric one, the asymmetric solution blows up in the negative $x$ direction before the double-amplitude tail could appear. Nevertheless, $u_-$ plays an important role in the analytical calculation of the minimal tail amplitude.
# VI. ASYMPTOTIC MATCHING ON THE COMPLEX PLANE

## A. Complex extension

Since the amplitude $\alpha$ of the tail is exponentially small in terms of $\epsilon$, it cannot be determined by the direct expansion and WKB methods that we have used for the calculation of the phase $\delta$. To deal with the amplitude we have to extend the functions to the complex $x$ plane and study them near the singularity closest to the real $x$ axis, applying the method of matched asymptotic expansions [10,21]. We will apply the Laplace transform to solve the inner problem [22].

We extend analytically Eq. (6) and its solution $u$ to the complex $x$ plane. The value of the function $u$ at the point $x = x_r - ix_i$ has to be the complex conjugate of the value at $x = x_r + ix_i$, according to the Schwarz reflection principle. If the original real function is symmetric at $x = 0$, then the extension will naturally satisfy $u(x_r + ix_i) = u(-x_r - ix_i)$. In this case the value of $u$ at $x = -x_r + ix_i$ will be the complex conjugate of the value at $x = x_r + ix_i$, and the function must take real values on the imaginary axis. The extension of an antisymmetric function, satisfying $u(-x) = -u(x)$, has the property that $u(-x_r + ix_i)$ is -1 times the complex conjugate of $u(x_r + ix_i)$, and it must be purely imaginary on the imaginary axis.

We also extend to the complex $x$ plane the asymptotic expansion (31) of $u$, where $u_n$ depends on $x$ according to (38), with already calculated coefficients $u_{nj}$. All $x$ dependence is through powers of the function $\operatorname{sech}^2(\gamma x)$, which is singular at the points $x = (2n + 1)i\pi/(2\gamma)$, for integer $n$. The closest singularity above the real axis is at $x = i\pi/(2\gamma)$. In a neighborhood of this singularity we define the rescaled $q$ complex coordinate by

$$x = \frac{i\pi}{2\gamma} + \epsilon q. \qquad (72)$$

Near the singularity $\operatorname{sech}^2(\gamma x)$ has a Laurent series expansion starting with a $\epsilon^{-2}q^{-2}$ term. The expansion of $u$ can be obtained by adding various powers of this using some algebraic manipulation software. Since $u$ is growing as $\epsilon^{-2}$ near the singularity, it is natural to write the result in terms of the rescaled function

$$v = \epsilon^2 u. \qquad (73)$$

Substituting (72) into (31), expanding in powers of $1/q$, and then also in powers of $\epsilon$, we obtain that

$$v = \sum_{n=0}^{\infty} \gamma^{2n} \epsilon^{2n} v_n, \qquad (74)$$

where the expansions of the first three functions are

$$v_0 = -\frac{2}{q^2} + \frac{30}{q^4} - \frac{930}{q^6} + \frac{49662}{q^8} - \frac{28918350}{7q^{10}} + \cdots, \qquad (75)$$

$$v_1 = \frac{2}{3}, \qquad (76)$$

$$v_2 = -\frac{2q^2}{15} + \frac{2}{3} + \frac{64}{q^2} + \frac{5856}{5q^4} - \frac{827520}{7q^6} + \cdots. \qquad (77)$$

Since there are only even powers of $q$ with real coefficients, any truncated versions of the above series correspond to the complex extension of symmetric functions. Apart from the exactly known constant $v_1$ all $v_n$ functions are given in terms of asymptotic expansions in $1/q$. We use these expansions in some domain where both $\epsilon q$ and $1/q$ are small; hence $\epsilon^{2n}v_n$ can be small too.

## B. Inner problem

Using the rescaled function $v$ and the complex coordinate $q$ Eq. (6) becomes

$$v_{qqqq} + v_{qq} + 3v^2 - \epsilon^2 cv = 0, \qquad (78)$$

where $c$ also depends on $\epsilon$ according to (14). If we search for solutions of this equation as an expansion of the form (74), then taking the various $\epsilon^n$ contributions we obtain differential equations for $v_n$. The first two equations are

$$v_{0qqqq} + v_{0qq} + 3v_0^2 = 0, \qquad (79)$$

$$v_{1qqqq} + v_{1qq} + 2v_0(3v_1 - 2) = 0. \qquad (80)$$

The equation for $v_2$ will be studied in Sec. VII. Independently of $v_0$, the second equation can always be solved by $v_1 = 2/3$. Since this also agrees with the result obtained in (76), we will use this solution for $v_1$ from now on. These equations for $v_n$ are the $n$th order representations of the so-called inner problem. The expansion solutions (75)–(77) of the outer problem will be used as matching conditions for $v_n$. These asymptotic expansions provide valid boundary conditions when $Imq < 0$ is fixed and $Req \to +\infty$, determining unique inner solutions, which we denote by $v_n^{(-)}$. The functions $v_n^{(-)}$ do not have any oscillating tail in the positive direction $Req > 0$. They can be associated with the unique asymmetric solution $u_-$ of the original Eq. (6) for which there is no tail in the positive $x$ direction.

### 1. Zeroth order inner problem

An appropriate precise solution of the $\epsilon$ independent Eq. (79) can be used to determine the tail amplitude of the solution $u$ of (6) to leading order in $\epsilon$. An expansion solution for $v_0$ consistent with (75) can be searched for in the form
$$v_0 = \sum_{n=1}^{\infty} b_n^{(0)} q^{-2n}. \qquad (81)$$

Equation (79) gives the following equation for the coefficients:

$$(2n-2)(2n-1)(2n)(2n+1)b_{n-1}^{(0)} + (2n)(2n+1)b_n^{(0)} + 3 \sum_{j=1}^{n} b_j^{(0)} b_{n-j+1}^{(0)} = 0. \qquad (82)$$

For $n=1$ there is no $b_{n-1}^{(0)}$ term, and the nonzero solution is $b_1^{(0)} = -2$. For $n \ge 2$ we can obtain a recursion relation by taking the first and the last terms out of the summation,

$$(2n-3)(2n+4)b_n^{(0)} = -(2n-2)(2n-1)(2n)(2n+1)b_{n-1}^{(0)} - 3 \sum_{j=2}^{n-1} b_j^{(0)} b_{n-j+1}^{(0)}. \qquad (83)$$

The obtained coefficients are the same as those in (75). For large $n$ the coefficients diverge as $b_n^{(0)} \sim (-1)^n (2n-1)!$.

Following the method introduced in [22], we look for the solution of (79) in the form of a Laplace transform of a function $V'_0(s)$,

$$v_0 = \int_{\Gamma} I_0(s) ds, \quad I_0(s) = \exp(-sq)V'_0(s), \qquad (84)$$

where the contour $\Gamma$ is from $s=0$ to infinity, satisfying $\text{Re}(sq) > 0$. The prime is here because $V'_0(s)$ is the derivative of a function $V_0(s)$ that can be defined in terms of the Borel transformation [21]. However, we will only use the $V'_0(s)$ defined by (84) in the following. Using the identity for the Laplace transform of powers of $s$ we can see that $V'_0(s)$ can be expanded as

$$V'_0(s) = \sum_{n=0}^{\infty} a_n^{(0)} s^{2n+1}, \qquad (85)$$

where

$$a_n^{(0)} = \frac{b_{n+1}^{(0)}}{(2n+1)!}. \qquad (86)$$

This can also be written as $b_n^{(0)} = (2n-1)!a_{n-1}^{(0)}$. Substituting into (82) we get

$$a_{n-1}^{(0)} + a_n^{(0)} + \frac{3}{(2n+3)!} \sum_{j=0}^{n} (2j+1)!(2n-2j+1)!a_j^{(0)} a_{n-j}^{(0)} = 0. \qquad (87)$$

Separating the first and last terms in the summation and using that $a_0^{(0)} = -2$, one can obtain a recursion relation which is valid for $n \ge 1$,

$$\frac{(n+3)(2n-1)}{(n+1)(2n+3)} a_n^{(0)} = -a_{n-1}^{(0)} - \frac{3}{(2n+3)!} \sum_{j=1}^{n-1} (2j+1)!(2n-2j+1)! a_j^{(0)} a_{n-j}^{(0)}. \qquad (88)$$

The singularities of the function $V'_0(s)$ will be determined by the large $n$ behavior of the coefficients. The leading order behavior is $a_n^{(0)} \approx K(-1)^n$, where $K \approx 19.97$. The constant $K$ will be very important in the following, since it will determine the radiation amplitude. Hence we intend to calculate $K$ to several digits precision. Unfortunately, no fully analytical method is known for this calculation. We look for the large $n$ behavior of the coefficients in the form

$$a_n^{(0)} = (-1)^n \tilde{a}_n^{(0)}, \quad \tilde{a}_n^{(0)} = K G_0(n), \quad G_0(n) = 1 + \sum_{j=1}^{\infty} \frac{g_j}{n^j}. \qquad (89)$$

The constants $g_j$ should be determined by substituting this expansion into (87). We assume that $n$ is large, but if we are interested in a finite number of $g_j$ constants, we do not have to take into account all the $n+1$ terms in the summation in (87). We substitute (89) into the equation

$$\tilde{a}_n^{(0)} - \tilde{a}_{n-1}^{(0)} + \sum_{j=0}^{j_m} (-1)^j W_{n,j}^{(0)} \tilde{a}_{n-j}^{(0)} = 0, \qquad (90)$$

where

$$W_{n,j}^{(0)} = \frac{6}{(2n+3)!} (2j+1)!(2n-2j+1)! a_j^{(0)}, \qquad (91)$$

and $j_m$ is some positive integer.

The more terms we intend to determine for $G_0(n)$ in (89), the higher $j_m$ we should choose. However, the constants $W_{n,j}^{(0)}$ only involve $a_j^{(0)}$ with small $j$, so they can be calculated explicitly. We can use algebraic manipulation software to substitute a truncated version of the expansion (89) into (90). Taking the coefficients in increasing powers of $1/n$, we can determine the constants $g_j$. The first seven terms of the result yield

$$G_0(n) = 1 - \frac{3}{n} + \frac{39}{4n^2} - \frac{69}{2n^3} + \frac{1929}{16n^4} - \frac{3381}{8n^5} + \frac{46041}{32n^6} - \frac{1089483}{224n^7} + \dots \qquad (92)$$
To get the correct coefficients up to this order one has to set at least $j_m = 3$. This means that we use at least the first four and last four terms from the summation in (87), but can neglect the others in between. The importance of (92) is that it allows us to determine the constant $K$ to several digits precision by calculating the concrete coefficients from the recursion up to moderately high $n$ values. Using some algebraic manipulation program, the calculation of $a_n^{(0)}$ by (88) can be made faster by using floating point arithmetic valid to hundred digits precision instead of using exactly represented but very long rational numbers. In this way the first few thousand coefficients can be calculated in a couple of minutes. The approximation for the proportionality constant can be calculated using an appropriately truncated version of $G_0(n)$ as $K \approx a_n^{(0)}(-1)^n/G_0(n)$. The result up to 22 digits is

$$K = -19.96894735876096051827. \quad (93)$$

The high precision will be useful because we intend to study higher order $\epsilon$ corrections in the following sections.

The series (85) for $V'_0(s)$ is convergent for $|s| < 1$ and can be analytically extended for larger $|s|$, showing that $V'_0(s)$ is unique. On the other hand, several different $v_0$ functions can be obtained from it using the Laplace transform (84), depending on how the path $\Gamma$ is located with respect to the singularities of $V'_0(s)$. The function $V'_0(s)$ satisfies an integral equation presented in Eq. (31) of [22]. All singularities are located on the imaginary axis at $s = \pm ni$, where $n$ is any positive integer. There is no singularity at $s = 0$.

Calculating $v_0$ by (84) in a domain where $\arg(q)$ is close to zero, it is natural to choose the positive real $s$ axis as the curve $\Gamma$. The solution obtained in this way is the asymmetric $v_0^{(-)}$, which is determined by the boundary condition (75) for $\text{Req} > 0$. We can extend the $v_0^{(-)}$ defined by the Laplace transform integral to the domain $-\pi < \arg(q) \le 0$, but for $\arg(q) \le -\pi/2$ the contour $\Gamma$ cannot remain on the real axis, and it should move into the region where $0 < \arg(s) < \pi/2$. However, if the function $v_0^{(-)}$ is defined smoothly on the positive part of the real $q$ axis, it cannot be extended by the integral (84) to the negative part of the axis, where $\arg(q) = -\pi$. That would require crossing the singularities on the upper half of the imaginary axis by the contour $\Gamma$. The solution $v_0^{(-)}$ defined in this way will correspond to the asymmetric solution $u_-$ of (6) which tends to zero exponentially for $x \to \infty$ in the positive direction. This solution has a core that is very close to the core of the minimal tail symmetric solution, but it generally diverges in the negative $x$ direction. There is a similar $v_0^{(+)}$ function that can be calculated using (84), which is valid for $-\pi \le \arg(q) < 0$. Then the contour $\Gamma$ has to be chosen in the quadrant $\pi/2 < \arg(s) \le \pi$. This solution is obviously the conjugated mirror image of the previous one with

respect to the imaginary $q$ axis. The difference of $v_0^{(-)}$ and $v_0^{(+)}$ can be calculated using the residue theorem. The dominant contribution to the difference will be given by the singularity at $s = i$, so we need to determine the behavior of $V'_0(s)$ close to there. This can be inferred from the large $n$ behavior of the coefficients $a_n^{(0)}$ in (85).

The residue at $s = i$ of the function $I_0(s)$ will be determined by substituting the leading order part of (89), which corresponds to $G_0(n) = 1$. The result can be summed,

$$ \begin{aligned} I_0(s) &\approx \exp(-sq) K \sum_{n=0}^{\infty} (-1)^n s^{2n+1} \\ &= \exp(-sq) \frac{Ks}{1+s^2}. \end{aligned} \quad (94) $$

Since the residue of $s/(1+s^2)$ is $1/2$, it follows that

$$ \operatorname{Res}_{s=i} I_0(s) = \exp(-iq) \frac{K}{2}. \quad (95) $$

The residues at the singularities $s = ni$, where $n \ge 1$ integer, will be proportional to $\exp(-inq)$; hence we can neglect their contributions. The residue theorem can be applied by choosing a curve going from $s = 0$ to infinity in the domain $0 < \arg(s) < \pi/2$ and coming back in the region $\pi/2 < \arg(s) \le \pi$. The difference of the two functions is

$$ v_0^{(-)} - v_0^{(+)} = \pi i K \exp(-iq), \quad (96) $$

which is valid for any $q$ satisfying $\text{Imq} < 0$.

Grimshaw and Joshi also define a third function, $v_0^{(m)}$, by setting the contour $\Gamma$ exactly as the upper half of the imaginary $s$ axis, running through all the singularities there. This way the integral will take half of the pole contributions, yielding

$$ v_0^{(-)} - v_0^{(m)} = \frac{1}{2} \pi i K \exp(-iq). \quad (97) $$

Since all the $a_n^{(0)}$ are real, from (84) and (85) follows that the function $v_0^{(m)}$ has no imaginary part on the imaginary $q$ axis; hence it corresponds to a symmetric $u$ solution on the real $x$ axis. Actually, since $K$ is real, it belongs to the one with minimal tail amplitude. Larger tail amplitude symmetric solutions could be obtained by adding $\alpha_r \exp(-iq)$ to $v_0^{(-)}$ with arbitrary real $\alpha_r$. Taking the imaginary part of (97) on the lower half of the imaginary axis we obtain that

$$ \operatorname{Im} v_0^{(-)} = \frac{1}{2} \pi K \exp(-iq) \quad \text{for } \text{Req} = 0, \quad \operatorname{Im} q < 0. \quad (98) $$
The behavior of $Imv_0^{(-)}$ along the lower part of the imaginary axis for large $|q|$ is determined by the constant K given in (93). This result will allow us to determine the minimal tail amplitude for symmetric solutions in the next sections.

### C. Complex extension of the linear correction

All functions $A_n$ defined in Sec. V contains various powers of sech($\gamma x$); hence the linear correction $u_w$ is singular at the same places as the original u solution. Close to the singularity at $x = i\pi/(2\gamma)$ we substitute (72) for $x$ into (63) to obtain a function depending on $q$. Expressing the sine function as a difference of two exponentials, we can neglect the small term proportional to $\exp[-k\pi/(2\gamma\epsilon)]$; hence

$$ u_w = \frac{i\beta}{2} \exp\left(\sum_{\substack{n=2 \\ \text{even}}}^{\infty} A_n \epsilon^n\right) \\ \times \exp\left(\frac{k\pi}{2\gamma\epsilon} - ikq + i\delta_w + i\sum_{\substack{n=1 \\ \text{odd}}}^{\infty} \tilde{A}_n \epsilon^n\right). \quad (99) $$

Using (54) to bring back the original $A_n$, we can write this into the simpler form, using a single summation for all $n$,

$$ u_w = \frac{i\beta}{2} \exp\left(\frac{k\pi}{2\gamma\epsilon}\right) \exp(-iq + i\delta_w) \\ \times \exp\left(-\frac{(k-1)\pi}{2\gamma\epsilon} + \sum_{n=1}^{\infty} A_n \epsilon^n\right). \quad (100) $$

We keep $k$ as it is in the first exponential term, but substituting the form (16) of $k$ into the third term we see that $(k-1)/\epsilon$ is small. Hence we proceed first by expanding the argument of the third exponential in powers of $1/q$, and then we expand the exponential of the result in $\epsilon$. The result correct up to order $\epsilon^4$ is

$$ u_w = \frac{i\beta}{2} \exp\left(\frac{k\pi}{2\gamma\epsilon}\right) \exp(-iq + i\delta_w) \\ \times [(1 + 5\gamma^2\epsilon^2)Q_0(q) + \gamma^4\epsilon^4 Q_2(q) \dots], \quad (101) $$

where

$$ Q_0(q) = 1 + \frac{6i}{q} - \frac{33}{q^2} - \frac{237i}{q^3} + \frac{1890}{q^4} + \frac{17028i}{q^5} \\ - \frac{167733}{q^6} + \dots, \quad (102) $$

$$ Q_2(q) = -\frac{2iq^3}{15} - \frac{q^2}{5} + \frac{39iq}{5} - 25 + \frac{234i}{q} - \frac{14343}{5q^2} \\ - \frac{181119i}{5q^3} + \dots. \quad (103) $$

We should like to point out that Eq. (101) disagrees with Eq. (48) of Ref. [22], where the factor $(1 + 5\gamma^2\epsilon^2)$ is missing. The truncated versions of the asymptotic expansions $Q_n(q)$ correspond to complex extensions of symmetric functions. Hence, for $\delta_w = 0$ the function $u_w$ is purely imaginary on the imaginary $q$ axis, corresponding to the complex extension of an antisymmetric function. If $\delta_w = \pi/2$, then $u_w$ is real on the imaginary axis, and it corresponds to a symmetric function. Remember that $\beta$ and $\delta_w$ may also have $\epsilon$ dependence.

### D. Amplitude up to third order

There are two important solutions of Eq. (6), the symmetric solution $u_m$ that has minimal tail, and the asymmetric solution $u_-$ that has no tail for $x > 0$. They both can be calculated by numerical methods, but no accurate analytical solutions are known for either of them. However, their difference can be determined very precisely by the earlier presented higher order WKB method. The difference $u_m - u_-$ is exponentially small in terms of $\epsilon$, not just in the tail but also in the core domain, and even on the complex $x$ plane, including the matching region near the singularity. This allows us to represent the difference with the linearized solution $u_w$ in (64) with appropriate tail-amplitude $\beta$ and phase $\delta_w$. The symmetric solution $u_m$ has the oscillating tail given by (68), with amplitude $\alpha_m$ and phase $\delta_m$. The tail for $x > 0$ can be compensated by $u_w$ to obtain the asymmetric solution,

$$ u_{-} = u_{m} - u_{w} \quad \text{for } \beta = \alpha_{m}, \quad \delta_{w} = 0. \quad (104) $$

This equation is also valid on the complex $q$ plane where $u_w$ is given by (101).

Although the difference $u_w$ is relatively small with respect to $u_m$ in the core region and on its complex extension, there is a place where it can be clearly observed. Because of its symmetry, $u_m$ is purely real on the imaginary axis, on both the complex $x$ and $q$ planes. Hence it follows from (104) that the imaginary parts of $u_-$ and $-u_w$ have to agree there,

$$ \mathrm{Im}u_{-} = -\mathrm{Im}u_{w} \quad \text{for } \mathrm{Re}q = 0, \quad \mathrm{Im}q < 0, \quad (105) $$

where $\beta = \alpha_m$ and $\delta_w = 0$ in the form (101) of $u_w$. This equality is valid to any order in $\epsilon$ and $1/q$. The solution of the inner problem discussed earlier can be used to obtain the imaginary part of $u_-$ on the imaginary axis. The importance of (105) is that after this we can directly obtain the minimal tail amplitude $\alpha_m = \beta$ as well.

To leading order in $\epsilon$ we can use (98) to determine the imaginary part of $u_-$ on the imaginary axis. Using $v = \epsilon^2 u$ in (73) and $v \approx v_0$ in (74) we obtain that
$$ \mathrm{Im}u_{-} = \frac{\pi K}{2\epsilon^2} \exp(-iq)(1 + O(\epsilon^2)) \quad \text{for } \mathrm{Req} = 0, \quad \mathrm{Im}q < 0. \tag{106} $$

Comparing this with $Imu_w$ according to (105), considering the leading order part of (101), since $\alpha_m = \beta$, we obtain the expression for the minimal tail amplitude,

$$ \alpha_m^{(k,1)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{k\pi}{2\gamma\epsilon}\right). \tag{107} $$

The notation $k$ in the upper index indicates that according to (16) we keep $k = \sqrt{1 + 4\gamma^2\epsilon^2}$ inside the exponential. Writing 1 in the upper index shows that this result is precise to linear order in $\epsilon$, since the expression is valid up to a factor of $(1 + O(\epsilon^2))$.

To leading order in $\epsilon$ we can substitute $k = 1$, and we obtain the result first derived by Pomeau et al. in [21],

$$ \alpha_m^{(0)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{\pi}{2\gamma\epsilon}\right). \tag{108} $$

In Eq. (21) of [21] the double tail amplitude has been calculated for the asymmetric solution which tends to zero exponentially for $x \to -\infty$. Furthermore, there is an additional unnecessary factor of 2 there because of a 1/2 lost earlier, and an obviously missing fraction slash in the exponential. In Fig. 4 we compare the various order analytic results $\alpha_m^{(j)}$ in this section to the numerical amplitude $\alpha_m$ calculated by the high precision spectral method. We plot logarithmically the relative difference $\Delta\alpha_m^{(j)} = |(\alpha_m - \alpha_m^{(j)})/\alpha_m|$ as a function of $\epsilon$. Since the numerical result is much more accurate in these cases, $\Delta\alpha_m^{(j)}$ shows the relative error of the jth order analytic expansion result.

<figure><img src="image_7.png" /><figcaption>FIG. 4. Log-log plot of $\Delta\alpha_m^{(j)} = |(\alpha_m - \alpha_m^{(j)})/\alpha_m|$, showing the relative difference of the numerically calculated $\alpha_m$ from its various order analytic approximations $\alpha_m^{(j)}$ up to order three in $\epsilon$.</figcaption></figure>

With the $k$ factor included in the exponential, the result $\alpha_m^{(k,1)}$ for the tail amplitude in (107) corresponds to Eq. (56) of Grimshaw and Joshi [22]. Since

$$ \exp\left(-\frac{k\pi}{2\gamma\epsilon}\right) = \exp\left(-\frac{\pi}{2\gamma\epsilon}\right) \left[ 1 - \pi\gamma\epsilon + \frac{\pi^2}{2}\gamma^2\epsilon^2 - \left(\frac{\pi^2}{6} - 1\right)\pi\gamma^3\epsilon^3 + \left(\frac{\pi^2}{24} - 1\right)\pi^2\gamma^4\epsilon^4 - \left(\frac{\pi^4}{120} - \frac{\pi^2}{2} + 2\right)\pi\gamma^5\epsilon^5 + O(\epsilon^6) \right], \tag{109} $$

we can see that $\alpha_m^{(k,1)}$ contains odd powers of $\epsilon$ in its expansion. This shows that the inclusion of $k$ in the exponential improves the result to make it valid to linear $\epsilon$ order. Since (74) and (101) contain only even powers of $\epsilon$, odd powers in the expansion of $\alpha_m$ can only come from $k$ in the exponential. However, we will show below that contrary to the claim in [22], the result (107) is not valid to $\epsilon^2$ order.

Substituting the expansion (109) into (107) we can obtain an alternative first-order result for the amplitude,

$$ \alpha_m^{(1)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{\pi}{2\gamma\epsilon}\right) (1 - \pi\gamma\epsilon). \tag{110} $$

Surprisingly, according to Fig. 4, $\alpha_m^{(1)}$ has much lower relative error than the other first order result $\alpha_m^{(k,1)}$. The reason for this is that the $\epsilon$ expansion of $\alpha_m^{(k,1)}$ contains an $\epsilon^2$ term with coefficient $\pi^2\gamma^2/2 \approx 4.93\gamma^2$ which is much larger than the correct coefficient, which turns out to be $\pi^2\gamma^2/2 - 5\gamma^2 \approx -0.07\gamma^2$ as we will see a bit later.

As we have already seen, Eq. (80) for the second order inner problem has an appropriate exact solution, $v_1 = 2/3$. This has no imaginary part on the imaginary $q$ axis. To this approximation (74) gives $v \approx v_0 + \gamma^2\epsilon^2 v_1$; hence (106) is also valid to $\epsilon^2$ order,

$$ \mathrm{Im}u_{-} = \frac{\pi K}{2\epsilon^2} \exp(-iq)(1 + O(\epsilon^4)) \quad \text{for } \mathrm{Req} = 0, \quad \mathrm{Im}q < 0. \tag{111} $$

On the other hand, Eq. (101) clearly has an $\epsilon^2$ part, because of the $(1 + 5\gamma^2\epsilon^2)$ factor. In order to make (105) valid to order $\epsilon^2$ we have to cancel this contribution by an $\epsilon$ dependent factor in the $\beta = \alpha_m$ amplitude. This way we obtain a higher order generalization of (107) for the minimal amplitude,

$$ \alpha_m^{(k,3)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{k\pi}{2\gamma\epsilon}\right) (1 - 5\gamma^2\epsilon^2). \tag{112} $$

This approximation is correct to $\epsilon^3$ order, since the next correction would be an $\epsilon^4$ term in the factor multiplying the
exponential term. As we can see in Fig. 4, the result $\alpha_m^{(k,3)}$ is significantly more precise than the lower order approximations. For the lowest $\epsilon$ values considered, it gives the amplitude to seven digits of precision. The amplitude $\alpha_m^{(k,3)}$ is the corrected version of the result, Eq. (56), in Ref. [22]. That expression is valid only to $\epsilon^1$ order. The correctness of this additional $5\gamma^2\epsilon^2$ term in (112) is clearly supported by the numerical simulations. Actually, the reported differences between the analytical amplitude in [22] and the numerical results of Boyd in [24] are due to this missing term.

Using the expansion (109) we can obtain the following second and third order results for the amplitude:

$$ \alpha_m^{(2)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{\pi}{2\gamma\epsilon}\right) \left[1 - \pi\gamma\epsilon + \left(\frac{\pi^2}{2} - 5\right)\gamma^2\epsilon^2\right], \quad (113) $$

$$ \alpha_m^{(3)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{\pi}{2\gamma\epsilon}\right) \left[1 - \pi\gamma\epsilon + \left(\frac{\pi^2}{2} - 5\right)\gamma^2\epsilon^2 - \left(\frac{\pi^2}{6} - 6\right)\pi\gamma^3\epsilon^3\right]. \quad (114) $$

Strangely, as we can see in Fig. 4, $\alpha_m^{(2)}$ gives a slightly less precise approximation than the lower order result $\alpha_m^{(1)}$ given earlier in (110). The reason for this is that for the presented range of $\epsilon$ the neglected third order term is still larger than the second order one. Their orders become reversed only below $\epsilon \approx 0.00477$, which is not accessible by our numerical code. However, when the $\epsilon^3$ term is included, we obtain $\alpha_m^{(3)}$, which, according to the figure, is almost as precise as $\alpha_m^{(k,3)}$ in (112). As can be expected, all $\Delta\alpha_m^{(j)}$ and $\Delta\alpha_m^{(k,j)}$ tend to zero as $\epsilon^{j+1}$ in the figure, except $\Delta\alpha_m^{(1)}$ which decays faster in the presented $\epsilon$ interval. In the next section we will determine the $\epsilon^4$ order correction to the amplitude in (112). The fact that we can obtain even higher order results gives a very strong support that the results in the present section are correct.

The minimal amplitude tail has the asymptotic phase $\delta_m$. As we have already seen, the tail amplitude $\alpha$ of general symmetric solutions with asymptotic tail phase $\delta$ can be calculated from the minimal amplitude $\alpha_m$ by (70).

# VII. FOURTH ORDER

## A. Inner problem

In this case the equations we have to solve follow from the substitution of the expansion (74) of $v$ into the inner equation (78). The first two equations have already been given in (79) and (80). After substituting $v_1 = 2/3$ the next equation is

$$ v_{2qqqq} + v_{2qq} + 6v_0v_2 - 16v_0 - \frac{4}{3} = 0. \quad (115) $$

The function $v_0$ is assumed to be known here. Although we do not know an exact solution for $v_0$, we can use the power series asymptotic expansion (81) which starts according to (75). There may also be an exponentially small correction, corresponding to (97), and given in more detail in (A2) of Appendix.

Since (77) shows that the expansion of $v_2$ consists only of even powers of $q$, and starts with a $q^2$ term, we look for solutions in the form

$$ v_2 = \sum_{j=-1}^{\infty} b_j^{(2)} q^{-2j}. \quad (116) $$

Substituting into (115) we obtain that $b_{-1}^{(2)} = -2/15$, and for $n \ge 0$ we get

$$ (2n)(2n+1)b_n^{(2)} + (2n-2)(2n-1)(2n)(2n+1)b_{n-1}^{(2)} + 6 \sum_{j=-1}^{n} b_j^{(2)} b_{n-j+1}^{(0)} - 16b_{n+1}^{(0)} = 0. \quad (117) $$

Separating the last term from the summation we can obtain a recursion relation, similar to (83). The coefficients $b_n^{(2)}$ that we get agree with those already given in (77). This shows that we match the inner problem appropriately to the outer problem.

## B. Laplace transform

We intend to apply the Laplace transform method that we have already used for the zeroth order in (84). Our aim is to calculate the imaginary part of the asymmetric function $v_2^{(-)}$ on the imaginary $q$ axis. However, for the $j \ge 0$ integer we cannot obtain $q^j$ as the Laplace transform of a smooth function. Hence we have to separate those terms, defining $\tilde{v}_2$ by

$$ v_2 = -\frac{2}{15}q^2 + \frac{2}{3} + \tilde{v}_2. \quad (118) $$

We can now look for $\tilde{v}_2$ in the form of the Laplace transform of a function $V_2'(s)$,

$$ \tilde{v}_2 = \int_{\Gamma} I_2(s) ds, \quad I_2(s) = \exp(-sq)V_2'(s), \quad (119) $$

where the contour $\Gamma$ is from $s = 0$ to infinity, satisfying $\text{Re}(sq) > 0$. Similar to (85) and (86), from the identity for the Laplace transform of powers of $s$ it follows that

$$ V_2'(s) = \sum_{n=0}^{\infty} a_n^{(2)} s^{2n+1}, \quad (120) $$

where
$$a_n^{(2)} = \frac{b_{n+1}^{(2)}}{(2n+1)!}. \qquad (121)$$

Since $b_1^{(2)} = 64$, it follows that $a_0^{(2)} = 64$. Substituting into (117), for $n \ge 1$ we get

$$a_n^{(2)} + a_{n-1}^{(2)} + \frac{6}{(2n+3)!} \sum_{j=0}^{n} (2j+1)!(2n-2j+1)! a_j^{(2)} a_{n-j}^{(0)} \\ - \frac{4}{5}(2n+5)(2n+4)a_{n+2}^{(0)} - 12a_{n+1}^{(0)} = 0. \qquad (122)$$

The recursion relation is

$$\begin{aligned} & \frac{(n+3)(2n-1)}{(n+1)(2n+3)} a_n^{(2)} \\ & = -a_{n-1}^{(2)} - \frac{6}{(2n+3)!} \sum_{j=0}^{n-1} (2j+1)!(2n-2j+1)! a_j^{(2)} a_{n-j}^{(0)} \\ & \quad + \frac{4}{5}(2n+5)(2n+4)a_{n+2}^{(0)} + 12a_{n+1}^{(0)}. \end{aligned} \qquad (123)$$

For large $n$, the leading order behavior is $a_n^{(2)} \sim n^3(-1)^n$. It follows that the series for $V'_2(s)$ is convergent. The singularities are again located at $s = \pm ni$, where $n$ is any positive integer, but not at $s = 0$. To study the function near the $s = i$ singularity we need to study the large $n$ behavior of the coefficients $a_n^{(2)}$. We look for the large $n$ behavior in the form

$$a_n^{(2)} = (-1)^n \tilde{a}_n^{(2)}, \qquad \tilde{a}_n^{(2)} = \sum_{l=-3}^{\infty} \frac{g_l^{(2)}}{n^l}. \qquad (124)$$

For $\tilde{a}_n^{(0)}$ we use the expansion (89) with the already calculated coefficients given in (92). We substitute these into the equation that follows from (122),

$$\tilde{a}_n^{(2)} - \tilde{a}_{n-1}^{(2)} + \sum_{j=0}^{j_m} (-1)^j W_{n,j}^{(0)} \tilde{a}_{n-j}^{(2)} + \sum_{j=0}^{j_m} (-1)^j W_{n,j}^{(2)} \tilde{a}_{n-j}^{(0)} \\ - \frac{4}{5}(2n+5)(2n+4)\tilde{a}_{n+2}^{(0)} + 12\tilde{a}_{n+1}^{(0)} = 0, \qquad (125)$$

where $W_{n,j}^{(0)}$ is defined in (91), $j_m$ is the positive integer used there, and

$$W_{n,j}^{(2)} = \frac{6}{(2n+3)!} (2j+1)!(2n-2j+1)! a_j^{(2)}. \qquad (126)$$

Equation (125) has a homogeneous part that can be obtained by setting $\tilde{a}_n^{(0)} = 0$. Note that the value of $W_{n,j}^{(0)}$ is nonzero even in this case. The homogeneous part agrees exactly with Eq. (90); hence, it has the solution $G_0(n)$

multiplied by an arbitrary constant. It follows that the general solution of (125) can be written as

$$\tilde{a}_n^{(2)} = K(G_2(n) + K_2 G_0(n)), \qquad (127)$$

where $K_2$ is an arbitrary constant and $G_2(n)$ is a particular solution. Its expansion can be obtained using some algebraic manipulation software, similar to (92),

$$G_2(n) = \frac{16}{15}n^3 + \frac{28}{5}n^2 + \frac{368}{15}n - \frac{132}{n} + \frac{4122}{5n^2} - \frac{50833}{10n^3} \\ + \frac{3144915}{112n^4} + \dots \qquad (128)$$

We have chosen the unique particular solution $G_2(n)$ for which there is no $n^0$ term. We can obtain the value of the constant $K_2$ by calculating $a_n^{(2)}$ using (123) for large $n$, similar to what we did for $K$ in (93). Using $K_2 \approx [a_n^{(2)}(-1)^n/K - G_2(n)]/G_0(n)$ we get

$$K_2 = -36.544068193583744293. \qquad (129)$$

### C. Asymmetric solution

Our aim is to calculate the imaginary part of the asymmetric $v_2 \equiv v_2^{(-)}$ solution on the lower part of the imaginary $q$ axis. Obviously, this will be the same as the imaginary part of $\tilde{v}_2 \equiv \tilde{v}_2^{(-)}$, since they differ only by two terms according to (118). We apply the same reasoning as we did for $v_0^{(-)}$ in the paragraph before Eq. (94). The function $\tilde{v}_j^{(-)}$ can be calculated in the domain $-\pi < \arg(q) \le 0$ by the integral (119), with a contour in the region satisfying $0 < \arg(s) < \pi/2$. The imaginary part of $\tilde{v}_2^{(-)}$ on the imaginary axis can be determined to leading order using the residue of the function $I_2(s)$ at $s = i$. Substituting (120), (124), and (127) we get

$$I_2(s) \approx \exp(-sq) K \sum_{n=0}^{\infty} (-1)^n [G_2(n) + K_2 G_0(n)] s^{2n+1}. \qquad (130)$$

From the expansion of $G_0(n)$ and $G_2(n)$ only terms of the form $(-1)^n n^j s^{2n+1}$ with $j \ge 0$ will give contributions to this residue. These terms can be summed, for example, for $j = 1$,

$$\sum_{n=0}^{\infty} (-1)^n n s^{2n+1} = - \frac{s^3}{(1+s^2)^2}. \qquad (131)$$

The sum can be calculated for other concrete $j \ge 0$ values, but we could not find a general formula valid for arbitrary $j$. However, the general residue turns out to be very simple,
$$ \mathrm{Res}_{s=i}^{\infty} (-1)^n n^j s^{2n+1} = \frac{1}{2} (-1)^j \quad \text{for } j \ge 0. \quad (132) $$

Hence from (92) and (128) we get

$$ \begin{aligned} \mathrm{Res}_{s=i} I_2(s) &= \exp(-iq) \frac{K}{2} \left( -\frac{16}{15} + \frac{28}{5} - \frac{368}{15} + K_2 \right) \\ &= \exp(-iq) \frac{K}{2} (K_2 - 20). \end{aligned} \quad (133) $$

Next we define a symmetric function $v_2^{(m)}$ by taking the path $\Gamma$ along the upper half of the imaginary axis. The difference of $v_2^{(-)}$ and $v_2^{(m)}$ can be calculated similar to Eq. (97). The residue theorem can be applied by choosing a curve going from $s=0$ to infinity in the domain $0 < \arg(s) < \pi/2$ and coming back along the upper half of the imaginary $s$ axis. This way we have to take into account half of the residue at $s=i$, so we get

$$ v_2^{(-)} - v_2^{(m)} \approx \frac{1}{2} \pi i K \exp(-iq) (K_2 - 20). \quad (134) $$

Here we write $\approx$ because we have not considered yet the change in $v_2$ due to the modification of $v_0$ from $v_0^{(-)}$ to $v_0^{(m)}$ in (115). This will provide a finite number of terms proportional to $q^n \exp(-iq)$ with $n > 0$ integers. The coefficients of those terms are completely fixed already by the leading order constant $K$, which determines $v_0^{(-)} - v_0^{(m)}$ according to (97). The detailed analysis of the perturbations of Eq. (115) is given in Appendix. The change in $v_2$ due to the modification of $v_0$ can be chosen in a way that it does not influence the terms $q^0 \exp(-iq)$ calculated in (134).

Analogously to (98), taking the imaginary part we obtain that on the lower part of the imaginary axis

$$ \mathrm{Im} v_2^{(-)} \approx \frac{1}{2} \pi K \exp(-iq) (K_2 - 20). \quad (135) $$

Here and in the following few equations the approximate equation sign indicates that we only consider the $q^0 \exp(-iq)$ parts of the expressions. Using $v = \epsilon^2 u$ and (74) we obtain a higher order generalization of (111),

$$ \mathrm{Im} u_{-} \approx \frac{\pi K}{2\epsilon^2} \exp(-iq) [1 + (K_2 - 20)\gamma^4 \epsilon^4 + O(\epsilon^6)] \quad (136) $$

for $\mathrm{Re}q = 0$, $\mathrm{Im}q < 0$.

### D. Tail amplitude up to fifth order

We generalize the expression for minimal amplitude $\alpha_m$ in (112) by allowing a fourth order contribution,

$$ \alpha_m^{(k,5)} = -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{k\pi}{2\gamma\epsilon}\right) (1 - 5\gamma^2\epsilon^2 - \xi_2\gamma^4\epsilon^4), \quad (137) $$

where $\xi_2$ is a constant that will be determined below. This approximation is correct to $\epsilon^5$ order, since the next correction would be proportional to $\epsilon^6$. Substituting $\beta = \alpha_m^{(k,5)}$ and $\delta_w = 0$ into (101) and only considering the $q^0 \exp(-iq)$ terms, we get

$$ u_w \approx -\frac{i\pi K}{2\epsilon^2} \exp(-iq) [1 - (\xi_2 + 50)\gamma^4 \epsilon^4 + O(\epsilon^6)]. \quad (138) $$

According to (105), we compare the imaginary part on the axis with (136) to obtain $K_2 - 20 = -\xi_2 - 50$, which gives

$$ \xi_2 = -K_2 - 30 \approx 6.544068193583744293. \quad (139) $$

In Fig. 5 we show the relative error of the analytic results $\alpha_m^{(j)}$ obtained in this section when compared to the precise numerical amplitude $\alpha_m$. As can be expected, the relative error of the jth order results tend to zero as $\epsilon^{j+1}$.

Using the expansion (109) we can obtain the following alternative fifth order result for the amplitude:

$$ \begin{aligned} \alpha_m^{(5)} = & -\frac{\pi K}{\epsilon^2} \exp\left(-\frac{\pi}{2\gamma\epsilon}\right) \left[1 - \pi\gamma\epsilon + \left(\frac{\pi^2}{2} - 5\right)\gamma^2\epsilon^2 \right. \\ & \left. - \left(\frac{\pi^2}{6} - 6\right)\pi\gamma^3\epsilon^3 + \left(\frac{\pi^4}{24} - \frac{7\pi^2}{2} - \xi_2\right)\gamma^4\epsilon^4 \right. \\ & \left. - \left(\frac{\pi^4}{120} - \frac{4\pi^2}{3} + 7 - \xi_2\right)\pi\gamma^5\epsilon^5 \right]. \end{aligned} \quad (140) $$

We can also define a fourth order amplitude $\alpha_m^{(4)}$ by dropping the $\epsilon^5$ term from $\alpha_m^{(5)}$. The relative error of these approximations are also shown in Fig. 5. The two different

<figure><img src="image_8.png" /><figcaption>FIG. 5. Log-log plot of $\Delta\alpha_m^{(j)} = |(\alpha_m - \alpha_m^{(j)})/\alpha_m|$, showing the relative difference of the precise numerical $\alpha_m$ from various order analytic results $\alpha_m^{(j)}$, up to order five in $\epsilon$.</figcaption></figure>
fifth order results are so close to each other that $\Delta\alpha_m^{(5)}$ and $\Delta\alpha_m^{(k,5)}$ are indistinguishable in the logarithmic figure. Actually, $\Delta\alpha_m^{(k,5)}$ is larger by less than 2% than $\Delta\alpha_m^{(5)}$. Both $\alpha_m^{(5)}$ and $\alpha_m^{(k,5)}$ are correct to nine digits of precision for the smallest $\epsilon$ value shown.

## VIII. CONCLUSIONS

We have considered the problem to compute corrections to the beyond-all-orders small amplitude standing wave tails of weakly localized soliton solutions in a fifth order KdV equation—the KdV equation with an added fifth order dispersion term, $\epsilon^2\partial_x^5$. The fKdV equation is not only of some genuine physical interest, but at the same time it also serves for us as a simplified model to prepare the ground for the significantly more complicated problem for oscillon-quasibreathers in field theories. For $\epsilon \ll 1$ the simplest stationary solution is a bounded, one-parameter family of KdV-type solitons which are weakly localized due to an asymptotic standing wave tail, tending to the KdV one-soliton for $\epsilon \to 0$. These solutions are symmetric, and they are characterized by the asymptotic phase, $\delta$, of their tail. Our main analytical result for the physically relevant minimal amplitude wave tail can be succinctly presented as

$$ \alpha_{\min} = \frac{\lambda}{\epsilon^2} e^{\frac{k(\epsilon)\Gamma}{2\epsilon}} (1 - 5\gamma^2\epsilon^2 - \zeta_4\gamma^4\epsilon^4 + \mathcal{O}(\epsilon^6)), \\ \lambda \approx 19.9689\pi, \quad \zeta_4 \approx 6.5441, $$

where $\gamma$ parametrizes the speed of the unperturbed KdV soliton and the wave number, $k(\epsilon) = (1 + 4\gamma^2\epsilon^2)^{1/2}$, has to be expanded up to $\mathcal{O}(\epsilon^5)$ [see Eq. (140)]. Our paper resolves a long-standing discrepancy for the $\mathcal{O}(\epsilon^2)$ coefficient between the result of Ref. [22] and the numerical results of Ref. [24].

The asymptotic phase of the minimal amplitude tail, $\delta_{\min}$, has been determined in a WKB approximation to rather high (~100) orders in $\epsilon$. We have also found that for a given $\epsilon$, the tail amplitude, $\alpha$, of any member of the one-parameter family, characterized with phase $\delta$, is related to the minimal tail amplitude $\alpha_{\min}$ by $\alpha = \alpha_{\min}/\cos(\delta - \delta_{\min})$.

We have also developed an efficient, arbitrary precision pseudospectral code to solve the fKdV equation, to investigate the asymptotic standing wave tails to high numerical precision. Our numerical code is fast enough so that we can employ numerical minimization procedures to find the value of the phase, $\delta_{\min}$, of the minimal amplitude solution very precisely, even for quite small values of $\epsilon$. The numerically obtained value for $\delta_{\min}$ agrees very well with the optimally truncated asymptotic series result (see Fig. 3). The remarkably good agreement of our higher order perturbative results for $\alpha_{\min}$ with our numerical calculations (see Figs. 4 and 5) gives very strong support that both our

analytic and numerical considerations are correct and reliable.

## ACKNOWLEDGMENTS

The research of G. F. has been supported in part by the National Research Development and Innovation Office (NKFIH) OTKA Grants No. K 138277 and No. K 142423.

## APPENDIX: LINEAR CORRECTION TO THE INNER PROBLEM

The complex extension of the exponentially small linear correction $u_w$ to the outer solution has been calculated in Sec. VI C. Similar small corrections can also be determined to $v_n$ at each order of the inner problem. We proceed by substituting $v_n \to v_n + w_n$ for all $n$ and linearizing for $w_n$.

From (79) we obtain

$$ w_{0qqqq} + w_{0qq} + 6v_0w_0 = 0. \quad (A1) $$

For $v_0$ we can use the asymptotic expansion (81) that starts with the terms given in (75). We search the solution for $w_0$ in the form

$$ w_0 = \frac{i}{2} b_0 \exp(-iq) F_0(q), \\ F_0(q) = 1 + \sum_{j=1}^{\infty} \frac{w_0^{(j)}}{q^j}, \quad (A2) $$

where $b_0$ and $w_0^{(j)}$ are complex constants. Substituting into (A1), the coefficients of the various powers of $1/q$ determine the constants $w_0^{(j)}$. The first few terms give

$$ F_0(q) = 1 + \frac{6i}{q} - \frac{33}{q^2} - \frac{237i}{q^3} + \frac{1890}{q^4} + \frac{17028i}{q^5} \\ - \frac{167733}{q^6} + \dots \quad (A3) $$

This naturally agrees with the function $Q_0(q)$ given in (102). Since Eq. (A1) is linear, the constant $b_0$ can be arbitrary. Real $b_0$ corresponds to a function antisymmetric on the real axis, while purely imaginary $b_0$ to a symmetric one.

The linearization of the $\epsilon^2$ order inner equation (80) gives

$$ w_{1qqqq} + w_{1qq} + 6v_0w_1 + 2(3v_1 - 2)w_0 = 0. \quad (A4) $$

For $v_1 = 2/3$ this equation becomes the same as (A1), and the solutions are

$$ w_1 = \frac{i}{2} b_1 \exp(-iq) F_0(q), \quad (A5) $$

where $b_1$ is an arbitrary complex constant. Since $v_1 = 2/3$ satisfies the matching condition to the outer solution, and it
is symmetric, without an imaginary part for Req = 0, we do not need to perturb $v_1$, we can set $b_1 = 0$. On the other hand, the $\epsilon^2$ part of (101) contains a perturbation corresponding to (A5) with $b_1 = 5b_0$. We had to cancel these terms by a $1-5\gamma^2\epsilon^2$ factor in the amplitude $\alpha_m$ in (112). Linearizing (115) we obtain

$$w_{2qqqq} + w_{2qq} + 6v_0w_2 + 2(3v_2 - 8)w_0 = 0. \quad (A6)$$

The homogeneous part agrees with Eq. (A1) for $w_0$; hence we can always add an arbitrary constant times $\exp(-iq)F_0(q)$ to the solution. Searching the solution in the form

$$w_2 = \frac{i}{2} \exp(-iq) \sum_{j=-3}^{\infty} \frac{w_2^{(j)}}{q^j}, \quad (A7)$$

from the coefficients of the various powers of q we obtain

$$w_2 = \frac{i}{2} \exp(-iq) [b_0 F_2(q) + b_2 F_0(q)], \quad (A8)$$

where $b_0$ is the constant in (A2), $b_2$ is an arbitrary complex constant, the function $F_0(q)$ is given in (A3), and

$$F_2(q) = -\frac{2i}{15}q^3 - \frac{1}{5}q^2 + \frac{39i}{5}q + \frac{384i}{q} - \frac{18468}{5q^2} - \frac{210744i}{5q^3} + \dots \quad (A9)$$

The function $F_2(q)$ was made unique here by setting the coefficient of the $q^0$ term zero. The $\epsilon^4$ part of (101), which is proportional to the $Q_2(q)$ given in (103), can be obtained from (A8) by setting $b_2 = -25b_0$.

[1] E. J. Copeland, M. Gleiser, and H.-R. Müller, Oscillons: Resonant configurations during bubble collapse, Phys. Rev. D **52**, 1920 (1995).

[2] E. P. Honda and M. W. Choptuik, Fine structure of oscillons in the spherically symmetric $\phi^4$ Klein-Gordon model, Phys. Rev. D **65**, 084037 (2002).

[3] D. Cuncynates and T. Giurgica-Tiron, Structure of the oscillator: The dynamics of attractive self-interaction, Phys. Rev. D **103**, 116011 (2021).

[4] G. Fodor, A review on radiation of oscillons and oscillatons, arXiv:1911.03340.

[5] L. Visinelli, Boson stars and oscillatons: A review, Int. J. Mod. Phys. D **30**, 2130006 (2021).

[6] H.-Y. Zhang, M. A. Amin, E. J. Copeland, P. M. Saffin, and K. D. Lozanov, Classical decay rates of oscillons, J. Cosmol. Astropart. Phys. **07** (2020) 055.

[7] J. Ollé, O. Pujolàs, and F. Rompineve, Recipes for oscillator longevity, J. Cosmol. Astropart. Phys. **09** (2021) 015.

[8] G. Fodor, P. Forgács, Z. Horváth, and M. Mezei, Computation of the radiation amplitude of oscillons, Phys. Rev. D **79**, 065002 (2009).

[9] G. Fodor, P. Forgács, Z. Horváth, and M. Mezei, Radiation of scalar oscillons in 2 and 3 dimensions, Phys. Lett. B **674**, 319 (2009).

[10] H. Segur and M. D. Kruskal, Nonexistence of Small-Amplitude Breather Solutions in $\phi^4$ Theory, Phys. Rev. Lett. **58**, 747 (1987).

[11] G. Fodor, P. Forgács, P. Grandclément, and I. Rácz, Oscillons and quasibreathers in the $\phi^4$ Klein-Gordon model, Phys. Rev. D **74**, 124003 (2006).

[12] P. M. Saffin and A. Tranberg, Oscillons and quasi-breathers in $D+1$ dimensions, J. High Energy Phys. **01** (2007) 030.

[13] T. Kakutani and H. Ono, Weak non-linear hydromagnetic waves in a cold collision-free plasma, J. Phys. Soc. Jpn. **26**, 1305 (1969).

[14] T. Kawahara, Oscillatory solitary waves in dispersive media, J. Phys. Soc. Jpn. **33**, 260 (1972).

[15] J. K. Hunter and J. Scheurle, Existence of perturbed solitary wave solutions to a model equation for water waves, Physica (Amsterdam) **32D**, 253 (1988).

[16] E. Benilov, R. Grimshaw, and E. Kuznetsova, The generation of radiating waves in a singularly-perturbed Korteweg-de Vries equation, Physica (Amsterdam) **69D**, 270 (1993).

[17] C. J. Amick and J. B. McLeod, A singular perturbation problem in water-waves, Stab. Appl. Anal. Contin. Media **1**, 127 (1991).

[18] B. T. N. Gunney, Y. A. Li, and P. J. Olver, Solitary waves in the critical surface-tension model, J. Eng. Math. **36**, 99 (1999).

[19] J. P. Boyd, A numerical calculation of a weakly non-local solitary wave: The $\phi^4$ breather, Nonlinearity **3**, 177 (1990).

[20] J. P. Boyd, Weakly Nonlocal Solitary Waves and Beyond-All-Orders Asymptotics: Generalized Solitons and Hyper-asymptotic Perturbation Theory (Springer US, Boston, MA, 1998).

[21] Y. Pomeau, A. Ramani, and B. Grammaticos, Structural stability of the Korteweg-de Vries solitons under a singular perturbation, Physica (Amsterdam) **31D**, 127 (1988).

[22] R. Grimshaw and N. Joshi, Weakly nonlocal solitary waves in a singularly perturbed Korteweg-De Vries equation, SIAM J. Appl. Math. **55**, 124 (1995).

[23] J. P. Boyd, Weakly non-local solitons for capillary-gravity waves: Fifth-degree Korteweg-de Vries equation, Physica (Amsterdam) **48D**, 129 (1991).
[24] J. P. Boyd, Multiple precision pseudospectral computations of the radiation coefficient for weakly nonlocal solitary waves: Fifth-order Korteweg–DeVries equation, *Comput. Phys. **9***, 324 (1995).

[25] J. Boyd, *Chebyshev and Fourier Spectral Methods: Second Revised Edition*, Dover Books on Mathematics (Dover Publications, New York, 2013).

[26] CLN—Class Library for Numbers, https://www.ginac.de/CLN/.

[27] Arb—a C library for arbitrary-precision ball arithmetic, https://arblib.org/.

[28] F. Johansson, Arb: Efficient arbitrary-precision midpoint-radius interval arithmetic, *IEEE Trans. Comput. **66***, 1281 (2017).

[29] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, *Numerical Recipes 3rd Edition: The Art of Scientific Computing* (Cambridge University Press, Cambridge, England, 2007).

[30] S. M. Sun, On the oscillatory tails with arbitrary phase shift for solutions of the perturbed Korteweg–de Vries equation, *SIAM J. Appl. Math. **58***, 1163 (1998).