# Stability of multi-parameter solitons: Asymptotic approach

## I Introduction

Solitary waves (’solitons’) can appear when an initial excitation applied to a medium is strong enough to cause nonlinear response. Formally solitons are solutions of some nonlinear partial differential equations and their dynamics generally is a complex phenomenon, which can be described exactly only in the very special integrable situations . The problems of soliton stability and instability induced dynamics in nonintegrable Hamiltonian models have paramount importance for understanding of a wide range of physical phenomena covering such fields as propagation of electromagnetic, water and plasma waves, condensed matter physics and classical field theory . Several analytical approaches to the stability problem are known. For instance, in the nearly integrable situations the perturbation theory based on the inverse scattering transform can be used . Far from the integrable limit variety of methods can be applied. Among them asymptotic stability theory , method of adiabatically varying soliton parameters , Lyapunov  and Evans  methods.

Generally, stability of a solitary wave in a Hamiltonian model can be lost due to bifurcations involving appearance of a positive eigenvalue (stationary instability) in the soliton spectrum or a pair of complex conjugate eigenvalues with positive real parts (oscillatory instability) . Both types of these instabilities have been extensively studied in the different solitonic contexts proving their ubiquitousness and fundamental importance, see Refs.  and Refs. , respectively, for the stationary and oscillatory instabilities. In the most of the known cases the loss of stability is associated with the collisions of the purely imaginary eigenvalues corresponding to the so called internal modes  of the soliton spectrum (see  for interesting exceptions).

Applying the above mentioned methods it was shown that in many cases a threshold of stationary instability of multi-parameter solitons is given by the zero of the determinant of the Jacobi matrix $`J_{ij}=_{\kappa _j}Q_i`$, where $`\kappa _j`$ are the soliton parameters and $`Q_i`$ are the associated motion integrals . The condition $`det(J_{ij})=0`$ is, in fact, the compatibility condition of the problem arising in the leading (zero) order of the asymptotic solution of the eigenvalue problem governing stability of the soliton . To find expressions for the eigenvalues it is necessary to proceed further and solve problems arising in higher (at least first) orders. Up to now this was done only for the specific class of model equations having single parameter soliton families . For stationary bifurcations of two-parameter solitons adiabatic method has been applied in Refs. . Linear approximation of this method actually gives an expression for eigenvalues, for more details see Section IV. However, all known developments of this method fail to give criterion indicating transition to the oscillatory instability, i.e. instability with complex eigenvalues. It is also difficult to extend this method beyond its first order because of the rather involved calculations.

The purpose of this work is to formulate a general asymptotic approach to stability of multi-parameter solitons in Hamiltonian models, to show how it can be used to find expressions for the instability growth rates with arbitrary accuracy and to formulate criterion for the oscillatory instability of solitons.

## II Model equations and symmetries

We will consider Hamiltonian equations in the form

$$i\frac{E_n}{z}=\frac{\delta H}{\delta E_n^{}},n=1,2\mathrm{}N,$$
(1)   
which describes a wide range of physical phenomena related with self-action and interaction of slowly varying wave envelopes in a variety of nonlinear media , for general review of the Hamiltonian formalism see . Here $`E_n`$ are complex fields, $`z`$ the propagation direction of the interacting waves, $`x`$ the coordinate characterizing dispersion or diffraction, $`H=H(_xE_n,E_n,_xE_n^{},E_n^{})`$ is the Hamiltonian and  means complex conjugation. We will assume that $`H`$ is invariant with respect to a set of $`(L1)`$ phase transformations:

$$E_nE_n\mathrm{exp}(i\gamma _{nl}\varphi _l),l=1,2,\mathrm{}(L1),$$
(2)   
$`\varphi _l`$ are arbitrary real phases and $`\gamma _{nl}`$ are some constants. Because $`H`$ does not depend on $`x`$ explicitly, Eqs. (1) are also invariant with respect to arbitrary translations along $`x`$:

