format: text

text:
# I Introduction

## I Introduction

Soliton solutions to nonlinear partial differential equations (PDEs) in (1+1) dimensions have been of great interest ever since the very discovery of completely integrable PDEs. There is a wide range of literature concerning integrable nonlinear PDEs and their soliton solutions (see, for instance, Refs.  and the references therein). Much is known about the behavior of solitons and their interactions in various integrable systems (soliton scattering, breather solutions, soliton bound states, etc). Such knowledge is very valuable not only for the underlying integrable systems, but also for nearly integrable systems which can be studied analytically by soliton perturbation theories.

It is an astonishing fact, however, that notwithstanding the almost three decades of advances in the study of soliton dynamics, still there are substantial gaps in our knowledge of soliton solutions in integrable nonlinear PDEs. Indeed, the reader familiar with the basics of the inverse scattering transform method knows that it is poles of the reflection coefficient (or, in modern terms, zeros of the Riemann-Hilbert problem) that give rise to the soliton solutions. The soliton solutions are usually derived by using one of the several well-known techniques, such as the dressing method  or the Riemann-Hilbert problem approach . However, in most publications (except, to our knowledge, Refs. , which are discussed below) only soliton solutions from simple poles are considered. It is usually assumed that a multiple-pole solution can be obtained in a straightforward way by coalescing several distinct poles (see, for instance, Ref. ) which describe multisoliton solutions. This would be indeed the case if such coalescing were a regular limit. However, this limit is obviously a singular one. Indeed, the soliton dressing matrix corresponding to a multisoliton solution is a rational matrix function which has distinct simple poles, while the coalescing procedure must produce multiple poles. Obviously, a more careful examination of this issue is necessary. This is the main subject of the present paper.

Soliton solutions corresponding to multiple poles, i.e., the higher-order solitons, have been investigated in the literature before. A soliton solution to the nonlinear Schrödinger (NLS) equation corresponding to a double pole was first given in Ref.  but without much analysis. The double- and triple-pole soliton solutions to the KdV equation were examined in Ref.  and the general $`N`$-pole soliton solution to the sine-Gordon equation was extensively studied in Ref.  using the associated Gelfand-Levitan-Marchenko equation. In Refs. , higher-order soliton solutions to the NLS equation were studied by employing the dressing method. Finally, in , higher order solitons in the Kadomtsev-Petviashvili I equation were derived by the inverse scattering method.

In this article, we consider the higher-order zeros of the Riemann-Hilbert problem for the $`N`$-wave system and study the corresponding soliton matrices. We derive the soliton dressing matrices for the simplest class of higher-order zeros – the elementary higher-order zeros. We call the $`n`$-th order zero $`k_1`$ of the Riemann-Hilbert problem elementary if the soliton matrix evaluated at $`k_1`$ has only one vector in the kernel, i.e. the geometric multiplicity of the zero is 1 (see Definition 1 in Sec. 3 for more details). The corresponding higher-order soliton solutions to the $`N`$-wave system are determined. Then we apply our theory to the physically important three-wave interaction model and derive the simplest higher-order soliton solution. In the generic case, this solution describes the breakup of a higher-order pumping wave into two higher-order elementary waves, and the reverse process. In non-generic cases, this solution could describe ($`i`$) the merger of a pumping sech wave and an elementary sech wave into two elementary waves (one sech and the other one higher-order); ($`ii`$) the breakup of a higher-order pumping wave into three sech waves — one pumping wave and two elementary waves. The higher-order soliton solution could also reproduce fundamental soliton solutions as a special case.

In general, one needs to consider zeros with the geometric multiplicities taking values from 1 to $`N1`$, where $`N`$ is the matrix dimension of the Riemann-Hilbert problem. The present work is the first step towards the solution of this general case. The point is that the soliton matrices derived here for the elementary zeros provide the building blocks for the most general case. We plan to address the general problem in the next paper. Thus, for the 3-wave interaction model, the higher-order soliton solutions we derived in this paper (which correspond to elementary higher-order zeros) may not be the most general higher-order soliton solutions in this wave system.

It is noted that the soliton dressing matrices for the higher-order soliton solutions were already a subject of interest in Refs. , where an ansatz for the soliton matrices was proposed for the $`2\times 2`$ Zakharov-Shabat spectral problem and expressions for the higher-order soliton solutions of the nonlinear Schrödinger equation were obtained. In the present paper, we study higher-order solitons in the $`N`$-wave system. Even though the evolution equations considered in  and this article are different, their Riemann-Hilbert formulations have a lot in common. In this paper, we place the emphasis on the derivation of the soliton matrices for the higher-order zeros and fill some gaps in the approach of Ref. . We show that the ansatz of  is precisely the soliton matrix for an elementary higher-order zero in a $`N\times N`$ spectral problem. For the $`2\times 2`$ Zakharov-Shabat spectral problem, any higher-order zero is elementary. Thus the ansatz of  is the general soliton matrix of the Zakharov-Shabat problem. Consequently, higher-order soliton solutions obtained in  are the general higher-order soliton solutions in the nonlinear Schrödinger equation. But in a $`N\times N`$ spectral problem with $`N>2`$, a higher-order zero is non-elementary in general. In that case, the ansatz of  will not be the general form of soliton matrices. In this paper, we also uncover the invariance property of the soliton matrix for a higher-order zero, which is necessary for the vector-parametrization of the soliton matrix to be self-consistent. The invariance property obtained in  is shown to be just a special case. Lastly, we point out that our derivation of higher-order solitons in the $`N`$-wave system is made for general dispersion laws and arbitrary matrix dimensions. In addition, our results can be generalized to the most general case of non-elementary zeros of the $`N\times N`$ spectral problem.

The paper is organized as follows. A summary on the Riemann-Hilbert problem is placed in section II. Section III is the central section of the paper. There we present the theory of soliton matrices corresponding to the higher-order zeros of the Riemann-Hilbert problem. For an elementary higher-order zero of the $`N\times N`$ Riemann-Hilbert problem (see definition in the text), the general soliton matrix is derived. This soliton matrix is similar to the ansatz as proposed in . In section IV, the theory is applied to derive the simplest higher-order soliton solution in the three-wave interaction model.

## II The Riemann-Hilbert problem approach: summary

The integrable nonlinear PDEs in 1+1 dimensions are associated with the matrix Riemann-Hilbert problem (consult, for instance, Refs. ). The matrix Riemann-Hilbert problem (below we work in the space of $`N\times N`$ matrices) is the problem of finding the holomorphic factorization, denoted below by $`\mathrm{\Phi }_+(k)`$ and $`\mathrm{\Phi }_{}^1(k)`$, in the complex plane of a nondegenerate matrix function $`G(k)`$ given on an oriented curve $`\gamma `$:

$`\mathrm{\Phi }_{}^1(k,x,t)\mathrm{\Phi }_+(k,x,t)`$ $`=`$ $`G(k,x,t)`$ (1)
$``$ $`\mathrm{exp}[\mathrm{\Lambda }(k)x\mathrm{\Omega }(k)t]G(k,0,0)\mathrm{exp}[\mathrm{\Lambda }(k)x+\mathrm{\Omega }(k)t],k\gamma .`$ (2)
Here the matrix functions $`\mathrm{\Phi }_+(k)`$ and $`\mathrm{\Phi }_{}^1(k)`$ are holomorphic in the two complementary domains of the complex $`k`$-plane: $`C_+`$ to the left and $`C_{}`$ to the right from the curve $`\gamma `$, respectively. The matrices $`\mathrm{\Lambda }(k)`$ and $`\mathrm{\Omega }(k)`$ are called the dispersion laws. In this paper, we require the dispersion laws to be diagonal (we have accounted for this on the right-hand side of equation (2) by writing the explicit $`(x,t)`$-dependence in the form of an exponent). The Riemann-Hilbert problem requires an appropriate normalization condition. Usually the curve $`\gamma `$ contains the infinite point $`k=\mathrm{}`$ of the complex plane and the normalization condition is formulated as

$$\mathrm{\Phi }_\pm (k,x,t)I,\text{as}k\mathrm{}.$$
(3)
This normalization condition is called the canonical normalization. Setting the normalization condition to an arbitrary nondegenerate matrix function $`S(x,t)`$ leads to the gauge equivalent integrable nonlinear PDE, e.g., the Landau-Lifshitz equation in the case of the NLS equation . Obviously, the new solution $`\widehat{\mathrm{\Phi }}_\pm (k,x,t)`$ to the Riemann-Hilbert problem, normalized to $`S(x,t)`$, is related to the canonical solution by the following transformation

$$\widehat{\mathrm{\Phi }}_\pm (k,x,t)=S(x,t)\mathrm{\Phi }(k,x,t).$$
(4)
Thus, without any loss of generality, we confine ourselves to the Riemann-Hilbert problem under the canonical normalization.

For physically applicable nonlinear PDEs the Riemann-Hilbert problem possesses the involution properties, which reduce the number of the dependent variables (complex fields). The $`N`$-wave interaction model admits the following involution property of the associated Riemann-Hilbert problem

$$\mathrm{\Phi }_+^{}(k)=\mathrm{\Phi }_{}^1(\overline{k}),\overline{k}k^{}.$$
(5)
Here the superscript “$``$” represents the Hermitian conjugate, and “\*” the complex conjugate. However, our approach can be trivially extended to the general case without such involution. To keep our treatment general, we will use the overlined quantities where applicable. The reduction to the involution is then done by associating the overline with the Hermitian conjugation in the case of vectors and matrices and with the complex conjugation in the case of scalar quantities.

To solve the Cauchy problem for the integrable nonlinear PDE posed on the whole axis $`x`$, one usually constructs the associated Riemann-Hilbert problem starting with the linear spectral equation

$$_x\mathrm{\Phi }(k,x,t)=\mathrm{\Phi }(k,x,t)\mathrm{\Lambda }(k)+U(k,x,t)\mathrm{\Phi }(k,x,t),$$
(6)
whereas the $`t`$-dependence is given by a similar equation

$$_t\mathrm{\Phi }(k,x,t)=\mathrm{\Phi }(k,x,t)\mathrm{\Omega }(k)+V(k,x,t)\mathrm{\Phi }(k,x,t).$$
(7)
The nonlinear integrable PDE corresponds to the compatibility condition of the system (6) and (7):

$$_tU_xV+[U,V]=0.$$
(8)
The essence of the approach based on the Riemann-Hilbert problem lies in the fact that the evolution governed by the complicated nonlinear PDE (8) is mapped to the evolution of the spectral data given by simpler equations such as (2) and (II). For details, consult Refs. .

Let the evolution equations for the spectral data be given. In our case, these are equation (2) for $`G`$ and equation (II) (see below) for the discrete data. Then the matrices $`U(k,x,t)`$ and $`V(k,x,t)`$ describing the evolution of $`\mathrm{\Phi }_\pm `$ can be retrieved from the Riemann-Hilbert problem. In our case, the potentials $`U(k,x,t)`$ and $`V(k,x,t)`$ are completely determined by the (diagonal) dispersion laws $`\mathrm{\Lambda }(k)`$ and $`\mathrm{\Omega }(k)`$ and the Riemann-Hilbert solution $`\mathrm{\Phi }\mathrm{\Phi }_\pm (k,x,t)`$. Indeed, let us assume that the dispersion laws are polynomial functions, i.e.,

$$\mathrm{\Lambda }(k)=\underset{j=0}{\overset{J_1}{}}A_jk^j,\mathrm{\Omega }(k)=\underset{j=0}{\overset{J_2}{}}B_jk^j.$$
(9)
Then using similar arguments as in Ref.  we get:

$$U=𝒫\{\mathrm{\Phi }\mathrm{\Lambda }\mathrm{\Phi }^1\},V=𝒫\{\mathrm{\Phi }\mathrm{\Omega }\mathrm{\Phi }^1\}.$$
(10)
Here the matrix function $`\mathrm{\Phi }(k)`$ is expanded into the asymptotic series,

$`\mathrm{\Phi }(k)=I+k^1\mathrm{\Phi }^{(1)}+k^2\mathrm{\Phi }^{(2)}+\mathrm{},k\mathrm{},`$

and the operator $`𝒫`$ cuts out the polynomial asymptotics of its argument as $`k\mathrm{}`$. An important property of matrices $`U`$ and $`V`$ is that

$$\text{Tr}U(k,x,t)=\text{Tr}\mathrm{\Lambda }(k),\text{Tr}V(k,x,t)=\text{Tr}\mathrm{\Omega }(k),$$
(11)
which evidently follows from equation (10). Below, let us consider the three-wave interaction system as an example . Set $`N=3`$,

