In 1968 one of the great discoveries in mathematical physics took place. Its authors: P. Lax, C. S. Gardner, J. M. Greene, M. D. Kruskal, R. M. Miura and N. J. Zabusky after several years of analysis proved that the KdV equation can be exactly solved by the inverse scattering method (ISM). This was the first and for some time, the only NLEE that could be solved exactly. Soon after that it was demonstrated that the KdV is completely integrable infinite-dimensional Hamiltonian system; its action-angle variables were found by Zakharov and Faddeev [137]. The whole story is well described by N. J. Zabusky in his review paper [128].

The second big step in this direction followed in 1971 by the seminal paper of Zakharov and Shabat who discovered the second equation integrable by the ISM: the nonlinear Schrödinger (NLS) equation [132]; in 1973 the same authors demonstrated that the NLS equation is integrable also under nonvanishing boundary conditions [133]. Both versions of NLS equations described interesting and important physical applications in nonlinear optics, plasma physics, hydrodynamics and others. This inspired many scientists, mathematicians and physicists alike to join the scientific community interested in the study of soliton equations. As a result new soliton equations started to appear one after another. Here we only mention the modified KdV (mKdV) equation [125], the N-wave equations [131], the Manakov system known also as the vector NLS equation [92] and many others. Many of them have already been included in monographs: see e.g. [1,10,18,64,104] and the numerous references therein.

The first few NLEE were related to the algebra $sl\left(2\right)$, so the corresponding inverse scattering problem could be solved using the famous Gelfand-Levitan-Marchenko (GLM) equation. For the Manakov system it was necessary to use $2\times 2$ block matrix Lax operator, so the GLM equation was naturally generalized also for that case. However, the ISM for the N-wave system came up to be substantially more difficult. Indeed, for the $2\times 2$ (block) matrix Lax operator the Jost solutions possess analyticity properties which are basic for the GLM eq. However the Lax pair for the N-wave system is the generalized $sl\left(n\right)$ Zakharov-Shabat system: where $a_k$ and $b_k$ are real constants such that $\mathrm{tr}J=0$, $\mathrm{tr}K=0$. Without restrictions we can assume that $a_1>a_2>\cdots >a_n$. In this case only the first and the last column of the corresponding Jost solutions allow analytic extension in the spectral parameter $\lambda$. This, however, was not enough to derive GLM equation. It was Shabat who discovered the way out of this difficulty [106,107]. He was able to modify the integral equations for the Jost solutions into integral equations that provide the fundamental analytic solutions (FAS) ${\chi }^{+}(x,t,\lambda )$ and ${\chi }^{+}(x,t,\lambda )$ of L which allowed analytic extensions for $\mathrm{Im}\lambda >0$ and $\mathrm{Im}\lambda <0$ respectively. As a result the interrelation between the FAS and the sewing function $G_0(x,t,\lambda )$: can be reformulated as Riemann-Hilbert problem (RHP). Now we can solve the ISP for L by using the RHP with canonical normalization: The normalization condition (1.4) ensures that the RHP has unique regular solution [104]. It also allowed Zakharov and Shabat to develop their dressing method, which enables one to calculate the N-soliton solutions not only for the N-wave system, but also for the whole hierarchy of NLEE related to L (1.1); see [134,135] and also [32,104]. In short, the dressing Zakharov-Shabat factor came up as one of the most effective methods for a) constructing soliton solutions, and b) understand that the dressed Lax operator has some additional discrete eigenvalues, as well as the explicit form of the dressed FAS.

The GZS system has natural reductions after which the potential $Q(x,t)$ and J belong either to $so\left(n\right)$ or to $sp\left(2r\right)$ algebra. Other reductions were proposed by Mikhailov [96] which substantially enlarged the classes of integrable NLEE. Some of these reductions require that J has complex-valued eigenvalues. Constructing FAS for such systems poses additional difficulties, which were overcome by Beals and Coifman [7] for the GZS systems related to $sl\left(n\right)$ algebra. Later the results of [7] were extended first for the systems related to $so\left(n\right)$ or to $sp\left(2r\right)$ algebra, (see [64]) as well as to Mikhailov’s reductions [57,64].

The FAS play an important role in soliton theory. Indeed, they can be used to introduce

The above arguments lead us to the hypothesis that we could use more effective approach to the integrable NLEE which starts from the RHP rather than from a specific Lax operator. In the first part of this paper we will demonstrate that FAS could be constructed and used also for quadratic pencils. We also formulate explicitly the corresponding RHP. For quadratic pencils we have additional natural symmetry which maps $\lambda \to -\lambda$. This symmetry is also inherent in the contour in the complex $\lambda$-plane, on which the RHP is defined.