$$E_n(x,z)E_n(xx_0,z).$$
(3)   
Symmetry properties (2), (3) together with Hamiltonian nature of our problem imply presence of $`L`$ conserved quantities, see, e.g., , which are the $`(L1)`$ energy invariants

$$Q_l=𝑑x\underset{n=1}{\overset{N}{}}\gamma _{nl}|E_n|^2,l=1,2,\mathrm{}(L1),$$
(4)   
and the momentum

$$Q_L=\frac{1}{2i}𝑑x\underset{n=1}{\overset{N}{}}(E_n^{}_xE_nE_n_xE_n^{}).$$
(5)   
Another important consequence of the invariances (2), (3) is that a certain class of solutions of Eqs. (1) can be sought in a form when $`x_0`$ and $`\varphi _l`$ are linear functions of $`z`$, i.e. $`x_0=\kappa _Lz`$ and $`\varphi _l=\kappa _lz`$, then

$$E_n(x,z)=a_n(x\kappa _Lz)\mathrm{exp}(i\underset{l=1}{\overset{L1}{}}\gamma _{nl}\kappa _lz),$$
(6)   
where $`\{\kappa _l\}_{l=1}^{L1}`$ and $`\kappa _L`$ are real parameters characterizing, respectively, phase velocities of the interacting waves and the soliton group velocity. Functions $`a_n(\tau )`$ obey a system of the ordinary differential equations

$$(i\kappa _L_\tau +\alpha _n)a_n=\frac{\delta H_a}{\delta a_n^{}},$$
(7)   
where $`H_aH(_\tau a_n,a_n,_\tau a_n^{},a_n^{})`$, $`\tau =x\kappa _Lz`$ and $`\alpha _n=_{l=1}^{L1}\gamma _{nl}\kappa _l`$. We assume now that in a certain domain of the parameter space $`(\kappa _1,\kappa _2,\mathrm{}\kappa _L)`$ Eqs. (7) have a family of the solitary solutions such that $`|a_n|0`$ for $`\tau \pm \mathrm{}`$.

## III Asymptotic stability analysis

To study stability of the solitons we seek solutions of Eqs. (1) in the form

$$E_n=(a_n(\tau )+\epsilon _n(\tau ,z))\mathrm{exp}(i\underset{l=1}{\overset{L1}{}}\gamma _{nl}\kappa _lz),$$
(8)   
where $`\epsilon _n(\tau ,z)`$ are small complex perturbations. Linearizing Eqs. (1) and assuming that $`\epsilon _n(\tau ,z)=\xi _n(\tau )e^{\lambda z}`$, $`\epsilon _n^{}(\tau ,z)=\xi _{n+N}(\tau )e^{\lambda z}`$ we get the following nonselfadjoint eigenvalue problem (EVP)

$$i\lambda \stackrel{}{\xi }=\widehat{}\stackrel{}{\xi }\left(\begin{array}{cc}\widehat{S}& \widehat{R}\\ \widehat{R}^{}& \widehat{S}^{}\end{array}\right)\stackrel{}{\xi },$$
(9)   
where $`\stackrel{}{\xi }=(\xi _1,\mathrm{}\xi _N,\xi _{N+1}\mathrm{}\xi _{2N})^T`$, and $`\widehat{R}`$, $`\widehat{S}`$ are $`N\times N`$ matrix operators with elements given by

$$\widehat{s}_{nl}=\delta _{nl}(\alpha _n+i\kappa _L_\tau )+\frac{\delta ^2H_a}{\delta a_n^{}\delta a_l},\widehat{r}_{nl}=\frac{\delta ^2H_a}{\delta a_n^{}\delta a_l^{}},$$

here $`\delta _{nl}`$ is the Kroneker symbol. Note, that the operator $`\widehat{S}`$ is a selfadjoint one, i.e. $`\widehat{S}=\widehat{S}^{}`$, and $`\widehat{R}`$ is a symmetric operator, i.e. $`\widehat{R}=\widehat{R}^T`$.