$$\mathrm{\Lambda }(k)=ikA,A=\left(\begin{array}{ccc}a_1& 0& 0\\ 0& a_2& 0\\ 0& 0& a_3\end{array}\right),\mathrm{\Omega }(k)=ikB,B=\left(\begin{array}{ccc}b_1& 0& 0\\ 0& b_2& 0\\ 0& 0& b_3\end{array}\right),$$
(12)
where $`a_j`$ and $`b_j`$ are real with the elements of $`A`$ being ordered: $`a_1>a_2>a_3`$. From equation (10) we get

$$U=\mathrm{\Lambda }(k)+i[A,\mathrm{\Phi }^{(1)}],V=\mathrm{\Omega }(k)+i[B,\mathrm{\Phi }^{(1)}].$$
(13)
Setting

$$u_1=\sqrt{a_1a_2}\mathrm{\Phi }_{12}^{(1)},u_2=\sqrt{a_2a_3}\mathrm{\Phi }_{23}^{(1)},u_3=\sqrt{a_1a_3}\mathrm{\Phi }_{13}^{(1)},$$
(14)
assuming the involution (5), and using equation (13) in (8) we get the three-wave system:

$`_tu_1+v_1_xu_1+i\epsilon \overline{u}_2u_3`$ $`=`$ $`0,`$ (16)
$`_tu_2+v_2_xu_2+i\epsilon \overline{u}_1u_3`$ $`=`$ $`0,`$ (17)
$`_tu_3+v_3_xu_3+i\epsilon u_1u_2`$ $`=`$ $`0.`$ (18)
Here

$$v_1=\frac{b_2b_1}{a_1a_2},v_2=\frac{b_3b_2}{a_2a_3},v_3=\frac{b_3b_1}{a_1a_3},$$
(19)
$$\epsilon =\frac{a_1b_2a_2b_1+a_2b_3a_3b_2+a_3b_1a_1b_3}{[(a_1a_2)(a_2a_3)(a_1a_3)]^{1/2}}.$$
(20)
The group velocities satisfy the following condition

$$\frac{v_2v_3}{v_1v_3}=\frac{a_1a_2}{a_2a_3}<0.$$
(21)
The three-wave system (14) can be interpreted physically. It describes the interaction of three wave packets with complex envelopes $`u_1`$, $`u_2`$ and $`u_3`$ in a medium with quadratic nonlinearity.

It is often desirable to relate the inverse-scattering parameters $`a_j`$ and $`b_j`$ $`(j=1,2,3)`$ to the physical parameters $`\epsilon `$ and $`v_j(j=1,2,3)`$. This relation can be easily found from (19) and (20) as

$$a_1a_2=\frac{\epsilon ^2}{(v_1v_2)(v_1v_3)},a_2a_3=\frac{\epsilon ^2}{(v_1v_2)(v_3v_2)}.$$
(22)
The other parameters $`a_1a_3`$ and $`b_j(j=1,2,3)`$ can be readily obtained from equations (22) and (19). Note that the inverse-scattering parameters are not uniquely determined. In fact, one of $`a_j`$ and one of $`b_j`$ $`(j=1,2,3)`$ are free parameters. It is an invariance in the inverse-scattering formulation of the 3-wave system and it does not affect the physical solution in any way.

In general, the Riemann-Hilbert problem (2)-(3) has multiple solutions. Different solutions are related to each other by the rational matrix functions $`\mathrm{\Gamma }(k)`$ (which also depend on the variables $`x`$ and $`t`$) :

$$\stackrel{~}{\mathrm{\Phi }}_\pm (k,x,t)=\mathrm{\Phi }_\pm (k,x,t)\mathrm{\Gamma }(k,x,t).$$
(23)
The rational matrix $`\mathrm{\Gamma }(k)`$ must satisfy the canonical normalization condition: $`\mathrm{\Gamma }(k)I`$ for $`k\mathrm{}`$ and must have poles only in $`C_{}`$ (the inverse function $`\mathrm{\Gamma }^1(k)`$ then has poles in $`C_+`$ only). Such a rational matrix $`\mathrm{\Gamma }(k)`$ will be called the soliton matrix below, since it gives the soliton part of the solution to the integrable nonlinear PDE.

To specify a unique solution to the Riemann-Hilbert problem the set of the Riemann-Hilbert data must be given. These data are also called the spectral data. The full set of the spectral data comprises the matrix $`G(k,x,t)`$ on the right-hand side of equation (2) and the appropriate discrete data related to the zeros of $`det\mathrm{\Phi }_+(k)`$ and $`det\mathrm{\Phi }_{}^1(k)`$. We will confine ourselves to the case of the Riemann-Hilbert problem with zero index, i.e., when $`det\mathrm{\Phi }_+(k)`$ and $`det\mathrm{\Phi }_{}^1(k)`$ have equal number of zeros (counting the multiplicity). For instance, in the case of involution (5) the Riemann-Hilbert problem has zero index because the zeros appear in complex conjugate pairs: $`\overline{k}_j=k_j^{}`$. It is known  (see also Ref. ) that in the generic case the spectral data include simple (distinct) zeros $`k_1,\mathrm{},k_n`$ of $`det\mathrm{\Phi }_+(k)`$ and $`\overline{k}_1,\mathrm{},\overline{k}_n`$ of $`det\mathrm{\Phi }_{}^1(k)`$, in their holomorphicity domains, and the null vectors $`|v_1,\mathrm{},|v_n`$ and $`\overline{v}_1|,\mathrm{},\overline{v}_n|`$ from the respective kernels:

$$\mathrm{\Phi }_+(k_j)|v_j=0,\overline{v}_j|\mathrm{\Phi }_{}^1(\overline{k}_j)=0.$$
(24)
Using the property (11) one can verify that the zeros do not depend on the variables $`x`$ and $`t`$. The $`(x,t)`$-dependence of the null vectors can be easily derived by differentiation of (24) and use of the linear spectral equations (6)-(7). This dependence reads:

$`|v_j`$ $`=`$ $`\mathrm{exp}\{\mathrm{\Lambda }(k_j)x\mathrm{\Omega }(k_j)t\}|v_j^{(0)},`$ (26)
$`\overline{v}_j|`$ $`=`$ $`\overline{v}_j^{(0)}|\mathrm{exp}\{\mathrm{\Lambda }(\overline{k}_j)x+\mathrm{\Omega }(\overline{k}_j)t\},`$ (27)
where $`|v_j^{(0)}`$ and $`\overline{v}_j^{(0)}|`$ are some constant vectors.

The vectors in equation (II) together with the zeros constitute the full set of the discrete data necessary to specify the soliton matrix $`\mathrm{\Gamma }(k,x,t)`$ and, hence, unique solution to the Riemann-Hilbert problem (2)-(3). Indeed, by constructing the soliton matrix $`\mathrm{\Gamma }(k)`$ such that the following matrix functions

$$\varphi _+(k)=\mathrm{\Phi }_+(k)\mathrm{\Gamma }^1(k),\varphi _{}^1(k)=\mathrm{\Gamma }(k)\mathrm{\Phi }_{}^1(k)$$
(28)
are nondegenerate and holomorphic in the domains $`C_+`$ and $`C_{}`$, respectively, we reduce the Riemann-Hilbert problem with zeros to another one without zeros and hence uniquely solvable (for details see, for instance, Refs. ). Below by matrix $`\mathrm{\Gamma }(k)`$ we will imply the matrix from equation (28) which reduces the Riemann-Hilbert problem (2)-(3) to the one without zeros. The corresponding solution to the integrable PDE (8) is obtained by using the asymptotic expansion of the matrix $`\mathrm{\Phi }(k)`$ as $`k\mathrm{}`$ in the linear equation (6). In the $`N`$-wave interaction model it is given by formula (13). The pure soliton solutions are obtained by using the rational matrix $`\mathrm{\Phi }=\mathrm{\Gamma }(k)`$.

## III Soliton matrices for multiple zeros

In this section we consider the soliton solution corresponding to a single multiple zero of arbitrary order in the case of an arbitrary matrix dimension $`N`$. Such soliton solutions will be referred to as the higher-order solitons. We will derive the general formulae for the soliton matrices corresponding to an elementary higher-order zero (see the definition below) starting from the usual elementary soliton matrices of the Riemann-Hilbert problem. Our formulae for the soliton matrices corresponding to an elementary higher-order zero are similar to the previously proposed ansatz for the $`2\times 2`$ Zakharov-Shabat spectral problem . However, in our approach some essential invariance properties and simple evolution formulae for the vector parameters in the soliton matrices are given, which were not known before. Thus we simplify the ansatz of Refs.  and put it on the rigorous footing. Although we work in the case of involution (5), usual for applications in nonlinear physics, our approach is valid for the general Riemann-Hilbert problem with zero index. Moreover, we present our formulae in a form transferable without any changes to that general case.

Let $`\mathrm{\Phi }_+(k)`$ and $`\mathrm{\Phi }_{}^1(k)`$ from (2) each have but one zero of order $`n`$, $`k_1`$ and $`\overline{k}_1`$, respectively:

$$det\mathrm{\Phi }_+(k)=(kk_1)^n\phi (k),det\mathrm{\Phi }_{}^1(k)=(k\overline{k}_1)^n\overline{\phi }(k),$$
(29)
where $`det\phi (k_1)0`$ and $`det\overline{\phi }(\overline{k}_1)0`$. The geometric multiplicity of $`k_1`$ ($`\overline{k}_1`$) is defined as the number of the null vectors in the kernel of $`\mathrm{\Phi }_+(k_1)`$ ($`\mathrm{\Phi }_{}^1(\overline{k}_1)`$), see (24). It can be easily shown that the order of a zero is always greater or equal to its geometric multiplicity. It is also obvious that the geometric multiplicity of a zero is less than the matrix dimension. Before we proceed with the construction of the soliton matrix $`\mathrm{\Gamma }(k)`$ corresponding to the multiple zero of order $`n`$, two important properties must be pointed out. It is convenient to formulate them in the form of two lemmas.

###### Lemma 1

Suppose vectors $`|v_j(1jm)`$ are in the kernel of matrix $`\mathrm{\Phi }_+(k_1)`$, i.e.,

$$\mathrm{\Phi }_+(k_1)|v_j=0,j=1,\mathrm{},m,$$
(30)
where $`m`$ is less or equal to $`k_1`$’s geometric multiplicity. Define the new matrix $`\stackrel{~}{\mathrm{\Phi }}_+(k)\mathrm{\Phi }(k)_+\chi ^1(k)`$ where

$$\chi (k)=I\frac{k_1\overline{k}_1}{k\overline{k}_1}P,$$
(31)
$$P=\underset{i,j=1}{\overset{m}{}}|v_i(K^1)_{ij}\overline{v}_j|,K_{ij}=\overline{v}_i|v_j,$$
(32)
and vectors $`\overline{v}_j|(1jm)`$ are arbitrary but they make matrix $`K`$ invertible. Then matrix $`\stackrel{~}{\mathrm{\Phi }}_+(k)`$ is also holomorphic in the upper half plane. In addition, if a new vector $`|w`$ is in the kernel of $`\stackrel{~}{\mathrm{\Phi }}_+(k_1)`$ and is orthogonal to $`\overline{v}_j|(1jm)`$, i.e.,

$$\stackrel{~}{\mathrm{\Phi }}_+(k_1)|w=0,\overline{v}_j|w=0,j=1,\mathrm{},m,$$
(33)
then

$$\mathrm{\Phi }_+(k_1)|w=0,$$
(34)
i.e., $`|w`$ is also in the kernel of $`\mathrm{\Phi }_+(k_1)`$. Furthermore, $`|w`$ is linearly independent of $`|v_j(1jm)`$. Similar results exist for matrix $`\mathrm{\Phi }_{}^1(k)`$ where the multiplication is one the left.

Remark: it is easy to see that in order for $`K`$ to be invertible, it is necessary that vectors $`\overline{v}_j|(1jm)`$ be linearly independent. But this condition is not sufficient. However, if $`\overline{v}_j|=|v_j^{}(j=1,\mathrm{},m)`$, then it can be shown that $`K`$ is invertible.

Proof. The matrix $`P`$ is clearly a projector matrix, thus

$$\chi ^1(k)=I+\frac{k_1\overline{k}_1}{kk_1}P.$$
(35)
Then, expanding the holomorphic function $`\mathrm{\Phi }^+(k)`$ into the Taylor series and recalling equation (30), we see that

$`\mathrm{\Phi }^+(k)\chi ^1(k)=`$ $`\left\{\mathrm{\Phi }^+(k_1)+(kk_1){\displaystyle \frac{d\mathrm{\Phi }^+(k_1)}{dk}}+(kk_1)^2{\displaystyle \frac{d^2\mathrm{\Phi }^+(k_1)}{2!dk^2}}+\mathrm{}\right\}(1+{\displaystyle \frac{k_1\overline{k}_1}{kk_1}}P)`$ (36)
$`=`$ $`\mathrm{\Phi }^+(k_1)+(k_1\overline{k}_1){\displaystyle \frac{d\mathrm{\Phi }^+(k_1)}{dk}}P+(k_1\overline{k}_1)(kk_1){\displaystyle \frac{d^2\mathrm{\Phi }^+(k_1)}{2!dk^2}}P+\mathrm{},`$ (37)
which is clearly holomorphic.