In Section 2 we first demonstrate that the well known methods for analysis of NLEE can be generalized also for Lax operators quadratic in $\lambda$, see eq. (2.28) below. On this level we for the first time meet with purely algebraic problem for constructing two commuting quadratic pencils. For polynomial pencils of higher orders those problems will be more and more difficult to solve. In particular we outline the construction of the FAS ${\chi }^{+}(x,t,\lambda )$ and ${\chi }^{-}(x,t,\lambda )$ which are analytic in ${\mathsf{\Omega }}_0\cup {\mathsf{\Omega }}_3$ and ${\mathsf{\Omega }}_2\cup {\mathsf{\Omega }}_4$ respectively, see Figure 4 below. As a result, the continuous spectrum of these Lax operators fill up the union of the real and imaginary axis of the complex $\lambda$-plane. As a consequence contours of the corresponding RHP will be $\mathbb{R}\cup i\mathbb{R}$ unless an additional factor complicates the picture. Thus we see that the symmetries of the NLEE or of its Lax pair determine the contour of the RHP.

In Section 3 we remind the notion of Mikhailov’s reduction group and briefly outline characteristic contours of the relevant RHP. We also demonstrate that Zakharov–Shabat theorem is valid for a larger class of Lax operators than it was proved before. It is important to request that the RHP is canonically normalized. This ensures that the RHP has unique regular solution, which is important for the application of the Zakharov-Shabat dressing method.

Another important factor in formulating the RHP is Mikhailov’s reduction group. In Section 4 we outline some of the obvious effects which the reduction group may have on the contour of the RHP. Therefore it is not only the order of the polynomial in $\lambda$, but also the symmetry (the reduction group) which determine the contour of the RHP. For example, if we add an additional Mikhailov’s symmetry that maps $\lambda \to 1/{\lambda }^{*}$ then the corresponding Lax operator will be polynomial in $\lambda$ and $1/\lambda$ which in turn will require adjustments in the techniques for deriving the dressing factors and soliton solutions.

In Section 5 we propose a parametrization of the solutions $\xi (x,t,\lambda )$ of RHP for the class of RHP related to homogeneous spaces, see eq. (5.1). Here we require that the coefficients $Q_s(x,t)$ provide local coordinates of the corresponding homogeneous space. Thus we are able to derive a new systems of N-wave equations, see also [36,52,53,56,74]. We also demonstrate that the dressing Zakharov-Shabat method [100,129,134,135,136] can be naturally extended to derive the soliton solutions of these new N-wave equations. At the same time the structure of the dressing factors depends substantially on the symmetries of the Lax operators. Thus even for the one-soliton solutions we need to solve linear block-matrix equations. The situation when we have two involutions: the Hermitian one ${\left({\chi }^{+}(x,t,\lambda )\right)}^{†}={\left({\chi }^{-}(x,t,{\lambda }^{*})\right)}^{-1}$ and the $\lambda \to -\lambda$ symmetry $C_0{\chi }^{\pm }(x,t,\lambda )C_0^{-1}={\chi }^{\pm }(x,t,-\lambda )$ are typical for Lax operators L related to the algebras $sl(n,\mathbb{C})$. But if we request in addition that L is related to symplectic or orthogonal Lie algebra then we have to deal with three involutions, and the corresponding linear equations get more involved. That is why we focus first on the one-soliton solutions. The derivation of the N soliton solutions is discussed later.

Section 6 is devoted to the MNLS equations which require the use of symmetric spaces, see Refs. [21,40,67,92,94,118,119,122,126]. We start again with the parametrization of the RHP which now must be compatible with the structure of the symmetric spaces. To us it was natural to limit ourselves to the four classes Hermitian symmetric spaces related to the non-exceptional Lie algebras, see [67]. Again we parametrize the coefficients $Q_s(x,t)$ as local coordinates of the corresponding symmetric spaces. In fact $\mathcal{Q}(x,t,\lambda )$ must have the same grading as the symmetric space, but we were able to apply additional reductional requesting $Q_{2s}=0$, see eq. (6.11) below. Thus we formulate the typical MNLS equations related to the four classes of symmetric spaces.

In Section 7 we derive the one soliton solutions of MNLS. Again, like in Section 5, we treat separately the MNLS related to A.III type symmetric spaces, because the corresponding FAS have only two involutions. The MNLS related to C.I and D.III symmetric spaces possess three involutions; the corresponding linear equations are similar to the ones for the class of N-wave equations, but the solutions are different. The symmetric spaces of BD.I class are treated separately, because their typical representation is provided by $3\times 3$ block matrices, so many of the calculations are indeed different. At the end of this Section we derive the soliton interactions for the BD.I class of MNLS [41]. More precisely we use the asymptotic of the dressing factor for $x\to \pm \infty$ applying it to the two-soliton solution in order to calculate the center-of-mass and the phase shifts of the solitons.

In Section 8 we introduce the resolvent of the Lax operators in terms of the FAS. The diagonal of the resolvent after a regularization can be expressed in terms of the solution of the RHP by ${\xi }^{\pm }(x,t,\lambda )J{\widehat{\xi }}^{\pm }(x,t,\lambda )$; here by 'hat' we denote the inverse matrix. It can be viewed as generating functional of the integrals of motion.

We end the paper by discussion and conclusions. Some technical aspects in the calculations such as the structure of the symmetric spaces, and the root systems of the simple Lie algebras as well as the Gauss decompositions of the elements of the simple Lie groups are given in the appendices.