To solve EVP (9) we apply the asymptotic approach, which relies on expansion of the unknown eigenvector $`\stackrel{}{\xi }`$ into an asymptotic series near either neutral eigenmodes , i.e. zero-eigenvalue modes, of the operator $`\widehat{}`$, or modes of continuum , or both of them . The neutral modes can be generated by infinitesimal variations of the free parameters of the soliton and thus always be presented as explicit functions of the soliton solution. At the same time continuum eigenmodes are explicitly known in the very rare, normally in integrable, situations . This is an important fact which makes the asymptotic expansion near the neutral modes by the very practical tool of the stability theory. However, as any approximate method, it has a certain limitation. Namely, it describes only eigenvalues $`\lambda `$ corresponding to a specific class of the perturbations which in the zero approximation can be expressed as a linear superposition of the neutral eigenmodes. Thus, generally speaking, on the basis of this approach one can get only sufficient conditions for soliton instability or, in other words, necessary conditions for soliton stability. Therefore presence of other instabilities which can be captured only numerically can always be expected .

By infinitesimal variation of $`\varphi _l`$ and $`x_0`$ it can be shown that

$`\stackrel{}{u}_l=(\gamma _{1l}a_1,\mathrm{}\gamma _{Nl}a_N,\gamma _{1l}a_1^{},\mathrm{}\gamma _{Nl}a_N^{})^T,\stackrel{}{u}_L={\displaystyle \frac{\stackrel{}{a}}{\tau }},`$ (10)   
$`\stackrel{}{a}(a_1,\mathrm{}a_N,a_1^{},\mathrm{}a_N^{})^T,l=1,\mathrm{}(L1)`$ (11)   
are neutral modes of $`\widehat{}`$, i.e. $`\widehat{}\stackrel{}{u}_l=0`$ $`(l=1,\mathrm{}L)`$. $`\widehat{}`$ also has $`L`$ associated vectors $`\stackrel{}{U}_l=\stackrel{}{a}/\kappa _l`$ such that $`\widehat{}\stackrel{}{U}_l=\stackrel{}{u}_l,l=1,\mathrm{}L`$.

It is straightforward to see that any solution of EVP (9) must obey $`L`$ solvability conditions

$$\stackrel{}{w}_l|\lambda \stackrel{}{\xi }=0,l=1,\mathrm{}L,$$
(12)   
where $`\stackrel{}{y}|\stackrel{}{z}=_{i=1}^{2N}𝑑xy_i^{}z_i`$ and $`\stackrel{}{w}_l`$ are the neutral modes of the operator $`\widehat{}^{}`$, $`\widehat{}^{}\stackrel{}{w}_l=0`$,

$`\stackrel{}{w}_l=(\gamma _{1l}a_1,\mathrm{}\gamma _{Nl}a_N,\gamma _{1l}a_1^{},\mathrm{}\gamma _{Nl}a_N^{})^T,\stackrel{}{w}_L=i{\displaystyle \frac{\stackrel{}{b}}{\tau }},`$ (13)   
$`\stackrel{}{b}=(a_1,\mathrm{}a_N,a_1^{},\mathrm{}a_N^{})^T,l=1,\mathrm{}(L1).`$ (14)   
Associated vectors of $`\widehat{}^{}`$ are $`\stackrel{}{W}_l=\stackrel{}{b}/\kappa _l`$ and they obey $`\widehat{}^{}\stackrel{}{W}_l=\stackrel{}{w}_l,l=1,\mathrm{}L`$.

Close to instability threshold it is naturally to assume that $`|\lambda |ϵ1`$. As it was already discussed above we will consider a special class of the perturbations which in the leading approximation can be presented as a linear combination of the neutral modes. Therefore we seek an asymptotic solution of EVP (9) in the following form