Next, if

$`\stackrel{~}{\mathrm{\Phi }}_+(k_1)|w=0,`$

recalling the definition of $`\stackrel{~}{\mathrm{\Phi }}_+(k)`$ and expanding $`\mathrm{\Phi }^+(k)`$ into the Taylor series, we get

$$\mathrm{\Phi }^+(k_1)|w+(k_1\overline{k}_1)\frac{d\mathrm{\Phi }^+(k_1)}{dk}P|w=0.$$
(38)
Since by assumption,

$`\overline{v}_j|w=0,j=1,\mathrm{},m,`$

thus,

$`P|w=0,`$

consequently,

$$\mathrm{\Phi }^+(k_1)|w=0.$$
(39)
Lastly, $`|w`$ is linearly independent of $`|v_j(j=1,\mathrm{},m)`$ because the matrix $`K`$ is invertible. Q.E.D.

Corollary 1 Suppose the kernel of $`\mathrm{\Phi }_+(k_1)`$ is spanned by vectors $`|v_j(1jm)`$ where $`m`$ is the geometric multiplicity of zero $`k_1`$. Define matrices $`\chi (k)`$, $`\stackrel{~}{\mathrm{\Phi }}_+(k)`$ and projector $`P`$ as in Lemma 1. Then, there exists no vector in the kernel of $`\stackrel{~}{\mathrm{\Phi }}_+(k_1)`$ which is simultaneously orthogonal to $`\overline{v}_j|(1jm)`$.

###### Lemma 2

Suppose that $`\mathrm{\Phi }_+(k_1)`$ has $`r`$ independent vectors in the kernel:

$$\mathrm{\Phi }_+(k_1)|v_j=0,j=1,\mathrm{},r,$$
(40)
i.e., $`\mathrm{rank}\mathrm{\Phi }_+(k_1)=Nr`$. Then the following matrix function $`\stackrel{~}{\mathrm{\Phi }}_+(k)\mathrm{\Phi }_+(k)\chi ^1(k)`$, where matrix $`\chi (k)`$ is as defined in equation (31) but with

$$P=\underset{i,j=1}{\overset{r}{}}|v_i(K^1)_{ij}\overline{v}_j|,K_{ij}=\overline{v}_i|v_j,$$
(41)
has at most $`r`$ vectors in the kernel at $`k=k_1`$, i.e., $`\mathrm{rank}\stackrel{~}{\mathrm{\Phi }}_+(k_1)Nr`$. Here vectors $`\overline{v}_j|(1jr)`$ are arbitrary but they make matrix $`K`$ invertible.

Proof. This lemma is easy to prove by contradiction. Suppose that there are at least $`r+1`$ independent vectors $`|u_1`$, …, $`|u_{r+1}`$ in the kernel of $`\stackrel{~}{\mathrm{\Phi }}_+(k_1)`$ defined above. Then one can find a non-zero vector $`|X`$ in the kernel of $`\stackrel{~}{\mathrm{\Phi }}_+(k_1)`$ such that

$$\overline{v}_j|X=0,j=1,\mathrm{},r.$$
(42)
Indeed, substitution of the expansion

$`|X={\displaystyle \underset{j=1}{\overset{r+1}{}}}C_j|u_j`$

into Eq. (42) leads to an underdetermined, hence, solvable system of equations

$`{\displaystyle \underset{j=1}{\overset{r+1}{}}}\overline{v}_i|u_jC_j=0,i=1,\mathrm{},r`$

which have non-zero solutions. But then, according to the second part of Lemma 1, $`|X`$ is also in the kernel of $`\mathrm{\Phi }_+(k_1)`$, thus

$$|X=\underset{j=1}{\overset{r}{}}C_j|v_j.$$
(43)
Substituting Eq. (43) into (42) and recalling that the matrix $`K`$ is invertible, we find that $`C_j=0,j=1,\mathrm{},r`$, hence $`X=0`$. Thus we have arrived at a contradiction. Q.E.D. (Note that a similar lemma is valid for $`\mathrm{\Phi }_{}^1(k)`$ at $`k=\overline{k}_1`$ with the multiplication on the left.)

To clarify the implications of Lemma 2 for the soliton matrix $`\mathrm{\Gamma }(k)`$ of the higher-order zeros, $`k=k_1`$ of $`det\mathrm{\Phi }_+(k)`$ and $`k=\overline{k}_1`$ of $`det\mathrm{\Phi }_{}^1(k)`$, let us examine the way such matrix is constructed. Starting from the solution $`\mathrm{\Phi }_\pm (k)`$ to the Riemann-Hilbert problem (2)-(3), one looks for the independent null vectors for the matrices $`\mathrm{\Phi }_+(k_1)`$ and $`\mathrm{\Phi }_{}^1(\overline{k}_1)`$:

$$\mathrm{\Phi }_+(k_1)|v_{i1}=0,\overline{v}_{i1}|\mathrm{\Phi }_{}^1(\overline{k}_1)=0,i=1,\mathrm{},s_1,$$
(44)
where $`s_1`$ is the smaller of $`k_1`$ and $`\overline{k}_1`$’s geometric multiplicities. Here we allow the two geometric multiplicities to be different in general, but they are always the same in the case of involution (5). Next, one constructs the elementary matrix

$$\chi _1(k)=I\frac{k_1\overline{k}_1}{k\overline{k}_1}P_1,$$
(45)
where

$$P_1=\underset{i,j}{\overset{s_1}{}}|v_{i1}(K^1)_{ij}\overline{v}_{j1}|,K_{ij}=\overline{v}_{i1}|v_{j1}.$$
(46)
It can be shown that $`\text{det}\chi _1=\left(\frac{kk_1}{k\overline{k}_1}\right)^{s_1}`$. If $`s_1<n`$, where $`n`$ is the order of the two zeros, then one considers the matrix functions $`\stackrel{~}{\mathrm{\Phi }}_+(k)=\mathrm{\Phi }_+(k)\chi _1^1(k)`$ and $`\stackrel{~}{\mathrm{\Phi }}_{}^1(k)=\chi _1(k)\mathrm{\Phi }_{}^1(k)`$. From Lemma 1, we know that matrices $`\stackrel{~}{\mathrm{\Phi }}_+(k)`$ and $`\stackrel{~}{\mathrm{\Phi }}_{}^1(k)`$ are also holomorphic in the respective half planes of the complex plane. In addition, $`k_1`$ ($`\overline{k}_1`$) is still a zero of $`\text{det}\stackrel{~}{\mathrm{\Phi }}_+(k)`$ ($`\text{det}\stackrel{~}{\mathrm{\Phi }}_{}^1(k)`$). Repeating the above steps one gets the elementary matrices $`\chi _1(k)`$, …, $`\chi _r(k)`$ such that $`s_1+s_2+\mathrm{}+s_r=n`$. Therefore,

$$\mathrm{\Gamma }(k)=\chi _r(k)\mathrm{}\chi _2(k)\chi _1(k),$$
(47)
where matrices $`\chi _l(k)`$ and projectors $`P_l`$ are as defined in equations (45) and (46) but the independent vectors $`|v_{il}`$ and $`\overline{v}_{il}|`$ ($`i=1,\mathrm{},s_l`$) are from the kernels of $`(\mathrm{\Phi }_+\chi _1^1\mathrm{}\chi _{l1}^1)(k_1)`$ and $`(\chi _{l1}\mathrm{}\chi _1\mathrm{\Phi }_{}^1)(\overline{k}_1)`$ respectively.

Lemma 2 indicates that in fact the sequence of ranks of the projectors $`P_l`$ in the matrix $`\mathrm{\Gamma }(k)`$ given by equation (47), i.e. built in the described way, is non-increasing:

$$\text{rank }P_r\text{rank }P_{r1}\mathrm{}\text{rank }P_1.$$
(48)
This result allows one to classify possible occurrences of a higher-order zero of the Riemann-Hilbert problem for arbitrary matrix dimension $`N`$. In general, for zeros of the same order $`n`$, different sequences of ranks in formula (48) give different classes of the higher-order soliton solutions. In the present paper we consider in detail only the higher-order zeros when the sequence of ranks (48) is the simplest possible: $`\text{rank }P_l=1`$, $`l=1,\mathrm{},n`$. We introduce the following definition.

Definition 1. In the soliton matrix (47) corresponding to a higher-order zero $`k_1`$ of a Riemann-Hilbert problem, if the ranks of all projectors $`P_l(1ln)`$ are 1, then we call this zero an elementary higher-order zero.

Remark 1: We observe from equation (48) that a higher-order zero (of arbitrary algebraic multiplicity) is elementary if and only if $`\text{rank}P_1=1`$, i.e., the geometric multiplicity of the zero is 1.

Remark 2: If the matrix dimension $`N=2`$ (as for the nonlinear Schrödinger equation), then all higher-order zeros are elementary since $`\text{rank}P_1`$ is always equal to 1.

Below we derive the soliton matrix $`\mathrm{\Gamma }(k)`$ and its inverse for an elementary higher-order zero. The results are presented in the following lemma.

###### Lemma 3

Consider a pair of elementary higher-order zeros of order $`n`$: $`k=k_1`$ in $`C_+`$ and $`k=\overline{k}_1`$ in $`C_{}`$. Then the corresponding soliton matrix $`\mathrm{\Gamma }(k)`$ and its inverse can be cast in the following form:

$$\mathrm{\Gamma }(k)=I+\underset{l=1}{\overset{n}{}}\underset{j=1}{\overset{l}{}}\frac{|\overline{q}_j\overline{p}_{l+1j}|}{(k\overline{k}_1)^{n+1l}}=I+(|\overline{q}_n,\mathrm{},|\overline{q}_1)\overline{D}(k)\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right),$$
(50)
$$\mathrm{\Gamma }^1(k)=I+\underset{l=1}{\overset{n}{}}\underset{j=1}{\overset{l}{}}\frac{|p_{l+1j}q_j|}{(kk_1)^{n+1l}}=I+(|p_1,\mathrm{},|p_n)D(k)\left(\begin{array}{c}q_n|\\ \mathrm{}\\ q_1|\end{array}\right),$$
(51)
where the matrices $`D(k)`$ and $`\overline{D}(k)`$ are defined as

$$\overline{D}(k)=\left(\begin{array}{cccc}\frac{1}{(k\overline{k}_1)}& 0& \mathrm{}& 0\\ \frac{1}{(k\overline{k}_1)^2}& \frac{1}{(k\overline{k}_1)}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& 0\\ \frac{1}{(k\overline{k}_1)^n}& \mathrm{}& \frac{1}{(k\overline{k}_1)^2}& \frac{1}{(k\overline{k}_1)}\end{array}\right),D(k)=\left(\begin{array}{cccc}\frac{1}{(kk_1)}& \frac{1}{(kk_1)^2}& \mathrm{}& \frac{1}{(kk_1)^n}\\ 0& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \frac{1}{(kk_1)}& \frac{1}{(kk_1)^2}\\ 0& \mathrm{}& 0& \frac{1}{(kk_1)}\end{array}\right),$$
(52)
and vectors $`|p_j,\overline{p}_j|,q_j|,|\overline{q}_j(j=1,\mathrm{},n)`$ are independent of $`k`$.

Remark: In , the ansatz of the form (3) was proposed for higher-order solitons in the nonlinear Schrödinger equation. The lemma above, together with Remark 2 below Definition 1, shows that their ansatz is in fact the most general soliton matrix for the nonlinear Schrödinger equation. If $`N>2`$, their ansatz then is just the soliton matrix for elementary higher-order zeros.

Proof. The representation (3) can be proved by induction. Consider, for instance, formula (50). Obviously, this formula is valid for $`n=1`$. In this case, $`\mathrm{\Gamma }(k)`$ reduces to an elementary matrix $`\chi (k)`$. Now, suppose that formula (50) is valid for $`n=m`$. Then we need to show that it is valid for $`n=m+1`$ as well. Indeed, denote the soliton matrices for $`n=m`$ and $`n=m+1`$ by $`\mathrm{\Gamma }(k)`$ and $`\stackrel{~}{\mathrm{\Gamma }}(k)`$ respectively. Then taking into account expression (47) and recalling our assumption of the elementary higher-order zero, we have

$$\stackrel{~}{\mathrm{\Gamma }}(k)=\chi _{m+1}(k)\mathrm{\Gamma }(k)=\left(I+\frac{|v_{m+1}\overline{v}_{m+1}|}{k\overline{k}_1}\right)\left(I+\underset{l=1}{\overset{m}{}}\underset{j=1}{\overset{l}{}}\frac{|\overline{q}_j\overline{p}_{l+1j}|}{(k\overline{k}_1)^{m+1l}}\right).$$
(53)
Here, for simplicity of the formulae below, we have normalized the vectors $`|v_{m+1}`$ and $`\overline{v}_{m+1}|`$ such that $`\overline{v}_{m+1}|v_{m+1}=\overline{k}_1k_1`$. Let us now multiply the two terms in the right-hand side of equation (53) and compute the coefficients at the poles:

$$\stackrel{~}{\mathrm{\Gamma }}(k)=I+\frac{\stackrel{~}{A}_1}{k\overline{k}_1}+\frac{\stackrel{~}{A}_2}{(k\overline{k}_1)^2}+\mathrm{}+\frac{\stackrel{~}{A}_{m+1}}{(k\overline{k}_1)^{m+1}},$$
(54)
where

$`\stackrel{~}{A}_{m+1}`$ $`=`$ $`|v_{m+1}\overline{v}_{m+1}|A_m=|v_{m+1}\overline{v}_{m+1}|\overline{q}_1\overline{p}_1|,`$
$`\stackrel{~}{A}_m`$ $`=`$ $`|v_{m+1}\overline{v}_{m+1}|A_{m1}+A_m`$
$`=`$ $`|v_{m+1}\overline{v}_{m+1}|\left(|\overline{q}_2\overline{p}_1|+|\overline{q}_1\overline{p}_2|\right)+|\overline{q}_1\overline{p}_1|`$
$`=`$ $`\left(|v_{m+1}\overline{v}_{m+1}|\overline{q}_2+|\overline{q}_1\right)\overline{p}_1|+|v_{m+1}\overline{v}_{m+1}|\overline{q}_1\overline{p}_2|,`$
$`\stackrel{~}{A}_{m1}`$ $`=`$ $`|v_{m+1}\overline{v}_{m+1}|A_{m2}+A_{m1}`$
$`=`$ $`|v_{m+1}\overline{v}_{m+1}|\left(|\overline{q}_3\overline{p}_1|+|\overline{q}_2\overline{p}_2|+|\overline{q}_1\overline{p}_3|\right)+|\overline{q}_2\overline{p}_1|+|\overline{q}_1\overline{p}_2|`$
$`=`$ $`\left(|v_{m+1}\overline{v}_{m+1}|\overline{q}_3+|\overline{q}_2\right)\overline{p}_1|+\left(|v_{m+1}\overline{v}_{m+1}|\overline{q}_2+|\overline{q}_1\right)\overline{p}_2|`$
$`+|v_{m+1}\overline{v}_{m+1}|\overline{q}_1\overline{p}_3|,`$
$`\mathrm{},`$
$`\stackrel{~}{A}_1`$ $`=`$ $`|v_{m+1}\overline{v}_{m+1}|+{\displaystyle \underset{j=1}{\overset{m}{}}}|\overline{q}_{m+1j}\overline{p}_j|.`$
Define new vectors:

$`|\stackrel{~}{\overline{q}}_1=|v_{m+1}\overline{v}_{m+1}|\overline{q}_1,|\stackrel{~}{\overline{q}}_j=|v_{m+1}\overline{v}_{m+1}|\overline{q}_j+|\overline{q}_{j1},j=2,\mathrm{},m,`$

$$|\stackrel{~}{\overline{q}}_{m+1}=|\overline{q}_m,\overline{p}_{m+1}|=\frac{\overline{v}_{m+1}|_{j=1}^{m1}\overline{v}_{m+1}|\overline{q}_{j+1}\overline{p}_{mj+1}|}{\overline{v}_{m+1}|\overline{q}_1}.$$
(55)
Then matrices $`\stackrel{~}{A}_1,\mathrm{},\stackrel{~}{A}_{m+1}`$ take the following representation:

$$\stackrel{~}{A}_{m+2l}=\underset{j=1}{\overset{l}{}}|\stackrel{~}{\overline{q}}_{l+1j}\overline{p}_j|,l=1,\mathrm{},m+1.$$
(56)
Thus formula (50) is valid for $`\stackrel{~}{\mathrm{\Gamma }}(k)`$ as well. It is noted that we must also show that the denominator in formula (55) is nonzero. This is easy to show, as $`\overline{v}_{m+1}|\overline{q}_1`$ is actually (up to a factor $`(\overline{k}_1k_1)^m`$) a product of inner products $`\overline{v}_{j+1}|v_j(1jm)`$, where $`|v_j`$ and $`\overline{v}_j|`$ are the projector vectors in matrix $`\chi _j`$:

$$\chi _j(k)=I\frac{k_1\overline{k}_1}{k\overline{k}_1}|v_j\overline{v}_j|$$
(57)
\[see equations (47) and (53)\]. If $`\overline{v}_{j+1}|v_j=0`$ for some $`j`$, then Lemma 1 indicates that the projector $`P_j`$ in matrix $`\chi _j`$ \[see (45) to (47)\] would have rank higher than 1, which contradicts our assumption of elementary higher-order zeros. Thus $`\overline{v}_{j+1}|v_j0`$ for all $`j`$, consequently, $`\overline{v}_{m+1}|\overline{q}_10`$. Expression for $`\mathrm{\Gamma }^1(k)`$ (51) can be proved in the same way. Q.E.D.

In the expressions for $`\mathrm{\Gamma }(k)`$ (50) and $`\mathrm{\Gamma }^1(k)`$ (51) there are twice as many vectors as in the elementary matrices (45) and (47). As the result, only half of the vector parameters, namely $`|p_1,\mathrm{},|p_n`$ and $`\overline{p}_1|,\mathrm{},\overline{p}_n|`$, are independent. To derive the formulae for the rest of the vector parameters in (3) we can use the identity $`\mathrm{\Gamma }(k)\mathrm{\Gamma }^1(k)=I`$. The poles of $`\mathrm{\Gamma }(k)\mathrm{\Gamma }^1(k)`$ at $`k=k_1`$, starting from the highest order pole, give:

$`\mathrm{\Gamma }(k_1)|p_1q_1|=0,`$
$`\mathrm{\Gamma }(k_1)\left(|p_2q_1|+|p_1q_2|\right)+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}\mathrm{\Gamma }(k_1)}{\text{d}k}}|p_1q_1|=0,`$
$`\mathrm{\Gamma }(k_1)\left(|p_3q_1|+|p_2q_2|+|p_1q_3|\right)+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}\mathrm{\Gamma }(k_1)}{\text{d}k}}\left(|p_2q_1|+|p_1q_2|\right)+{\displaystyle \frac{1}{2!}}{\displaystyle \frac{\text{d}^2\mathrm{\Gamma }(k_1)}{\text{d}k^2}}|p_1q_1|=0,`$
$`\mathrm{}.`$
Hence, we obtain:

$`\mathrm{\Gamma }(k_1)|p_1=0,`$ (59)
$`\mathrm{\Gamma }(k_1)|p_2+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}\mathrm{\Gamma }(k_1)}{\text{d}k}}|p_1=0,`$ (60)
$`\mathrm{\Gamma }(k_1)|p_3+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}\mathrm{\Gamma }(k_1)}{\text{d}k}}|p_2+{\displaystyle \frac{1}{2!}}{\displaystyle \frac{\text{d}^2\mathrm{\Gamma }(k_1)}{\text{d}k^2}}|p_1=0,`$ (61)
$`\mathrm{},`$ (62)
$`\mathrm{\Gamma }(k_1)|p_n+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}\mathrm{\Gamma }(k_1)}{\text{d}k}}|p_{n1}+\mathrm{}+{\displaystyle \frac{1}{(n1)!}}{\displaystyle \frac{\text{d}^{n1}\mathrm{\Gamma }(k_1)}{\text{d}k^{n1}}}|p_1=0.`$ (63)
Equations (III) can be written in a compact form for the following matrix $`𝚪(k)`$:

$$𝚪(k_1)\left(\begin{array}{c}|p_1\\ \mathrm{}\\ |p_n\end{array}\right)=0,𝚪(k)\left(\begin{array}{cccc}\mathrm{\Gamma }& 0& \mathrm{}& 0\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Gamma }& \mathrm{\Gamma }& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& 0\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Gamma }& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Gamma }& \mathrm{\Gamma }\end{array}\right).$$
(64)
Note that, as a block matrix, $`𝚪(k)`$ has (lower-triangular) Toeplitz form, i.e. along each diagonal it has the same (matrix) element.

In much the same way, by considering the poles at $`k=\overline{k}_1`$ in $`\mathrm{\Gamma }(k)\mathrm{\Gamma }^1(k)`$, one derives the following formula

$$(\overline{p}_1|,\mathrm{},\overline{p}_n|)\overline{𝚪}(\overline{k}_1)=0,\overline{𝚪}(k)=\left(\begin{array}{cccc}\mathrm{\Gamma }^1& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Gamma }^1& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Gamma }^1\\ 0& \mathrm{\Gamma }^1& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Gamma }^1\\ 0& \mathrm{}& 0& \mathrm{\Gamma }^1\end{array}\right).$$
(65)
Equations (64) and (65) allow us to find the expressions for the dependent vector parameters. For convenience of the presentation, let us introduce the following $`k`$-dependent vectors:

$$\overline{Z}_j(k)|=\underset{l=j}{\overset{n}{}}\frac{\overline{p}_{l+1j}|}{(k\overline{k}_1)^{n+1l}},|Z_j(k)=\underset{l=j}{\overset{n}{}}\frac{|p_{l+1j}}{(kk_1)^{n+1l}}.$$
(66)
Then, by reordering the summation in (3) we get

$$\mathrm{\Gamma }(k)=I+(|\overline{q}_n,\mathrm{},|\overline{q}_1)\left(\begin{array}{c}\overline{Z}_n(k)|\\ \mathrm{}\\ \overline{Z}_1(k)|\end{array}\right),$$
(68)
$$\mathrm{\Gamma }^1(k)=I+(|Z_n(k),\mathrm{},|Z_1(k))\left(\begin{array}{c}q_n|\\ \mathrm{}\\ q_1|\end{array}\right).$$
(69)
Let us now substitute the expression (68) into equation (64) and solve for $`|\overline{q}_1,\mathrm{},|\overline{q}_n`$. We have

$`|p_1+(|\overline{q}_n,\mathrm{},|\overline{q}_1)\left(\begin{array}{c}\overline{Z}_n(k_1)|p_1\\ \mathrm{}\\ \overline{Z}_1(k_1)|p_1\end{array}\right)=0,`$
$`|p_2+(|\overline{q}_n,\mathrm{},|\overline{q}_1)\left(\begin{array}{c}\overline{Z}_n(k_1)|p_2+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_n(k_1)|p_1\\ \mathrm{}\\ \overline{Z}_1(k_1)|p_2+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_1(k_1)|p_1\end{array}\right)=0,`$
$`|p_3+(|\overline{q}_n,\mathrm{},|\overline{q}_1)\left(\begin{array}{c}\overline{Z}_n(k_1)|p_3+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_n(k_1)|p_2+\frac{1}{2!}\frac{\text{d}^2}{\text{d}k^2}\overline{Z}_n(k_1)|p_1\\ \mathrm{}\\ \overline{Z}_1(k_1)|p_3+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_1(k_1)|p_2+\frac{1}{2!}\frac{\text{d}^2}{\text{d}k^2}\overline{Z}_1(k_1)|p_1\end{array}\right)=0,`$
$`\mathrm{}.`$
Hence

$$(|\overline{q}_n,\mathrm{},|\overline{q}_1)=(|p_1,\mathrm{},|p_n)\overline{𝒦}^1,$$
(73)
where

$$\overline{𝒦}=\left(\begin{array}{cccc}\overline{Z}_n(k_1)|p_1& \overline{Z}_n(k_1)|p_2+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_n(k_1)|p_1& \mathrm{}& \underset{l=1}{\overset{n}{}}\frac{1}{(nl)!}\frac{\text{d}^{nl}}{\text{d}k^{nl}}\overline{Z}_n(k_1)|p_l\\ \mathrm{}& \mathrm{}& & \mathrm{}\\ \overline{Z}_1(k_1)|p_1& \overline{Z}_1(k_1)|p_2+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{Z}_1(k_1)|p_1& \mathrm{}& \underset{l=1}{\overset{n}{}}\frac{1}{(nl)!}\frac{\text{d}^{nl}}{\text{d}k^{nl}}\overline{Z}_1(k_1)|p_l\end{array}\right).$$
(74)
Similarly, we get

$$\left(\begin{array}{c}q_n|\\ \mathrm{}\\ q_1|\end{array}\right)=𝒦^1\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right),$$
(75)
where