$$\stackrel{}{\xi }=\underset{m=0}{\overset{\mathrm{}}{}}ϵ^m\stackrel{}{\xi }_m(x),\stackrel{}{\xi }_0=\underset{l=1}{\overset{L}{}}C_l\stackrel{}{u}_l$$
(15)   
where constants $`C_l`$ and vector-functions $`\stackrel{}{\xi }_{m>0}`$ have to be defined. Here and below $`l=1,2,\mathrm{}L`$. Substitution (15) into EVP (9) gives a recurrent system of equations for $`\stackrel{}{\xi }_m`$

$$\stackrel{}{\xi }_{m>0}=\left[\frac{i\lambda }{ϵ}\widehat{}^1\right]^m\stackrel{}{\xi }_0.$$
(16)   
Substituting (15), (16) into conditions (12) one will find the homogeneous system of the $`L`$ linear algebraic equations

$$\lambda ^2\stackrel{}{w}_l|\underset{m=0}{\overset{\mathrm{}}{}}(\lambda ^2)^m\widehat{}^{2m}\underset{l=1}{\overset{L}{}}C_l\stackrel{}{U}_l=0$$
(17)   
for $`L`$ unknown constants $`C_l`$. System (17) has a nontrivial solution providing that the corresponding determinant is equal to zero. This determinant is an infinite-order polynomial with respect to $`\lambda ^2`$, which, in fact, is the asymptotic expansion of Evans function . Zeros of this polynomial define the spectrum of the solitary wave linked with the chosen class of the perturbations. Thus the equation specifying eigenvalues $`\lambda `$ is

$$\lambda ^{2L}\underset{j=0}{\overset{\mathrm{}}{}}(\lambda ^2)^jD_j=0,$$
(18)   
where $`D_j`$ are the real constants. Eq. (18) always has zero root of the $`2L`$-order. It indicates that each of the zero eigenvalues corresponding to the neutral modes $`\stackrel{}{u}_l`$ is doubly degenerate one. This degeneracy originates from the presence of the associated vectors $`\stackrel{}{U}_l`$.

To write the explicit expressions for $`D_j`$ it will be convenient to introduce vectors $`\stackrel{}{}_l^{(m)}=(_{l1}^{(m)}\mathrm{}_{lL}^{(m)})`$, where,

$$_{ll^{^{}}}^{(m)}=\stackrel{}{w}_l|\widehat{}^{2m}\stackrel{}{U}_l^{^{}},m=0,1,\mathrm{}\mathrm{}.$$

Now each $`D_j`$ can be presented as

$$D_j=\underset{m_1+\mathrm{}m_L=j}{}𝒟(\stackrel{}{}_1^{(m_1)},\mathrm{}\stackrel{}{}_L^{(m_L)}),$$
(19)   
where $`𝒟(\stackrel{}{}_1^{(m_1)},\mathrm{}\stackrel{}{}_L^{(m_L)})`$ is the determinant of the $`L\times L`$ matrix consisting of the rows $`\stackrel{}{}_l^{(m_l)}`$ and the sum is taken over all such combinations of $`(m_1,\mathrm{}m_L)`$ that $`_{l=1}^Lm_l=j`$. $`_{ll^{^{}}}^{(0)}`$ can be readily expressed via derivatives of the conserved quantities with respect to the soliton parameters:

$$_{ll^{^{}}}^{(0)}=\frac{Q_l}{\kappa _l^{^{}}},$$
(20)   
and practical calculation of $`_{ll^{^{}}}^{(m)}`$ for $`m>0`$ can be simplified: $`_{ll^{^{}}}^{(m)}=\stackrel{}{W}_l|\widehat{}^{(12m)}\stackrel{}{U}_l^{^{}}.`$ Note, that in most of the cases solitary solution itself can be found only numerically using any of the well established methods for solving the nonlinear ode’s. Recurrent calculations of $`\widehat{}^{(12m)}\stackrel{}{U}_l`$ can be readily reduced to the numerically even simpler problem of solving of the linear inhomogeneous ode’s.