$$𝒦=\left(\begin{array}{ccc}\overline{p}_1|Z_n(\overline{k}_1)& \mathrm{}& \overline{p}_1|Z_1(\overline{k}_1)\\ \overline{p}_2|Z_n(\overline{k}_1)+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{p}_1|Z_n(\overline{k}_1)& \mathrm{}& \overline{p}_2|Z_1(\overline{k}_1)+\frac{1}{1!}\frac{\text{d}}{\text{d}k}\overline{p}_1|Z_1(\overline{k}_1)\\ \mathrm{}& & \mathrm{}\\ \underset{l=1}{\overset{n}{}}\frac{1}{(nl)!}\frac{\text{d}^{nl}}{\text{d}k^{nl}}\overline{p}_l|Z_n(\overline{k}_1)& \mathrm{}& \underset{l=1}{\overset{n}{}}\frac{1}{(nl)!}\frac{\text{d}^{nl}}{\text{d}k^{nl}}\overline{p}_l|Z_1(\overline{k}_1)\end{array}\right).$$
(76)
In terms of the independent vector parameters, the soliton matrices (50) and (51) can be rewritten as

$$\mathrm{\Gamma }(k)=I(|p_1,\mathrm{},|p_n)\overline{𝒦}{}_{}{}^{1}\overline{D}(k)\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right),$$
(77)
$$\mathrm{\Gamma }^1(k)=I(|p_1,\mathrm{},|p_n)D(k)𝒦^1\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right),$$
(78)
where matrices $`𝒦`$ and $`\overline{𝒦}`$ are given in equations (76) and (74).

The soliton matrices given by (77) and (78) possess invariance properties. The invariance is the transformation of the independent vector parameters which preserves the form of the soliton matrices and equations defining the vector parameters, i.e. equations (64)-(65). Let us first consider transformations of vectors $`|p_j(j=1,\mathrm{},n)`$. Suppose these vectors are transformed as

$$(|p_1,\mathrm{},|p_n)=(|\stackrel{~}{p}_1,\mathrm{},|\stackrel{~}{p}_n)B,$$
(79)
where $`B`$ is a $`k`$-independent matrix which, in general, depends on $`(x,t)`$. Here the vectors $`\overline{p}_j|(j=1,\mathrm{},n)`$ remain intact. Simple calculations show that the new vectors $`|\stackrel{~}{p}_1,\mathrm{},|\stackrel{~}{p}_n`$ satisfy equation (64) if and only if the matrix $`B`$ has upper-triangular Toeplitz form,

$$B=\left(\begin{array}{cccccc}b_1& b_2& \mathrm{}& \mathrm{}& b_n& \\ 0& b_1& b_2& \mathrm{}& \mathrm{}& \\ \mathrm{}& 0& \mathrm{}& \mathrm{}& \mathrm{}& \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& b_2& \\ 0& \mathrm{}& \mathrm{}& 0& b_1& \end{array}\right).$$
(80)
Further, we note that under the transformation (79)-(80) the matrix $`\overline{𝒦}`$ transforms as

$$\overline{𝒦}=\stackrel{~}{\overline{𝒦}}B,$$
(81)
where matrix $`\stackrel{~}{\overline{𝒦}}`$ is as given by equation (74) but with vectors $`|p_j`$ replaced by the new vectors $`|\stackrel{~}{p}_j`$. From formulae (79) and (81) it is seen that the form (77) of matrix $`\mathrm{\Gamma }(k)`$ is preserved. We still need to show that for matrix $`B`$ of the form (80), the transformation (79) also preserves the form (78) of matrix $`\mathrm{\Gamma }^1(k)`$. Notice that matrix $`D(k)`$ also has upper-triangular Toeplitz form, thus $`D(k)`$ and $`B`$ are commutable. Utilizing this property, we can easily show that under the transformation (79), matrix $`𝒦`$ transforms as

$$𝒦=\stackrel{~}{𝒦}B,$$
(82)
where $`\stackrel{~}{𝒦}`$ is given by equation (76) but with $`|p_j`$ replaced by $`|\stackrel{~}{p}_j`$. Thus the form of matrix $`\mathrm{\Gamma }^1(k)`$ is also preserved. In short, soliton matrices (77) and (78) are invariant under the transformation (79) with matrix $`B`$ given by (80).

Similarly, we can show that soliton matrices (77) and (78) are also invariant under the transformation

$$\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right)=\overline{B}\left(\begin{array}{c}\stackrel{~}{\overline{p}}_1|\\ \mathrm{}\\ \stackrel{~}{\overline{p}}_n|\end{array}\right),$$
(83)
where the $`k`$-independent matrix $`\overline{B}`$ (which, in general, depends on $`(x,t)`$) has lower-triangular Toeplitz form,

$$\overline{B}=\left(\begin{array}{ccccc}\overline{b}_1& 0& \mathrm{}& \mathrm{}& 0\\ \overline{b}_2& \overline{b}_1& 0& \mathrm{}& \mathrm{}\\ \mathrm{}& \overline{b}_2& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& 0\\ \overline{b}_n& \mathrm{}& \mathrm{}& \overline{b}_2& \overline{b}_1\end{array}\right),$$
(84)
and vectors $`|p_j(j=1,\mathrm{},n)`$ remain intact.

Summarizing, we conclude that soliton matrices (77) and (78) are invariant under the triangular Toeplitz transformations

$$(|p_1,\mathrm{},|p_n)=(|\stackrel{~}{p}_1,\mathrm{},|\stackrel{~}{p}_n)B,\left(\begin{array}{c}\overline{p}_1|\\ \mathrm{}\\ \overline{p}_n|\end{array}\right)=\overline{B}\left(\begin{array}{c}\stackrel{~}{\overline{p}}_1|\\ \mathrm{}\\ \stackrel{~}{\overline{p}}_n|\end{array}\right),$$
(85)
of the independent vectors $`|p_j`$ and $`\overline{p}_j|(1jn)`$. Here $`B`$ and $`\overline{B}`$ are arbitrary lower and upper triangular Toeplitz matrices, respectively, in general $`(x,t)`$-dependent. Here we point out that the invariance transformation found in  is given by $`b_j=\overline{b}_j=0(2jn1)`$, i.e., only $`b_1,b_n,\overline{b}_1`$ and $`\overline{b}_n`$ being non-zero. Thus it is just a special case of the invariance property of the soliton matrices.

The invariance transformations indicate that arbitrary sets of vectors $`|p_1,\mathrm{},|p_n`$ and $`\overline{p}_1|,\mathrm{},\overline{p}_n|`$ satisfying equations (64) and (65) can be chosen as the independent vector parameters. This is, in fact, also a necessary condition for such vector parameterization of the soliton matrix to be self-consistent.

Now let us derive the $`(x,t)`$-dependence of the vector parameters which enter the soliton matrix. We can start with the fact that the soliton matrix $`\mathrm{\Gamma }(k,x,t)`$ must satisfy equations (6)-(7) with some potentials $`U(k,x,t)`$ and $`V(k,x,t)`$:

$`_x\mathrm{\Gamma }(k,x,t)`$ $`=`$ $`\mathrm{\Gamma }(k,x,t)\mathrm{\Lambda }(k)+U(k,x,t)\mathrm{\Gamma }(k,x,t),`$ (87)
$`_t\mathrm{\Gamma }(k,x,t)`$ $`=`$ $`\mathrm{\Gamma }(k,x,t)\mathrm{\Omega }(k)+V(k,x,t)\mathrm{\Gamma }(k,x,t).`$ (88)
The derivation is based on the use of equations (64) and (65) (quite similar to the derivation of equations (II) in section II). First of all we need to find the equations for the triangular block Toeplitz matrices $`𝚪`$ and $`\overline{𝚪}`$. To this goal one needs to differentiate equations (3) with respect to $`k`$ up to the $`(n1)`$-th order. It is easy to see that, for instance, the equations for the $`𝚪`$ have the same form as equations (3):

$`_x𝚪(k,x,t)`$ $`=`$ $`𝚪(k,x,t)𝚲(k)+𝐔(k,x,t)𝚪(k,x,t),`$ (90)
$`_t𝚪(k,x,t)`$ $`=`$ $`𝚪(k,x,t)𝛀(k)+𝐕(k,x,t)𝚪(k,x,t),`$ (91)
if we introduce the lower-triangular block Toeplitz matrices $`𝚲`$, $`𝛀`$, $`𝐔`$, and $`𝐕`$:

$$𝚲\left(\begin{array}{cccc}\mathrm{\Lambda }& 0& \mathrm{}& 0\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Lambda }& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{\Lambda }& 0\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Lambda }& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Lambda }& \mathrm{\Lambda }\end{array}\right),𝛀\left(\begin{array}{cccc}\mathrm{\Omega }& 0& \mathrm{}& 0\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Omega }& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{\Omega }& 0\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Omega }& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Omega }& \mathrm{\Omega }\end{array}\right),$$
(92)
$$𝐔\left(\begin{array}{cccc}U& 0& \mathrm{}& 0\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}U& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& U& 0\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}U& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}U& U\end{array}\right),𝐕\left(\begin{array}{cccc}V& 0& \mathrm{}& 0\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}V& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& V& 0\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}V& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}V& V\end{array}\right).$$
(93)
Indeed, this is due to the fact that the matrix multiplication in (III) exactly reproduces the Leibniz rule for higher-order derivatives of a product. Similarly, using the equations for $`\mathrm{\Gamma }^1`$, one finds that

$`_x\overline{𝚪}(k,x,t)`$ $`=`$ $`\overline{𝚲}(k)\overline{𝚪}(k,x,t)\overline{𝚪}(k,x,t)\overline{𝐔}(k,x,t),`$ (95)
$`_t\overline{𝚪}(k,x,t)`$ $`=`$ $`\overline{𝛀}(k)\overline{𝚪}(k,x,t)\overline{𝚪}(k,x,t)\overline{𝐕}(k,x,t),`$ (96)
for the upper-triangular block Toeplitz matrices $`\overline{𝚲}`$, $`\overline{𝛀}`$, $`\overline{𝐔}`$, and $`\overline{𝐕}`$:

$$\overline{𝚲}=\left(\begin{array}{cccc}\mathrm{\Lambda }& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Lambda }& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Lambda }\\ 0& \mathrm{\Lambda }& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Lambda }\\ 0& \mathrm{}& 0& \mathrm{\Lambda }\end{array}\right),\overline{𝛀}=\left(\begin{array}{cccc}\mathrm{\Omega }& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Omega }& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{\Omega }\\ 0& \mathrm{\Omega }& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{\Omega }\\ 0& \mathrm{}& 0& \mathrm{\Omega }\end{array}\right),$$
(97)
$$\overline{𝐔}=\left(\begin{array}{cccc}U& \frac{1}{1!}\frac{\text{d}}{\text{d}k}U& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}U\\ 0& U& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}U\\ 0& \mathrm{}& 0& U\end{array}\right),\overline{𝐕}=\left(\begin{array}{cccc}V& \frac{1}{1!}\frac{\text{d}}{\text{d}k}V& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}V\\ 0& V& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}V\\ 0& \mathrm{}& 0& V\end{array}\right).$$
(98)
The $`(x,t)`$-dependence of the vector parameters $`|p_1,\mathrm{},|p_n`$ and $`\overline{p}_1|,\mathrm{},\overline{p}_n|`$ can be found by differentiation of equations (64) and (65) with the help of equations (III) and (93). First, we note that for commuting matrices the corresponding block Toeplitz matrices as introduced above also commute with each other. Second, it is shown in the Appendix that for a diagonal matrix \[e.g. $`\mathrm{\Lambda }(k)x+\mathrm{\Omega }(k)t`$\] the operation of raising to the exponent commutes with the construction of the block Toeplitz matrix. Therefore, taking into account the invariance property, we find the $`(x,t)`$-dependence of the vector parameters as

$$\left(\begin{array}{c}|p_1\\ \mathrm{}\\ |p_n\end{array}\right)=\mathrm{exp}\left\{𝚲(k_1)x𝛀(k_1)t\right\}\left(\begin{array}{c}|p_1^{(0)}\\ \mathrm{}\\ |p_n^{(0)}\end{array}\right),$$
(100)
$$(\overline{p}_1|,\mathrm{},\overline{p}_n|)=(\overline{p}_1^{(0)}|,\mathrm{},\overline{p}_n^{(0)}|)\mathrm{exp}\left\{\overline{𝚲}(\overline{k}_1)x+\overline{𝛀}(\overline{k}_1)t\right\}.$$
(101)
Here the superscript “0” is used to denote constant vectors and the exponents stand for the triangular block Toeplitz matrices:

$$\mathrm{exp}\left\{𝚲(k_1)x𝛀(k_1)t\right\}=\left(\begin{array}{cccc}E(k_1)& 0& \mathrm{}& \mathrm{\hspace{0.33em}0}\\ \frac{1}{1!}\frac{\text{d}}{\text{d}k}E(k_1)& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& E(k_1)& \mathrm{\hspace{0.33em}0}\\ \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}E(k_1)& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}E(k_1)& E(k_1)\end{array}\right),$$
(103)
$$\mathrm{exp}\left\{\overline{𝚲}(\overline{k}_1)x+\overline{𝛀}(\overline{k}_1)t\right\}=\left(\begin{array}{cccc}E^1(\overline{k}_1)& \frac{1}{1!}\frac{\text{d}}{\text{d}k}E^1(\overline{k}_1)& \mathrm{}& \frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}E^1(\overline{k}_1)\\ 0& E^1(\overline{k}_1)& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \frac{1}{1!}\frac{\text{d}}{\text{d}k}E^1(\overline{k}_1)\\ 0& \mathrm{}& 0& E^1(\overline{k}_1)\end{array}\right),$$
(104)
where $`E(k)\mathrm{exp}\left\{\mathrm{\Lambda }(k)x\mathrm{\Omega }(k)t\right\}`$. After the temporal and spatial evolutions for vectors $`|p_j`$ and $`\overline{p}_j|`$ have been obtained as above, the corresponding higher-order soliton solution can be obtained from equations (8), (13), (50) and (73).

## IV Application to the three wave interaction model

Here we apply the theory developed in the previous section to the three-wave interaction model (14). The three-wave model has wide applications in nonlinear physics. For instance, under the additional constraint $`u_j=iq_j`$ where $`q_j`$ are real variables, it describes the “exact resonance” in parametric interaction of three wave packets, while under the reduction of the dispersion laws (12): $`a_3=a_1`$, $`a_2=0`$, $`b_3=b_1`$, $`b_2=0`$ and the condition $`u_2=u_1`$, it models the generation of second harmonics. The usual (fundamental) soliton solutions to the three-wave interaction model have been well studied (consult Ref. ). Such solitons approach sech profiles as $`t\pm \mathrm{}`$ on the characteristics $`xv_jt=const`$.

Let us consider the simplest higher-order solitons in the three-wave system: solitons which correspond to an elementary higher-order zero of order 2. Here we should take into account the involution property given by equation (5). For instance, we have

$`\overline{k}=k^{},\overline{p}_j|=|p_j^{}`$

(here and below the overline is associated with the Hermitian or, in the case of scalar quantities, complex conjugation). Then the soliton matrix reads

$$\mathrm{\Gamma }(k)=I(|p_1,|p_2)\overline{𝒦}^1\left(\begin{array}{c}\overline{Z}_2(k)|\\ \overline{Z}_1(k)|\end{array}\right),$$
(105)
where

$$\overline{𝒦}=\left(\begin{array}{cc}\overline{Z}_2(k_1)|p_1& \overline{Z}_2(k_1)|p_2+\frac{\text{d}}{\text{d}k}\overline{Z}_2(k_1)|p_1\\ \overline{Z}_1(k_1)|p_1& \overline{Z}_1(k_1)|p_2+\frac{\text{d}}{\text{d}k}\overline{Z}_1(k_1)|p_1\end{array}\right),$$
(106)
and

$`\overline{Z}_2(k)|={\displaystyle \frac{\overline{p}_1|}{k\overline{k}_1}},\overline{Z}_1(k)|={\displaystyle \frac{\overline{p}_2|}{k\overline{k}_1}}+{\displaystyle \frac{\overline{p}_1|}{(k\overline{k}_1)^2}}.`$

The $`(x,t)`$-dependence of the vector parameters $`|p_1,|p_2`$ has the following form

$$\left(\begin{array}{c}|p_1\\ |p_2\end{array}\right)=\left(\begin{array}{cc}E(k_1)& 0\\ \frac{\text{d}}{\text{d}k}E(k_1)& E(k_1)\end{array}\right)\left(\begin{array}{c}|p_1^{(0)}\\ |p_2^{(0)}\end{array}\right),E(k_1)=e^{ik_1(Ax+Bt)}.$$
(107)
We denote $`k_1=\xi +i\eta `$, where $`\xi `$ and $`\eta `$ are real numbers ($`\eta >0`$ since $`k_1`$ lies the upper half plane of the complex plane), and choose the following parameterization of the constant vectors $`|p_1^{(0)}`$ and $`|p_2^{(0)}`$:

$$|p_1^{(0)}=2i\eta \left(\begin{array}{c}\theta _1^{(1)}\\ \theta _2^{(1)}\\ \theta _3^{(1)}\end{array}\right),|p_2^{(0)}=\left(\begin{array}{c}\theta _1^{(2)}\\ \theta _2^{(2)}\\ \theta _3^{(2)}\end{array}\right),$$
(108)
where $`\theta _j^{(i)}`$’s are complex constants. It is noted that due to the invariance property (85), where the matrix $`B`$ contains two arbitrary complex constants, we have 2 free components in each vector in formula (108). Hence there are 10 free real parameters (including $`\xi `$ and $`\eta `$) in the higher-order soliton solution.

The $`(x,t)`$-dependence of the components of the vector parameters reads

$$p_{1j}=2i\eta \theta _j^{(1)}e^{f_j/2i\chi _j},p_{2j}=\left[\theta _j^{(2)}+f_j\theta _j^{(1)}\right]e^{f_j/2i\chi _j},$$
(109)
where

$$f_j=2\eta (a_jx+b_jt),\chi _j=\xi (a_jx+b_jt),j=1,2,3.$$
(110)
By simple calculations we obtain the elements of matrix $`\overline{𝒦}`$ as

$$\overline{𝒦}_{11}=2i\eta \underset{j=1}{\overset{3}{}}|\theta _j^{(1)}|^2e^{f_j},\overline{𝒦}_{12}=\underset{j=1}{\overset{3}{}}\left(\overline{\theta }_j^{(1)}\theta _j^{(2)}+(f_j1)|\theta _j^{(1)}|^2\right)e^{f_j},$$
(111)
$$\overline{𝒦}_{21}=\underset{j=1}{\overset{3}{}}\left(\theta _j^{(1)}\overline{\theta }_j^{(2)}+(f_j1)|\theta _j^{(1)}|^2\right)e^{f_j},\overline{𝒦}_{22}=\frac{1}{2i\eta }\underset{j=1}{\overset{3}{}}\left(|\theta _j^{(2)}+(f_j1)\theta _j^{(1)}|^2+|\theta _j^{(1)}|^2\right)e^{f_j}.$$
(112)
It is easy to verify that the determinant of $`\overline{𝒦}`$ is

$$det\overline{𝒦}=\underset{i,j=1}{\overset{3}{}}\left(|\theta _i^{(1)}\theta _j^{(1)}|^2+\frac{1}{2}|\theta _i^{(2)}\theta _j^{(1)}\theta _i^{(1)}\theta _j^{(2)}+(f_if_j)\theta _i^{(1)}\theta _j^{(1)}|^2\right)e^{f_i+f_j},$$
(113)
which is always non-zero:

For the soliton solution corresponding to the matrix (105) we need the first-order term of its asymptotics as $`k\mathrm{}`$:

$$\mathrm{\Gamma }^{(1)}=\frac{1}{det\overline{𝒦}}\left(\overline{𝒦}_{22}|p_1\overline{p}_1|+\overline{𝒦}_{11}|p_2\overline{p}_2|\overline{𝒦}_{12}|p_1\overline{p}_2|\overline{𝒦}_{21}|p_2\overline{p}_1|\right).$$
(114)
Using formulae (109) for $`|p_1`$ and $`|p_2`$ and (111, 112) for the elements of $`\overline{𝒦}`$ we get

$$\mathrm{\Gamma }_{lm}^{(1)}=\frac{2i\eta }{det\overline{𝒦}}e^{(f_l+f_m)/2i(\chi _l\chi _m)}\underset{j=1}{\overset{3}{}}C_{lmj}e^{f_j},$$
(115)
where

$`C_{lmj}=\left[\theta _j^{(1)}\theta _l^{(1)}\left(f_jf_l1\right)+\theta _j^{(2)}\theta _l^{(1)}\theta _l^{(2)}\theta _j^{(1)}\right]\left[\overline{\theta }_m^{(1)}\overline{\theta }_j^{(1)}\left(f_mf_j+1\right)+\overline{\theta }_m^{(2)}\overline{\theta }_j^{(1)}\overline{\theta }_j^{(2)}\overline{\theta }_m^{(1)}\right]`$

$$\theta _l^{(1)}\overline{\theta }_m^{(1)}|\theta _j^{(1)}|^2.$$
(116)
The three nonlinear waves $`u_1,u_2,u_3`$ are given by formula (14). Thus

$$u_1=\sqrt{a_1a_2}\mathrm{\Gamma }_{12}^{(1)},u_2=\sqrt{a_2a_3}\mathrm{\Gamma }_{23}^{(1)},u_3=\sqrt{a_1a_3}\mathrm{\Gamma }_{13}^{(1)},$$
(117)
where $`\mathrm{\Gamma }_{ij}^{(1)}`$ are given by Eqs. (115) and (116). To be explicit, our soliton solution corresponding to an elementary higher-order zero of order 2 in the three-wave interaction system is

$$u_1=\frac{2i\eta \sqrt{a_1a_2}}{det\overline{𝒦}}e^{(f_1+f_2)/2i(\chi _1\chi _2)}\underset{j=1}{\overset{3}{}}C_{12j}e^{f_j},$$
(118)
$$u_2=\frac{2i\eta \sqrt{a_2a_3}}{det\overline{𝒦}}e^{(f_2+f_3)/2i(\chi _2\chi _3)}\underset{j=1}{\overset{3}{}}C_{23j}e^{f_j},$$
(119)
$$u_3=\frac{2i\eta \sqrt{a_1a_3}}{det\overline{𝒦}}e^{(f_1+f_3)/2i(\chi _1\chi _3)}\underset{j=1}{\overset{3}{}}C_{13j}e^{f_j},$$
(120)
where $`det\overline{𝒦},C_{ijk},f_k`$ and $`\chi _k`$ are given by Eqs. (110), (113) and (116).

The above solutions are fairly complicated. But some information about them can be gained from considering the asymptotics as $`t\pm \mathrm{}`$. Evidently, the $`t`$-asymptotics is nonzero only on the characteristics:

$`z_1(f_1f_2)/2=\eta (a_1a_2)(xv_1t),`$

$`z_2(f_2f_3)/2=\eta (a_2a_3)(xv_2t),`$

$`z_3(f_1f_3)/2=\eta (a_1a_3)(xv_3t).`$

The asymptotic formulae depend on the relation between the velocities of the waves. For definiteness, let us choose $`v_2<v_1`$. This is equivalent to the condition $`\epsilon >0`$ in view that

$`\epsilon =\left({\displaystyle \frac{(a_1a_2)(a_2a_3)}{a_1a_3}}\right)^{1/2}(v_1v_2),`$

as follows from formulae (19)-(20). Then the condition (21) requires that $`v_3`$ lies between $`v_2`$ and $`v_1`$:

$$v_2<v_3<v_1.$$
(121)
Further, we notice that any solution $`\stackrel{~}{u}_1,\stackrel{~}{u}_2,\stackrel{~}{u}_3`$ of the three-wave interaction model (14) in the case of the opposite inequality, i.e. $`v_1<v_3<v_2`$, is mapped onto the solution satisfying the inequality (121) by the following transformation: $`\stackrel{~}{u}_j(x,t;v_1,v_2,v_3)=u_j(x,t;v_1,v_2,v_3)`$. Thus, the case of $`v_1<v_2`$ is easy to recover (it describes the reverse process to that of $`v_2<v_1`$).

The asymptotic formulae also depend on whether some of the components in vectors $`\theta ^{(1)}`$ and $`\theta ^{(2)}`$ are zero or not. We first consider the generic case when none of the parameters $`\theta _j^{(1)}`$ for $`j=1,2,3`$ is zero. Define the following real quantities:

$`\alpha _{lm}=\mathrm{ln}\left({\displaystyle \frac{|\theta _m^{(1)}|}{|\theta _l^{(1)}|}}\right),\varrho _{lm}+i\sigma _{lm}={\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\theta _l^{(2)}}{\theta _l^{(1)}}}{\displaystyle \frac{\theta _m^{(2)}}{\theta _m^{(1)}}}\right)+\alpha _{lm},`$ (122)
$`\phi _j^{(s)}=\text{arg}(\theta _j^{(s)}),\phi _{lm}=\phi _l^{(1)}\phi _m^{(1)}{\displaystyle \frac{\xi \alpha _{lm}}{\eta }},`$ (123)
and denote

$$z_{12}=z_1\alpha _{12},z_{23}=z_2\alpha _{23},z_{13}=z_3\alpha _{13}.$$
(124)
Then, simple calculations show that the asymptotics of the waves (118), (119) and (120) are as follows:

$$u_10,u_20,t\mathrm{};u_30,t\mathrm{},$$
(125)
$$u_12i\eta \sqrt{a_1a_2}\frac{(z_{12}+\varrho _{12})\mathrm{sinh}z_{12}(1+i\sigma _{12})\mathrm{cosh}z_{12}}{\mathrm{cosh}^2z_{12}+(z_{12}+\varrho _{12})^2+\sigma _{12}^2}e^{i(\phi _{12}\xi z_{12}/\eta )},t\mathrm{},$$
(126)
$$u_22i\eta \sqrt{a_2a_3}\frac{(z_{23}+\varrho _{23})\mathrm{sinh}z_{23}(1+i\sigma _{23})\mathrm{cosh}z_{23}}{\mathrm{cosh}^2z_{23}+(z_{23}+\varrho _{23})^2+\sigma _{23}^2}e^{i(\phi _{23}\xi z_{23}/\eta )},t\mathrm{},$$
(127)
$$u_32i\eta \sqrt{a_1a_3}\frac{(z_{13}+\varrho _{13})\mathrm{sinh}z_{13}(1+i\sigma _{13})\mathrm{cosh}z_{13}}{\mathrm{cosh}^2z_{13}+(z_{13}+\varrho _{13})^2+\sigma _{13}^2}e^{i(\phi _{13}\xi z_{13}/\eta )},t\mathrm{}.$$
(128)
We see that as $`t\mathrm{}`$, only the pumping wave $`u_3`$ is non-zero, while as $`t\mathrm{}`$ only the elementary waves $`u_1`$ and $`u_2`$ are nonzero. Thus in the generic case of higher-order solitons under condition (121), the solution describes the breakdown of the pumping higher-order soliton $`u_3`$ into the higher-order solitons of the elementary waves $`u_1`$ and $`u_2`$. \[For the opposite inequalities in formula (121), the solution describes the reverse process: merger of the two elementary waves $`u_1`$ and $`u_2`$ into the pumping wave $`u_3`$.\] These properties are identical to fundamental solitons (see, for instance, pp. 174-184 in Ref. ). However, differences between higher-order solitons and fundamental solitons are also obvious: none of the asymptotics (126) to (128) of the higher-order solitons is sech-shaped, while the asymptotics of fundamental solitons are all sech-shaped.

Asymptotics (125) to (128) are invalid in the non-generic cases when at least one of the parameters $`\theta _j^{(1)}`$, $`j=1,2,3`$, is zero. Consider first the case of $`\theta _1^{(1)}=0`$ and $`\theta _2^{(1)}`$ and $`\theta _3^{(1)}`$ being non-zero. There are two possibilities depending on whether $`\theta _1^{(2)}`$ is zero or not.

(a) If $`\theta _1^{(2)}0`$, then the asymptotics of the waves (118), (119) and (120) become

$$u_10,t\mathrm{};u_30,t\mathrm{};$$
(129)
$$u_1i\eta \sqrt{a_1a_2}e^{i(\phi _1^{(2)}\phi _2^{(1)}\xi z_1/\eta )}\text{sech}(z_1\beta _1),t\mathrm{},$$
(130)
$$u_2i\eta \sqrt{a_2a_3}e^{i(\phi _2^{(1)}\phi _3^{(1)}\xi z_2/\eta )}\text{sech}(z_2\beta _2),t\mathrm{},$$
(131)
$$u_22i\eta \sqrt{a_2a_3}\frac{(z_{23}+\varrho _{23})\mathrm{sinh}z_{23}(1+i\sigma _{23})\mathrm{cosh}z_{23}}{\mathrm{cosh}^2z_{23}+(z_{23}+\varrho _{23})^2+\sigma _{23}^2}e^{i(\phi _{23}\xi z_{23}/\eta )},t\mathrm{},$$
(132)
$$u_3i\eta \sqrt{a_1a_3}e^{i(\phi _1^{(2)}\phi _3^{(1)}\xi z_3/\eta )}\text{sech}(z_3\beta _3),t\mathrm{},$$
(133)
where parameters $`\beta _j(j=1,2,3)`$ are defined as

$$\beta _1=\mathrm{ln}\left(\frac{|\theta _2^{(1)}|}{|\theta _1^{(2)}|}\right),\beta _2=\mathrm{ln}\left(\frac{|\theta _3^{(1)}|}{|\theta _2^{(1)}|}\right),\beta _3=\mathrm{ln}\left(\frac{|\theta _3^{(1)}|}{|\theta _1^{(2)}|}\right).$$
(134)
The above asymptotics have two important features. One is that as $`t\mathrm{}`$, both $`u_3`$ and $`u_2`$ waves are non-zero. Thus, the higher-order soliton solution (118) to (120) in this non-generic case does not describe the breakdown of the pumping wave $`u_3`$. Instead, it describes a new process:

$$u_2+u_3u_1+u_2.$$
(135)
This is very different from fundamental solitons. The other feature is that the waves $`u_2`$ and $`u_3`$ as $`t\mathrm{}`$ and the wave $`u_1`$ as $`t\mathrm{}`$ all have sech profiles, but the wave $`u_2`$ as $`t\mathrm{}`$ is the higher-order soliton.

(b) If $`\theta _1^{(2)}=0`$, then the higher-order solitons (118) to (120) are degenerate:

$$u_1(x,t)=u_3(x,t)=0u_2(x,t)=u_2^{(0)}(xv_2t),$$
(136)
where $`u_2^{(0)}(x)`$ is the initial solution of $`u_2`$. This is a trivial solution.

If two components of the vector $`\theta ^{(1)}`$ are zero, then higher-order soliton solutions (118) to (120) reduce to fundamental-soliton solutions or trivial solutions. For instance, if $`\theta _1^{(1)}=\theta _2^{(1)}=0`$, $`\theta _3^{(1)}0`$, $`\theta _1^{(2)}0`$, and $`\theta _2^{(2)}0`$, then the asymptotics of the waves become

$$u_10,u_20,t\mathrm{};u_30,t\mathrm{};$$
(137)
$$u_1i\eta \sqrt{a_1a_2}e^{i(\phi _1^{(2)}\phi _2^{(2)}\xi z_1/\eta )}\text{sech}(z_1\stackrel{~}{\beta }_1),t\mathrm{},$$
(138)
$$u_2i\eta \sqrt{a_2a_3}e^{i(\phi _2^{(2)}\phi _3^{(1)}\xi z_2/\eta )}\text{sech}(z_2\stackrel{~}{\beta }_2),t\mathrm{},$$
(139)
$$u_3i\eta \sqrt{a_1a_3}e^{i(\phi _1^{(2)}\phi _3^{(1)}\xi z_3/\eta )}\text{sech}(z_3\beta _3),t\mathrm{},$$
(140)
where $`\beta _3`$ is defined in equation (134) and

$$\stackrel{~}{\beta }_1=\mathrm{ln}\left(\frac{|\theta _2^{(2)}|}{|\theta _1^{(2)}|}\right),\stackrel{~}{\beta }_2=\mathrm{ln}\left(\frac{|\theta _3^{(1)}|}{|\theta _2^{(2)}|}\right).$$
(141)
This is the fundamental soliton solution. If one or both of $`\theta _1^{(2)}`$ and $`\theta _2^{(2)}`$ is zero, the solution is trivial (similar to (136)). We note that when $`\theta _1^{(1)}=\theta _2^{(1)}=0`$, then $`\theta _3^{(1)}`$ can not be zero, because otherwise, the denominator $`det\overline{𝒦}`$ in the solution is zero.

It turns out that consideration of the case when $`\theta _3^{(1)}=0`$ is similar to the above case of $`\theta _1^{(1)}=0`$ with the only difference that now the elementary waves $`u_1`$ and $`u_2`$ are interchanged. For instance, when $`\theta _3^{(1)}=0`$, $`\theta _1^{(1)}0`$, $`\theta _2^{(1)}0`$, and $`\theta _3^{(2)}0`$ we have the following (also new) process:

$$u_1+u_3u_1+u_2,$$
(142)
where the waves $`u_1`$ and $`u_3`$ as $`t\mathrm{}`$ and the wave $`u_2`$ as $`t\mathrm{}`$ all have sech profiles, while the wave $`u_1`$ as $`t\mathrm{}`$ is the higher-order soliton \[see equation (127)\].

The only (different) case which is left to consider is the case of $`\theta _2^{(1)}=0`$ with $`\theta _1^{(1)}0`$ and $`\theta _3^{(1)}0`$. The asymptotics depends on whether $`\theta _2^{(2)}`$ is zero or not. In the former case, i.e. $`\theta _2^{(2)}=0`$, we have a degenerate solution,

$`u_1(x,t)=u_2(x,t)=0,u_3(x,t)=u_3^{(0)}(xv_3t),`$

which is similar to solution (136). If however $`\theta _2^{(2)}0`$, then the asymptotics of the waves (118)-(120) are as follows:

$$u_10,u_20,t\mathrm{},$$
(143)
$$u_32i\eta \sqrt{a_1a_3}\frac{(z_{13}+\varrho _{13})\mathrm{sinh}z_{13}(1+i\sigma _{13})\mathrm{cosh}z_{13}}{\mathrm{cosh}^2z_{13}+(z_{13}+\varrho _{13})^2+\sigma _{13}^2}e^{i(\phi _{13}\xi z_{13}/\eta )},t\mathrm{};$$
(144)
$$u_1i\eta \sqrt{a_1a_2}e^{i(\phi _1^{(1)}\phi _2^{(2)}\xi z_1/\eta )}\text{sech}(z_1\widehat{\beta }_1),t\mathrm{},$$
(145)
$$u_2i\eta \sqrt{a_2a_3}e^{i(\phi _2^{(2)}\phi _3^{(1)}\xi z_2/\eta )}\text{sech}(z_2\stackrel{~}{\beta }_2),t\mathrm{},$$
(146)
$$u_3i\eta \sqrt{a_1a_3}e^{i(\phi _1^{(1)}\phi _3^{(1)}\xi z_3/\eta )}\text{sech}(z_3\widehat{\beta }_3),t\mathrm{},$$
(147)
where $`\stackrel{~}{\beta }_2`$ is defined in equation (141) and

$$\widehat{\beta }_1=\mathrm{ln}\left(\frac{|\theta _2^{(2)}|}{|\theta _1^{(1)}|}\right),\widehat{\beta }_3=\mathrm{ln}\left(\frac{|\theta _3^{(1)}|}{|\theta _1^{(1)}|}\right).$$
(148)
These asymptotic formulae describe yet another new process:

$$u_3u_1+u_2+u_3,$$
(149)
where waves $`u_1`$, $`u_2`$, and $`u_3`$ as $`t\mathrm{}`$ all have sech profiles, while the pumping wave $`u_3`$ as $`t\mathrm{}`$ is more complicated. Thus, this process describes a breakup of the higher-order pumping wave into three sech waves, the two elementary waves and the pumping wave.

Lastly, we present the graphical pictures of the above higher-order solitons for both the generic and non-generic cases. In all figures, the common solution parameters are $`(a_1,a_2,a_3)=(2,1,1)`$, $`(b_1,b_2,b_3)=(0.5,2,1),`$ $`\xi =1,\eta =1`$, and $`\theta ^{(2)}=(1,1+i,2)`$. Only the vector $`\theta ^{(1)}`$ is different. It is easy to check that for these parameters, the inequality (121) holds, thus the asymptotics of these higher-order solitons have been described in the previous text. In all figures, the solid lines are $`|u_1|`$, the dashed lines are $`|u_2|`$, and the dashed-dotted lines are $`|u_3|`$.