Because $`|\lambda |`$ was assumed to be small, Eq. (18) has an asymptotic character. Therefore to make it work some additional assumptions must be made about orders of $`D_j`$. If these assumptions are satisfied then Eq. (18) describes correctly the soliton spectrum and predicts bifurcations of the soliton. The corresponding eigenvalues can be found using Eq. (18) with any degree of accuracy. For example, let us assume that $`D_0ϵ^2`$ and $`D_{j>0}O(1)`$. Then, presenting $`\lambda ^2`$ as

$$\lambda ^2=ϵ^2\underset{j=0}{\overset{\mathrm{}}{}}\zeta _j,\zeta _jϵ^{2j},$$
(21)   
in the first order Eq. (18) gives a linear equation for $`\zeta _0`$,

$$D_0ϵ^2\zeta _0D_1=0,$$
(22)   
which indicates a threshold of the stationary bifurcation at $`D_0=0`$. This is precisely the condition $`det(J_{ij})=0`$ discussed in the introduction. Continuing to the next order one obtains

$$\lambda ^2=\frac{D_0}{D_1}\left(1\frac{D_0D_2}{D_1^2}+O(ϵ^4)\right).$$
(23)   
If $`D_1ϵ^2`$ then the asymptotic expression (23) fails. To have a balanced equation for $`\zeta _0`$, we must now assume that $`D_0ϵ^4`$. However, in such a case the Eq. 22 for $`\zeta _0`$ changes from linear to quadratic:

$$D_0ϵ^2\zeta _0D_1+ϵ^4\zeta _0^2D_2=0.$$
(24)   
Eq. (24) gives two threshold conditions $`D_0=0`$ and $`D_1^2=4D_0D_2`$, see Fig. 1. The latter condition indicates onset of the oscillatory instability for

$$D_1^2<4D_0D_2.$$
(25)   
Thus we have formulated analytic criterion for the oscillatory instability. It is also clear that the point $`D_{0,1}=0`$ is a source for the novel stationary instability, see rightmost region $`D_1^2>4D_0D_2`$, $`D_1>0`$ in Fig. 1, where an eigenvalue which is positive throughout this region can not be predicted by Eq. (23).

It follows by recurrence that if $`D_{j^{^{}}>0}ϵ^2`$ then to have a balanced equation for $`\zeta _0`$ we must assume that $`D_{j<j^{^{}}}ϵ^{2(1+j^{^{}})}`$. In other words asymptotic expansion near the neutral modes can only describe the soliton spectrum in regions of the parameter space which are close to codimension-$`(j^{^{}}+1)`$ bifurcation. If $`j^{^{}}=0`$ then only one condition must be satisfied and our asymptotic approach predicts presence of either two purely imaginary or two purely real eigenvalues, which can collide at zero. If $`j^{^{}}=1`$ then two conditions must be satisfied and the asymptotic approach predicts presence of two pairs of eigenvalues which can be real, imaginary or complex. In this situation the soliton becomes oscillation unstable providing that two pairs of imaginary eigenvalues collided. For each further $`j^{^{}}`$ two new eigenvalues come into play.

## IV Discussion

General formulae (18),(19) giving soliton eigenvalues with any degree of accuracy and criterion for the oscillatory instability (25) are main novel results of this work. At the same time expressions for the eigenvalues near the stationary instability threshold, analogs of the formula $`\lambda ^2=D_0/D_1+\mathrm{}`$, have been earlier obtained in a number of papers. It is instructive now to give explicit expressions for $`D_j`$ in the two simplest situations of one- and two-parameter solitons and to compare them with previously reported results. For the one parameter solitons: $`D_0=_{\kappa _1}Q_1`$, $`D_1=\stackrel{}{W}_1|\widehat{}^1\stackrel{}{U}_1`$, $`D_2=\stackrel{}{W}_1|\widehat{}^3\stackrel{}{U}_1`$. Using these formulae one can show that in the case when $`D_1O(1)`$ the first term in Eq. (23) gives the same expression for $`\lambda ^2`$ which was obtained in Refs. , where generalised Nonlinear Shrödinger equation  and equations describing propagation in quadratically nonlinear media  have been investigated. If $`D_1D_2>0`$ then it can be concluded that the second term in Eq. (23) indicates saturation of the growth rate when the distance from the instability threshold, $`D_0=0`$, growthes, which agrees with numerical results .