First, we illustrate the generic solution with $`\theta ^{(1)}=(1,i,1)`$ in Fig. 1. As we can see from this figure as well as the asymptotics (125) to (128), as $`t\mathrm{}`$, only the pumping $`u_3`$ solution is non-zero. As $`t\mathrm{}`$, this $`u_3`$ wave breaks into elementary $`u_1`$ and $`u_2`$ waves. This process is similar to fundamental solitons. But there is a difference: the asymptotics of each wave in Fig. 1 has a complex structure which signals that it is a higher-order soliton instead of a fundamental soliton. Next, we let $`\theta _1^{(1)}`$ approach zero. Specifically, we let $`\theta _1^{(1)}=10^4`$, while the $`\theta _2^{(1)}`$ and $`\theta _3^{(1)}`$ values remain the same. The corresponding soliton solution is illustrated in Fig. 2. We see that in this case, the pumping $`u_3`$ wave at $`t\mathrm{}`$ splits into two sech pulses. As time moves on, the front $`u_3`$ sech pulse breaks into $`u_1`$ and $`u_2`$ sech pulses. Then this $`u_2`$ sech pulse and the back $`u_3`$ sech pulse interact. The final outcome is two $`u_1`$ sech pulses moving in the positive $`x`$ direction, and a higher-order $`u_2`$ wave moving in the negative $`x`$ direction. Thirdly, we consider the non-generic case where $`\theta _1^{(1)}=0`$, while $`\theta _2^{(1)}`$ and $`\theta _3^{(1)}`$ still do not change. This soliton solution is illustrated in Fig. 3. We see that as $`t\mathrm{}`$, both the $`u_3`$ and $`u_2`$ waves are non-zero and sech-shaped. After their interaction, the pumping $`u_3`$ wave is depleted, and a new $`u_1`$ sech wave and a higher-order $`u_2`$ wave are created. We note that this $`u_2+u_3u_1+u_2`$ process is novel, and it has not been carefully investigated before. Fourthly, we consider the non-generic case where $`\theta _1^{(1)}=\theta _2^{(1)}=0`$, and $`\theta _3^{(1)}`$ is still $`1`$. This solution is illustrated in Fig. 4. We see that it is the same as a fundamental soliton solution, and it describes the process of a pumping $`u_3`$ sech wave breaking into two elementary $`u_1`$ and $`u_2`$ sech waves. Thus, our higher-order soliton solution reduces to a fundamental soliton solution as a special case. Lastly, we consider the non-generic case where $`\theta _2^{(1)}=0`$ while $`\theta _1^{(1)}=1`$ and $`\theta _3^{(1)}=1`$ as in Fig. 1. This solution is shown in Fig. 5. As we can see, as $`t\mathrm{}`$, the only non-zero wave is the pumping wave $`u_3`$, which is a higher-order soliton. As $`t\mathrm{}`$, this pumping wave breaks up into a sech waves in each component. Thus, this is the new $`u_3u_1+u_2+u_3`$ process which we have presented in the text above.

We conclude this section with some comments on the soliton solutions to the three-wave model corresponding to the higher-order zeros of order $`n2`$. If the higher-order zero is elementary, i.e. when the sequence of ranks in formula (48) is rank$`P_j=1`$, $`j=1,\mathrm{},n`$, the corresponding soliton solutions can be derived using the soliton matrix (77). However, there are two other possible sequences of ranks in formula (48), namely

$`(a)\text{rank}P_j=2,j=1,\mathrm{},r,n=2r;`$

$`(b)\text{rank}P_j=2,j=1,\mathrm{},r;\text{rank}P_j=1,j=r+1,\mathrm{},r+s,n=2r+s.`$

We note that the soliton matrix for the sequence of ranks (a), which corresponds to the higher-order zero of order $`2r`$, has an equivalent soliton matrix corresponding to the elementary higher-order zero of order $`r`$. Indeed, let us consider the soliton matrix in the representation (47), where each $`\chi _j(k)`$ is defined similar to formula (45) \[$`P_j`$ substituted for $`P_1`$\] with rank$`P_j=2`$. Consider the following procedure. First, define new projectors $`Q_j=IP_j`$. Evidently rank$`Q_j=1`$. Second, multiply the soliton matrix $`\mathrm{\Gamma }(k)`$ (47) by a scalar quotient,

$$\stackrel{~}{\mathrm{\Gamma }}(k)=\left(\frac{k\overline{k}_1}{kk_1}\right)^r\mathrm{\Gamma }(k),$$
(150)
such that each $`\chi _j(k)`$, $`j=1,\mathrm{},r`$, gets a multiplier $`(k\overline{k}_1)/(kk_1)`$. We have

$`\stackrel{~}{\chi }_j(k){\displaystyle \frac{k\overline{k}_1}{kk_1}}\chi _j(k)={\displaystyle \frac{k\overline{k}_1}{kk_1}}\left({\displaystyle \frac{kk_1}{k\overline{k}_1}}I+{\displaystyle \frac{k_1\overline{k}_1}{k\overline{k}_1}}Q_j\right)=I+{\displaystyle \frac{k_1\overline{k}_1}{kk_1}}Q_j,`$

and

$$\stackrel{~}{\mathrm{\Gamma }}(k)=\stackrel{~}{\chi }_r(k)\mathrm{}\stackrel{~}{\chi }_2(k)\stackrel{~}{\chi }_1(k).$$
(151)
Evidently, the new matrix $`\stackrel{~}{\mathrm{\Gamma }}(k)`$ in (151) satisfies the linear system of equations (87)-(88) for the original matrix $`\mathrm{\Gamma }(k)`$. Furthermore, it corresponds to an elementary higher-order zero of order $`r`$, though now in the complementary half plane: $`k=\overline{k}_1`$. It is noted that in this section we considered a soliton matrix $`\mathrm{\Gamma }(k)`$ corresponding to a zero $`k_1=\xi +i\eta `$ lying in the upper half plane, i.e. with $`\eta >0`$. However, the case of $`\eta <0`$ is admissible as well. The only significant change would be in the asymptotic formulae, and the effect of this change is similar to reversing the time variable: $`tt`$. Thus, the sequence of ranks in case (a) brings no new higher-order soliton solutions as compared to the simple sequence of ranks, but the solution process is reversed. For the fundamental soliton solutions, a similar fact has been noted in Ref. , where it is mentioned that the fundamental soliton corresponding to the projector of rank 2 describes the three-wave interaction process which is reverse to that of the soliton solution corresponding to the projector of rank 1.

There is no transformation similar to (150) for case (b) (similar multiplication will produce a rational matrix function having poles in both half planes, thus such a matrix does not belong to the class of soliton matrices). Higher-order soliton solutions in this case require construction of the soliton matrices for non-elementary higher-order zeros and will be addressed in a forthcoming paper.

## V Conclusion

We have proposed a unified and systematic approach to study the higher-order soliton solutions of nonlinear PDEs integrable by the Riemann-Hilbert problem of arbitrary matrix dimension. We have derived the soliton dressing matrix for the elementary higher-order zeros in the $`N\times N`$-dimensional spectral problem, i.e., zeros having the geometric multiplicity 1. The associated higher-order solitons in the $`N`$-wave system have also been obtained. We have also clarified that the soliton dressing ansatz proposed in  is the general soliton matrix for the nonlinear Schrödinger equation (where $`N=2`$), thus the soliton solutions obtained in  are the most general higher-order solitons in the the nonlinear Schrödinger equation. For $`N\times N`$-dimensional spectral problems the soliton dressing ansatz of  corresponds to the elementary higher-order solitons.

We have applied our theory to the three-wave interaction model, and the simplest higher-order soliton solution has been obtained. The generic case of this solution describes the process $`u_3u_1+u_2`$, similar to fundamental solitons. But each wave involved here is higher-order. The non-generic case of this solution could describe three new processes. The first two are similar to each other:

$`u_1+u_3u_1+u_2,`$
$`u_2+u_3u_1+u_2.`$
Here the waves on the left are all sech waves; the waves on the right are a sech wave and a higher-order wave. The third process reads

$`u_3u_1+u_2+u_3,`$

where the pumping wave on the left is a higher-order wave, and the waves on the right all have sech shape. The non-generic solutions could also reduce to fundamental solitons or trivial solutions as special cases.

We anticipate that the higher-order soliton solutions will have wide applications. First of all, the new processes they describe may find physical applications where three-wave interaction takes place. Second, as it has been mentioned in Ref. , the higher-order soliton solution describes a weak bound state of solitons, thus it may appear in the study of the train propagation of solitons with nearly equal amplitudes and velocities in nonlinear integrable PDEs. The usual approach in the analytical study of the soliton trains is reduction of the governing equations for the soliton parameters to the complex Toda chain (consult, for instance, Refs. ). The higher-order soliton approach may provide an alternative to this study. Thirdly, multi-hump solitary waves in the non-integrable nonlinear PDEs can be another field of application of the higher-order solitons. For instance, the so-called multisoliton complexes, or more precisely, oscillatory and stationary solitons observed in an oscillating water trough  and subsequently reproduced in numerical simulations  of the governing parametrically driven, damped NLS equation may have the same relation to the higher-order solitons as the usual solitary-wave solutions of the non-integrable PDEs to the fundamental solitons. Analytical study of the soliton complexes needs the perturbation theory for the higher-order solitons, just as the study of usual solitary-wave solutions needs the perturbation theory for the fundamental solitons. The perturbation theory for the higher-order solitons can be developed in a similar way as it is done for the fundamental solitons (see for instance Ref. ). Such a theory is left for future studies.

Lastly, we point out that the soliton matrices for the elementary zeros serve as the building blocks for the general case of zeros with arbitrary geometric multiplicity. This work is in progress and will be reported in a forthcoming paper. There the most general higher-order soliton solution for the $`N`$-wave system will be given.

###### Acknowledgements.

Helpful discussions with Professor Mark Ablowitz are gratefully acknowledged. The work by V.S. was sponsored by the University of Vermont. The warm hospitality of the Department of Mathematics and Statistics is also acknowledged with gratitude. The work by J.Y. was supported in part by AFOSR under contract USAF F49620-99-1-0174 and by NSF under grant DMS-9971712.

## A Exponent of the Toeplitz matrices

Here we show that for a diagonal matrix $`M(k)`$ the exponent of the block Toeplitz matrix, defined as in formula (92), and the block Toeplitz matrix of the exponent of $`M(k)`$ coincide. As the derivatives $`\frac{\text{d}^m}{\text{d}k^m}M(k)`$ commute with each other, it is enough to prove this statement for a scalar function. Consider, for example, the lower-triangular Toeplitz matrix of a scalar function $`f(k)`$:

$$𝐅=\left(\begin{array}{cccc}f& 0& \mathrm{}& 0\\ \frac{\text{d}}{\text{d}k}f& f& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{\text{d}^{n1}}{\text{d}k^{n1}}f& \frac{\text{d}^{n2}}{\text{d}k^{n2}}f& \mathrm{}& f\end{array}\right).$$
(A1)
It can be rewritten in the following form

$$𝐅=H_0f+H_1\frac{\text{d}f}{\text{d}k}+\mathrm{}+H_{n1}\frac{\text{d}^{n1}f}{\text{d}k^{n1}},\left(H_j\right)_{l,m}\delta _{l+j,m}.$$
(A2)
Note the product rule for the “diagonals”: $`H_jH_i=H_{j+i}`$ and that for $`j+i>n1`$ the product is zero. Therefore, the exponent of $`𝐅`$ is a finite sum of the diagonals $`H_j`$:

$$\mathrm{exp}(𝐅)=c_0H_0+c_1H_1+\mathrm{}+c_{n1}H_{n1},$$
(A3)
where $`c_0,\mathrm{},c_{n1}`$ are constants. Due to the formula $`H_j=H_1^j`$, $`j=0,\mathrm{},n1`$, computing the coefficients $`c_j`$ is equivalent to taking the finite sum of the first $`n`$ terms of the Taylor expansion of an equivalent scalar function:

$$\mathrm{exp}\left\{\underset{j=0}{\overset{n1}{}}\frac{\text{d}^jf(k)}{\text{d}k^j}ϵ^j\right\}=c_0+c_1ϵ+\mathrm{}+c_{n1}ϵ^{n1}+𝒪\{ϵ^n\},$$
(A4)
where $`ϵ`$ is the parameter of the Taylor expansion which represents $`H_1`$. On the other hand, computing the Taylor expansion reduces to taking derivatives with respect to $`k`$ of $`\mathrm{exp}\{f(k+ϵ)\}`$ at $`ϵ=0`$:

$`\mathrm{exp}\left\{{\displaystyle \underset{j=0}{\overset{n1}{}}}{\displaystyle \frac{\text{d}^jf(k)}{\text{d}k^j}}ϵ^j\right\}=\mathrm{exp}\left\{f(k+ϵ)\right\}+𝒪\{ϵ^n\}`$ (A5)
$`=\mathrm{exp}\{f(k)\}+{\displaystyle \frac{1}{1!}}{\displaystyle \frac{\text{d}}{\text{d}k}}\mathrm{exp}\{f(k)\}ϵ+\mathrm{}+{\displaystyle \frac{1}{(n1)!}}{\displaystyle \frac{\text{d}^{n1}}{\text{d}k^{n1}}}\mathrm{exp}\{f(k)\}ϵ^{n1}+𝒪\{ϵ^n\}.`$ (A6)
Therefore

$$c_j=\frac{1}{j!}\frac{\text{d}^j}{\text{d}k^j}\mathrm{exp}\{f(k)\},j=1,\mathrm{},n1,$$
(A7)
thus,

$$\mathrm{exp}\{𝐅\}=H_0\mathrm{exp}\{f(k)\}+H_1\frac{1}{1!}\frac{\text{d}}{\text{d}k}\mathrm{exp}\{f(k)\}+\mathrm{}+H_{n1}\frac{1}{(n1)!}\frac{\text{d}^{n1}}{\text{d}k^{n1}}\mathrm{exp}\{f(k)\}.$$
(A8)
Q.E.D.