For the two-parameter solitons:

$`D_0=\left|\begin{array}{cc}\frac{Q_1}{\kappa _1}& \frac{Q_1}{\kappa _2}\\ \frac{Q_2}{\kappa _1}& \frac{Q_2}{\kappa _2}\end{array}\right|,`$ (28)   
$`D_1=\left|\begin{array}{cc}\frac{Q_1}{\kappa _1}& \frac{Q_1}{\kappa _2}\\ _{21}^{(1)}& _{22}^{(1)}\end{array}\right|+\left|\begin{array}{cc}_{11}^{(1)}& _{12}^{(1)}\\ \frac{Q_2}{\kappa _1}& \frac{Q_2}{\kappa _2}\end{array}\right|,`$ (33)   
$`D_2=\left|\begin{array}{cc}_{11}^{(1)}& _{12}^{(1)}\\ _{21}^{(1)}& _{22}^{(1)}\end{array}\right|+\left|\begin{array}{cc}\frac{Q_1}{\kappa _1}& \frac{Q_1}{\kappa _2}\\ _{21}^{(2)}& _{22}^{(2)}\end{array}\right|+\left|\begin{array}{cc}_{11}^{(2)}& _{12}^{(2)}\\ \frac{Q_2}{\kappa _1}& \frac{Q_2}{\kappa _2}\end{array}\right|.`$ (40)   
The threshold condition $`D_0=0`$ has been previously found for two-parameter solitons in different physical contexts . However derivation of an accurate expression for the soliton eigenvalues near this threshold has remained a controvertial problem. Indeed, comparing eigenvalues given by Eqs. (23), (28), (33),(40) and eigenvalues which can be calculated from the ordinary differential equations for soliton parameters presented in  one will discover that results are slightly different . It has also been argued  that the sign of $`(_{11}^{(1)}_{22}^{(1)}_{12}^{(1)}_{21}^{(1)})`$, which is the first term in Eq. (40), plays an important role in stability of two-parameter solitons. However Eqs. (23),(28),(33),(40) apparently conflict with this finding.

Among open problems I would like to mention derivation of finite-dimensional normal forms describing dynamical evolution of the soliton parameters near the oscillatory instability threshold. A guideline for this work can be theory of G. Iooss  for the normal forms of the reversible ordinary differential equations  in vicinity of the codimension-2 bifurcation, wich is an equivalent of the our point $`D_0=D_1=0`$. The simplest case of the codimension-1 stationary instability, $`D_0=0`$, has only one homoclinic orbit separating regions of the periodic oscillations from the spreading or collapse . The vicinity of the codimension-2 point can contain the very reach dynamics, including multiple homoclinic orbits and stochastic regimes.

## V Summary

General form of the asymptotic approach to stability problem of multi-parameter solitons in Hamiltonian systems has been developed. It has been shown that the asymptotic study of the soliton stability reduces to the calculation of a certain sequence of determinants, where the famous determinant of the matrix consisting from the derivatives of the system invariants with respect to the soliton parameters  is just the first in the series. Knowledge of these determinants allows to calculate eigenvalues governing soliton instability with arbitrary accuracy. The most important consequence is that the presented approach gives first analytic criterion for the oscillatory instability of solitons in Hamiltonian systems.

## VI Acknowledgment

Author acknowledges useful discussions with W.J. Firth and D.E. Pelinovsky and financial support from the Royal Society of Edinburgh and British Petroleum.