Integrability and Linear Stability of Nonlinear Waves

It is well known that the linear stability of solutions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1+1$$\end{document}1+1 partial differential equations which are integrable can be very efficiently investigated by means of spectral methods. We present here a direct construction of the eigenmodes of the linearized equation which makes use only of the associated Lax pair with no reference to spectral data and boundary conditions. This local construction is given in the general \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times N$$\end{document}N×N matrix scheme so as to be applicable to a large class of integrable equations, including the multicomponent nonlinear Schrödinger system and the multiwave resonant interaction system. The analytical and numerical computations involved in this general approach are detailed as an example for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=3$$\end{document}N=3 for the particular system of two coupled nonlinear Schrödinger equations in the defocusing, focusing and mixed regimes. The instabilities of the continuous wave solutions are fully discussed in the entire parameter space of their amplitudes and wave numbers. By defining and computing the spectrum in the complex plane of the spectral variable, the eigenfrequencies are explicitly expressed. According to their topological properties, the complete classification of these spectra in the parameter space is presented and graphically displayed. The continuous wave solutions are linearly unstable for a generic choice of the coupling constants.

ties, the complete classification of these spectra in the parameter space is presented and graphically displayed. The continuous wave solutions are linearly unstable for a generic choice of the coupling constants.

Introduction
The problem of stability is central to the entire field of nonlinear wave propagation and is a fairly broad subject. Here, we are specifically concerned with the early stage of amplitude modulation instabilities due to quadratic and cubic nonlinearities, and we consider in particular dispersive propagation in a one-dimensional space, or diffraction in a two-dimensional space.
After the first observations of wave instability (Benjamin and Feir 1967;Rothenberg 1990Rothenberg , 1991; see also e.g. Zakharov and Ostrovsky 2009), the research on this subject has grown very rapidly because similar phenomena appear in various contexts such as water waves (Yuen and Lake 1980), optics (Agrawal 1995), Bose-Einstein condensation (Kevrekidis et al. 2007) and plasma physics (Kuznetsov 1977). Experimental findings were soon followed by theoretical and computational works. Predictions regarding short time evolution of small perturbations of the initial profile can be obtained by standard linear stability methods (see e.g. Skryabin and Firth 1999 and references therein). Very schematically, if u(x, t) is a particular solution of the wave equation, evolving in time t, and if u+δu is the perturbed solution of the same equation, then, at the first order of approximation, δu satisfies a linear equation whose coefficients depend on the solution u(x, t) itself, and therefore, they are generally non-constant. Consequently, solving the initial value problem δu(x, 0) → δu (x, t) in general is not tractable by analytical methods. It is only for special solutions u(x, t), such as nonlinear plane waves or solitary localized waves, see e.g. Skryabin and Firth (1999), that this initial value problem can be approached by solving an eigenvalue problem for an ordinary differential operator in the space variable x. In this way, the computational task reduces to constructing the eigenmodes, i.e. the eigenfunctions of an ordinary differential operator, while the corresponding eigenvalues are simply related to the proper frequencies. For very special solutions u(x, t), this procedure exceptionally leads to a linearized equation with constant coefficients which can be solved therefore via Fourier analysis. A simple and well-known example of this case is the linearization of the focusing nonlinear Schrödinger (NLS) equation iu t + u x x + 2|u| 2 u = 0 around its continuous wave (CW) solution u(x, t) = e 2it . Here and thereafter, subscripts x and t denote partial differentiation, unless differently specified. The computation of all the complex eigenfrequencies, in particular of their imaginary parts, yields the relevant information about the stability of u(x, t), provided the set of eigenmodes be complete in the functional space characterized by the boundary conditions satisfied by the initial perturbation δu(x, 0). Since the main step of this method is that of finding the spectrum of a differential operator, the stability property of the solution u(x, t) is also referred to as spectral stability. It is clear that this method applies to a limited class of solutions of the wave equation. Although herein we are concerned with linear stability only, quite a number of studies on other forms of nonlinear waves stability have been produced by using different mathematical techniques and aimed at various physical applications. For instance, variational methods to assess orbital stability have been applied to solitary waves and standing waves (e.g., see Maddocks and Sachs 1993;Georgiev and Ohta 2012).
An alternative and more powerful approach to stability originated from Ablowitz et al. (1974) shortly after the discovery of the complete integrability and of the spectral method to solve the Korteweg-de Vries (KdV) and NLS equations (e.g. see the textbooks Ablowitz and Segur 1981;Calogero and Degasperis 1982;Novikov et al. 1984). This method stems from the peculiar fact that the so-called squared eigenfunctions (see next section for their definition) are solutions of the linearized equation solved by the perturbation δu (x, t). Indeed, depending on boundary conditions, this technique yields a representation of the perturbation δu (x, t) in terms of such squared eigenfunctions. With respect to the spectral methods in use for non-integrable wave equations, the squared eigenfunctions approach to stability shows its power by formally applying to almost any solution u (x, t), namely also to cases where standard methods fail. Moreover, this method, with appropriate algebraic conditions, proves to be applicable (see Sect. 2) to a very large class of matrix Lax pairs and, therefore, to quite a number of integrable systems other than KdV and NLS equations (e.g. sine-Gordon, mKdV, derivative NLS, coupled NLS, three-wave resonant interaction, massive Thirring model and other equations of interest in applications). Its evident drawback is that its applicability is limited to the very special class of integrable wave equations. Notwithstanding this condition, it remains of important practical interest because several integrable partial differential equations have been derived in various physical contexts as reliable, though approximate, models (Dodd et al. 1982;Ablowitz and Segur 1981;Dauxois and Peyrard 2006). Moreover, the stability properties of particular solutions of an integrable wave equation provide a strong insight about similar solutions of a different non-integrable, but close enough, equation. Among the many properties defining the concept of integrability the one that we consider here is the existence of a Lax pair of two linear ordinary differential equations for the same unknown, one in the space variable x and the other in the time variable t (see next Section), and whose compatibility condition is just the wave equation. Thus, a spectral problem with respect to the variable x already appears at the very beginning of the integrability scheme. With appropriate specifications, as stated below, this observation leads to the construction of the eigenmodes of the linearized equation, in terms of the solutions of the Lax pair. Moreover, via the construction of the squared eigenfunctions one is able to compute the corresponding eigenfrequency ω, which gives the (necessary and sufficient) information to assess linear stability by the condition Im(ω) > 0. Explicit expressions of such eigenmodes have been obtained if the unperturbed wave amplitude u is a cnoidal wave (e.g. see Kuznetsov and Mikhailov 1974;Sachs 1983;Kuznetsov et al. 1984 for the KdV equation and Kuznetsov and Spector 1999 for the NLS equation), or if it is a soliton solution (Yang 2000) or, although only formally, an arbitrary solution (Yang 2002). Therefore, the computational strategy amounts to constructing the set of eigenmodes and eigenfrequencies. It should be pointed out that the integrability methods, in an appropriate functional space of the wave fields u (x, t), provide also the way of deriving the closure and completeness relations of the eigenmodes, see e.g. Kaup (1976a), Yang (2000Yang ( , 2002 for solutions which vanish sufficiently fast as |x| → 0. In this respect, a word of warning is appropriate. The boundary conditions imposed on the solutions u(x, t) play a crucial role in proving that the wave evolution be indeed integrable. Thus, in particular for the NLS equation, integrability methods have been applied so far to linear stability of wave solutions which, as |x| → ∞, either vanish as a localized soliton (Yang 2000), or go to a CW solution (see the lecture notes Degasperis and Lombardo 2016), or else are periodic, u(x, t) = u(x + L , t) (Bottman et al. 2009). In these cases, by solving the so-called direct spectral problem, to any solution u(x, t) one can associate a set of spectral data, the spectral transform, say the analogue of the Fourier transform in a nonlinear context. This correspondence allows to formally solve the initial value problem of the wave equation. As a by-product, this formalism yields also a spectral representation of the small perturbations δu (x, t) in terms of the corresponding small change of the spectral data. This connection is given by the squared eigenfunctions (see Ablowitz et al. 1974;Kaup 1976a;Yang and Kaup 2009 for the NLS equation and Calogero and Degasperis 1982 for the KdV equation) which play the role which the Fourier exponentials have in the linear context. Indeed, the squared eigenfunctions, which are computed by solving the Lax pair, are the eigenmodes of the linearized equation for δu (x, t). This result follows from the inverse spectral transform machinery. However, as we show below, the squared eigenfunctions' property of being solutions of the linearized equation is a local one, as it follows directly from the Lax pair without any need of the spectral transform. More than this, integrability allows to go beyond the linear stage of the evolution of small perturbations. This is possible by the spectral method of solving the initial value problem for the perturbed solution u + δu which therefore yields the long time evolution of δu beyond the linear approximation (see, for instance, Zakharov and Gelash 2013;Biondini and Mantzavinos 2016). However, this important problem falls outside the scope of the present work and it will not be considered here (for the initial value problem and unstable solutions of the NLS equation, see Santini 2017, 2018a, b) The stability properties of a given solution u(x, t) may depend on parameters. These parameters come from the coefficients of the wave equation and from the parameters (if any) which characterize the solution u(x, t) itself. This obvious observation implies that one may expect the parameter space to be divided into regions where the solution u(x, t) features different behaviours in terms of linear stability. Indeed, this is the case, and crossing the border of one of these regions by varying the parameters, for instance a wave amplitude, may correspond to the opening of a gap in the instability frequency band, so that a threshold occurs at that amplitude value which corresponds to crossing. The investigation of such thresholds is rather simple when dealing with scalar (one-component) waves. For instance, the KdV equation has no frozen coefficient, for a simple rescaling can set them equal to any real number, so that it reads u t + u x x x + uu x = 0. On the other hand, after rescaling, the NLS equation comes with a sign in front of the cubic term, distinguishing between defocusing and focusing self-interaction. These two different versions of the NLS equation lead to different phenomena such as modulation stability and instability of the continuous wave solution (for an introductory review, see Degasperis and Lombardo 2016). Wave propagation equations which model different physical systems may have more structural coefficients whose values cannot be simultaneously, and independently, rescaled. This is the case when two or more waves resonate and propagate while coupling to each other. In this case, the wave equations do not happen to be integrable for all choices of the coefficients. A well-known example, which is the focus of Sect. 3, is that of two interacting fields, u j , j = 1, 2, which evolve according to the coupled system of NLS equations where (s 1 |u 1 | 2 +s 2 |u 2 | 2 ) is the self-and cross-interaction term. This is integrable only in three cases (Zakharov and Shulmann 1982), namely (after appropriate rescaling): s 1 = s 2 = 1 (defocusing Manakov model), s 1 = s 2 = −1 (focusing Manakov model) (Manakov 1973) and the mixed case s 1 = −s 2 = 1. These three integrable systems of two coupled NLS (CNLS) equations are of interest in few special applications in optics (Menyuk 1987;Evangelides et al. 1992;Wang and Menyuk 1999) and in fluid dynamics (Onorato et al. 2010), while, in various contexts (e.g. in optics (Agrawal 1995) and in fluid dynamics (Yuen and Lake 1980;Ablowitz and Horikis 2015)), the coupling constants s 1 , s 2 take different values and the CNLS system happens to be non-integrable. Yet the analysis of the three integrable cases is still relevant in the study of the (sufficiently close) non-integrable ones (Yang and Benney 1996). The linear stability of CW solutions, |u j (x, t)| = constant, of integrable CNLS systems is of special interest not only because of its experimental observability, but also because it can be analysed via both standard methods and the squared eigenfunctions approach. As far as the standard methods are concerned, the linear stability of CW solutions has been investigated only for the focusing and defocusing regimes, but not for the mixed one (s 1 = −s 2 ), and only in the integrable cases, by means of the Fourier transform (Forest et al. 2000). Conversely, as far as the integrability methods are concerned, it has been partially discussed in (Ling and Zhao 2017) to mainly show that instability may occur also in defocusing media, in contrast to scalar waves which are modulationally unstable only in the focusing case. In the following we approach the linear stability problem of the CW solutions of (1) within the integrability framework to prove that the main object to be computed is a spectrum (to be defined below) as a curve in the complex plane of the spectral variable, together with the eigenmodes wave numbers and frequencies defined on it. In particular, we show that the spectrum which is relevant to our analysis is related to, but not coincident with, the spectrum of the Lax equation for . In addition, if λ is the spectral variable, the computational outcome is the wave number k(λ) and frequency ω(λ), so that the dispersion relation and also the instability band are implicitly defined over the spectrum through their dependence on λ. Since spectrum and eigenmodes depend on parameters, we explore the entire parameter space of the two amplitudes and coupling constants to arrive at a complete classification of spectra by means of numerically assisted, algebraic-geometric techniques. Our investigation in Sect. 3 illustrates how the linear stability analysis works within the theory of integrability. Our focus is on x-and t-dependent CWs, a case which is both of relevance to physics and is computationally approachable. This case is intended to be an example of the general method developed in Sect. 2. It is worth observing again that the linear stability of the CW solutions can indeed be discussed also by standard Fourier analysis, e.g., see Forest et al. (2000) for the CNLS systems in the focusing and defocusing regimes. However, such analysis is of no help to investigate the stability of other solutions. On the contrary, at least for the integrable CNLS system (1), our method relies only on the existence of a Lax pair, and as such, it has the advantage of being applicable also to other solutions as well. In particular, it can be applied to the CW solutions in all regimes (as we do it here), as well as to solutions such as dark-dark, bright-dark and higher-order solitons travelling on a CW background, to which the standard methods are not applicable. This article is organized as follows. In the next section (Sect. 2), we give the general (squared eigenfunctions) approach together with the expression of the eigenmodes of the linearized equation for the N × N matrix scheme, so as to capture a large class of integrable systems. There we define the x-spectrum in the complex plane of the spectral variable. In Sect. 3, we provide an example of application of the theory by specializing the formalism introduced in Sect. 2 to deal with the CNLS equations. We characterize the x-spectrum in the complex plane of the spectral variable according to their topological features, and we cover the entire parameter space according to five distinct classes of spectra. This characterization of the spectrum holds under the assumption that the small perturbation of the background CW solution is localized. Sect. 3.3 is devoted to discussing the classification of spectra and the corresponding stability features in the focusing, defocusing and mixed coupling regimes, in terms of the physical parameters, while a conclusion with open problems and perspectives of future work is the content of Sect. 4. Details regarding computational and numerical aspects of the problem are confined in Appendices.

Integrable Wave Equations and Small Perturbations
The integrable partial differential equations (PDEs) which are considered here are associated with the following pair of matrix ordinary differential equations (ODEs), also known as Lax pair (e.g. see Calogero and Degasperis 1982;Novikov et al. 1984;Ablowitz and Clarkson 1991), where , X and T are N × N matrix-valued complex functions. The existence of a fundamental (i.e. non singular) matrix solution = (x, t) of this overdetermined system is guaranteed by the condition that the two matrices X and T satisfy the differential equation We recall here that, unless differently specified, a subscripted variable means partial differentiation with respect to that variable, and [A , B] stands for the commutator AB − B A. In order to identify this condition (3) as an integrable partial differential equation for some of the entries of the matrix X , it is essential that both matrices X and T parametrically depend on an additional complex variable λ, known as the spectral parameter. In order to make this introductory presentation as simple as possible, we assume that X (λ) and T (λ) be polynomial in λ with degrees n and m, respectively. As a consequence, the matrix X t − T x + [X , T ] is as well polynomial in λ with degree n + m, and therefore, compatibility Eq. (3) yields n + m + 1 equations for the matrix coefficients of the polynomials X and T . If the pair X and T is a given solution of (3), we consider a new solution X → X + δ X , T → T + δT , which differs by a small change of the matrices X and T , with the implication that the pair of matrices δ X and δT , at the first order in this small change, satisfies the linearized equation Again, the left-hand side of this linearized equation has a polynomial dependence on λ and the vanishing of all its coefficients results in a number of algebraic or differential equations. These obvious observations lead us to focus on matrix linearized Eq. (4) itself, which reads, by setting A = δ X , B = δT , and to search for its solutions A(x, t, λ) and B(x, t, λ). Our main target is to find those solutions which are related to the fundamental matrix solution (x, t, λ) of Lax pair (2). To this purpose, we first note that the similarity transformation of a constant (i.e. x, t-independent) matrix M(λ) yields the transformed matrix which, for any given arbitrary matrix M, satisfies the pair of linear ODEs Equation (7) are compatible with each other because of (3). Then, for future reference, we point out the following observations.
Proposition 1 If the pair A , B solves linearized Eq. (5), then also the pair is a solution of the same linearized Eq. (5), namely This is a straight consequence of the Jacobi identity and of the assumption that the matrix be a solution of (7).

Proposition 2
The following expressions provide a solution of linearized Eq. (9).
The validity of this statement follows from the fact that the matrices obviously solve linearized Eq. (5) and from Proposition 1. At this point, we go back to the nonlinear matrix PDE which follows from condition (3). In view of the applications within the theory of nonlinear resonant phenomena that we have in mind (which include the multicomponent nonlinear Schrödinger system and the multiwave resonant interaction system), we assume that the polynomials X (λ) and T (λ), see (2), have degree one and, respectively, two, namely where , Q, T 0 , T 1 and T 2 are matrix-valued functions of x and t. The extension to higher degree polynomials results only in an increased computational effort. Moreover, before proceeding further, few preliminary observations and technicalities are required. First, we assume that the N × N matrix , see (12), be constant and Hermitian. Therefore, without any loss of generality, is set to be diagonal and real, namely, in block-diagonal notation, where the real eigenvalues α j , j = 1, . . . , L, satisfy α j = α k if j = k, while 1 j is the n j × n j identity matrix where n j is the (algebraic) multiplicity of the eigenvalue α j . Of course, L j=1 n j = N . Note that this matrix induces the splitting of the set of N × N matrices into two subspaces, namely that of block-diagonal matrices and that of block-off-diagonal matrices. Precisely, if M is any N × N matrix, then we adopt the following notation where is the block-diagonal part of M and is its block-off-diagonal part. Consistently with this notation, the entries M jk of an N × N matrix M are meant to be matrices themselves of dimension n j × n k with the implication that the "matrix elements" M jk may not commute with each other (i.e. the subalgebra of block-diagonal matrices, see (15a), is non-commutative). Moreover, one should keep in mind that the off-diagonal entries M jk are generically rectangular. We also note that the product of two generic N × N matrices A and B follows the rules: is neither block-diagonal nor block-off-diagonal. Next the matrix Q(x, t) in (12) is taken to be block-off-diagonal, whose entries are assumed to be functions of x, t only. Its required property is just differentiability up to sufficiently high order, while no relation among its entries is assumed.
The matrix T , see (12), satisfies compatibility condition (3), which entails the following expression of the coefficients with the following comments and definitions. The matrices C j , j = 0, 1, 2, are constant and block-diagonal, C (o) j = 0. In the following, we set C 0 = 0, because this matrix is irrelevant to our purposes. The notation (·) stands for the linear invertible map acting only on the subspace of the block-off-diagonal matrices (15b) according to the following definition and properties which show that also the matrix (M) is block-off-diagonal. Also the maps D j (·), j = 1, 2, act only on block-off-diagonal matrices according to the definitions Finally, the matrices I j , j = 1, 0, are block-diagonal and take the integral expression .
Because of (19), the matrix Q(x, t) evidently satisfies an integro-differential, rather than a partial differential, equation as a consequence of compatibility condition (3). Incidentally, we note that this kind of non-locality generated by Lax pair (2) has been already pointed out (Degasperis 2011), while its solvability via spectral methods remains an open question. Here, however, we focus our attention on local equations only. Therefore, the condition that Q(x, t) evolves according to a partial differential equation is equivalent to the vanishing of the matrices I j , j = 1, 0, which ultimately implies restrictions on the constant matrices C 1 and C 2 , as specified by the following Proposition 3 The matrices I 1 and I 0 identically vanish if, and only if, the blocks of C 1 and C 2 are proportional to the identity matrix, namely Through the rest of the paper, we maintain these locality conditions so that the resulting evolution equation for the matrix Q reads This equation can be linearized around a given solution Q(x, t) by substituting Q with Q + δ Q and by neglecting all nonlinear terms in the variable δ Q. In this way, we obtain the following linear PDE We are now in the position to formulate our next proposition which is the main result of this section.

Sketch of a Proof
The proof of this result is obtained by straight computation, which is rather long and tedious. Therefore, we skip detailing plain steps, but, for the benefit of the reader, we point out here just a few hints on how to go through the calculation which is merely algebraic.
The guiding idea is to eliminate the variable λ while combining two ODEs (7) so as to obtain a PDE with λ-independent coefficients for the matrix F = [ , ]. To this purpose, one starts by rewriting the first of ODEs (7) for the block-diagonal (d) and Next, one considers the time derivative F t by using expression (23), along with the right-hand-side commutator of the second relation in (7), and one replaces the matrix T (λ) with its expression as in (12) and (16). Then, all terms of the form λF, wherever they appear, should be replaced by the right-hand-side expression given in the first equation in (24), in order to eliminate the spectral variable λ. The outcome of this substitution is that all terms containing the matrix (d) do indeed cancel out, provided conditions (20) are satisfied. The remaining terms of the expression of F t can be rearranged as to coincide with those which appear in the righthand-side expression of the linearized PDE (22), by making use of algebraic matrix identities. Few remarks are in order to illustrate these findings. The matrix [see (23) and (6)], which has the same block-off-diagonal structure (15b) of the matrix δ Q, turns out to be a λ-dependent solution of linearized Eq. (22). Its λ-dependence originates possibly from the (still arbitrary) matrix M and certainly from the matrix solution of Lax pair (2). Indeed, F plays the same role as the exponential solution of linear equations with constant coefficients; namely, by varying λ over a spectrum (see below), it provides the set of the "Fourier-like" modes of linear PDE (22). In this respect, we note that the property of the matrix F of satisfying linearized PDE (22) does not depend on the boundary values of Q(x, t) at x = ±∞. In other words, this result applies as well to solutions Q(x, t) of integrable PDE (21) with vanishing and non-vanishing boundary values, or to periodic solutions, as required in various physical applications. This statement follows from the fact that Proposition 4 has been obtained by algebra and differentiation only. It is, however, clear that the boundary conditions affect the expression of the matrix F through its definition (25) in terms of the solution of the Lax pair of ODEs (2). Moreover, with appropriate reductions imposed on the matrix Q(x, t), integrable matrix PDE (21) associated with Lax pair (2) with (12) includes, as special cases, wave propagation equations of physical relevance. Indeed, the PDE corresponding to L = 2, N = 2 and C 1 = 0 yields the nonlinear Schrödinger equation, while for its multicomponent versions (such as the vector or matrix generalizations of the NLS equation) one may set L = 2, N > 2 and C 1 = 0. By setting instead L = N > 2 and C 2 = 0, one obtains the three-wave resonant interaction system for N = 3 (Kaup 1976b) and many-wave interaction type equations for N > 3 Ablowitz and Segur (1981) (e.g. see Calogero and Degasperis 2004;Degasperis andLombardo 2007, 2009).

x-Spectrum S x of the Solution Q(x, t)
The result stated by Proposition 4 in the previous section implies that any sum and/or integral of F(x, t, λ) over the spectral variable λ is a solution δ Q of linear Eq. (22). Here and in the following, we assume that such linear combination of matrices F(x, t, λ), which formally looks as the Fourier-like integral representation is such to yield a solution δ Q of (22) which is both bounded and localized in the x variable at any fixed time t (periodic perturbations δ Q(x, t) are not considered here). The boundedness of δ Q implies that the matrix F(x, t, λ) be itself bounded on the entire x-axis for any fixed value of t and for any of the values of λ which appear in integral (26). This boundedness condition of F(x, t, λ) defines a special subset of the complex λ-plane, over which the integration runs, which will be referred to as the x-spectrum S x of the solution Q(x, t). This spectrum obviously depends on the behaviour of the matrix Q(x) for large |x|. Indeed, if Q(x, t) vanishes sufficiently fast as |x| → ∞ (like, for instance, when its entries are in L 1 ), then S x coincides with the spectrum of the differential operator (2)). However, if instead Q(x, t) goes to a non-vanishing and finite value as |x| → ∞, this being the case for continuous waves, the spectrum S x associated with the solution Q(x, t) does not coincide generically with the spectrum of the differential operator d/dx − i λ − Q(x). In the present matrix formalism, this happens for N > 2, and it is due to the fact that the spectral analysis applies to the ODE x = [X, ], see (7), rather than to the Lax equation We illustrate this feature for N = 3 in the next section, where we consider the stability of CW solutions. The spectrum S x consists of a piecewise continuous curve and possibly of a finite number of isolated points. Note that at any point λ of the spectrum S x the matrices F(x, t, λ), for any fixed t, span a functional space of matrices whose dimension depends on λ. Here, we do not take up the issue of the completeness and closure of the set of matrices F(x, t, λ) and we rather consider solutions of (22) which have integral representation (26) and vanish sufficiently fast as |x| → ∞, this being a requirement on the matrix M(λ) which appears in definition (23) with (6). As in the standard linear stability analysis, the given solution Q(x, t) is then linearly stable if any initially small change δ Q(x, t 0 ) remains small as time grows, say t > t 0 . Thus, the basic ingredient of our stability analysis is the time dependence of the matrices F(x, t, λ) for any λ ∈ S x . As we show below, getting this information requires knowing the spectrum S x whose computation therefore is our main goal. two coupled NLS equations), the rest of the paper will focus on the matrix case N = 3, Here and below, the asterisk denotes complex conjugation and the four field variables u 1 , u 2 , v 1 , v 2 are considered as independent functions of x and t and are conveniently arranged as two two-dimensional vectors, that is While we conveniently stick to this notation in this section, these formulae will be eventually reduced to those which lead to two coupled NLS equations as discussed in Sect. 3.3. Indeed, we note that adopting non-reduced formalism (28), as we do it here, makes this presentation somehow simpler. In the present case, all matrices are 3 × 3, and T -matrix (16), with C 2 = 2i , C 1 = C 0 = 0, specializes to the expression Then, matrix PDE (21) becomes which is equivalent to the two vector PDEs Here, the dagger notation denotes the Hermitian conjugation (which takes column vectors into row vectors). In this simpler setting, if Q(x, t) is a given solution of Eq.
Moreover, Proposition 4 guarantees that the matrix F(x, t, λ) satisfies this same linear PDE, namely and, for λ ∈ S x , these solutions should be considered as eigenmodes of the linearized equation.
The spectral analysis based on Proposition 4 applies to a large class of solutions Q(x, t) of nonlinear wave equation (31). However, analytic computations are achievable if the fundamental matrix solution (x, t, λ) of the Lax pair corresponding to the solution Q(x, t) is explicitly known. Examples of explicit solutions of the Lax pair are known for particular Q(x, t), for instance multisoliton (reflectionless) solutions (see the stability analysis in Kapitula 2007), cnoidal waves and continuous waves. Here, we devote our attention to the stability of the periodic, continuous wave (CW) solution of (31), or of equivalent vector system (32), In these expressions a and b are arbitrary, constant and, with no loss of generality, real 2-dim vectors: The interest in system (32), or rather in its reduced version (see also Sect. 3.3), is motivated by both its wide applicability and by the fact that its NLS one-component version, for u 2 = v 2 = 0, v 1 = −u 1 , turns out to be a good model of the Benjamin-Feir (or modulational) instability which is of great physical relevance (Benjamin and Feir 1967;Hasimoto and Ono 1972). This kind of instability of the plane wave solution of the NLS equation occurs only in the self-focusing regime. In contrast, in the coupled NLS equations the cross-interaction and the counterpropagation (q = 0) introduce additional features (e.g. see Forest et al. 2000; Baronio et al. 2014) which have no analogues in the NLS equation. This is so if the instability occurs even when the selfinteraction terms have defocusing effects (see Sect. 3.3). Moreover, and by comparing the coupled case with the single field as in the NLS equation, we observe that this CW solution (35) depends on the real amplitudes a 1 , a 2 , b 1 , b 2 and, in a crucial manner (see below), on the real parameter q which measures the wave number mismatch between the two wave components u 1 and u 2 (or v 1 and v 2 ). The main focus of this section is understanding how the spectrum S x changes by varying the parameters a 1 , a 2 , b 1 , b 2 and q. The results we obtain here will be specialized to the CNLS system in Sect. 3.3. In matrix notation, see (28), this continuous wave solution (35) reads where the matrix has expression (27), while the matrix σ is and we conveniently introduce the real parameters which will be handy in the following. Next we observe that a fundamental matrix solution (x, t, λ) of the Lax equations has the expression where the x, t-independent matrices W and Z are found to be with the property that they commute, [W , Z ] = 0, consistently with compatibility condition (3). In order to proceed with plain arguments, we consider here the eigenvalues w j (λ) and z j (λ), j = 1, 2, 3, of W (λ) and, respectively, of Z (λ) as simple, as indeed they are for generic values of λ. In this case, both W (λ) and Z (λ) are diagonalized by the same matrix U (λ), namely Next we construct the matrix F(x, t, λ) via its definition, see (23) and (6), which, because of explicit expression (40), reads As for the matrix M(λ), it lies in a nine-dimensional linear space whose standard basis is given by the matrices B ( jm) , whose entries are where δ jk is the Kronecker symbol (δ jk = 1 if j = k and δ jk = 0 otherwise). However, the alternative basis V ( jm) , which is obtained via the similarity transformation where U (λ) diagonalizes W and Z (see (43)), is more convenient to our purpose. Indeed, expanding the generic matrix M(λ) in this basis as the scalar functions μ jm being its components, and inserting this decomposition into expression (45), leads to the following representation of F where we have introduced the x, t-independent matrices The advantage of expression (49) is to explicitly show the dependence of the matrix F on the six exponentials e i [(x(w j Proposition 4 stated in the previous section guarantees that, for any choice of the functions μ jm (λ), expression (49) be a solution of linearized Eq. (33), see (34). It is plain (see (26)) that the requirement that such solution δ Q(x, t) be localized in the variable x implies the necessary condition that the functions μ jm (λ) be vanishing for j = m, μ j j = 0, j = 1, 2, 3. The further condition that the solution δ Q(x, t) be bounded in x at any fixed time t results in integrating expression (49) with respect to the variable λ over the spectral curve S x of the complex λ-plane (see (26)): Here, according to Sect. 2.1, S x can be geometrically defined as follows: Definition 1 The x-spectrum S x , namely the spectral curve on the complex λ-plane, is the set of values of the spectral variable λ such that at least one of the three complex numbers k j = w j+1 − w j+2 , j = 1, 2, 3 (mod 3), or explicitly is real.
Observe that the k j 's play the role of eigenmode wave numbers (see (49)).
To the purpose of establishing the stability properties of continuous wave solution (35), we do not need to compute integral representation (51) of the solution δ Q of (33). Indeed, it is sufficient to compute the eigenfrequencies as suggested by the exponentials which appear in (49). Their expression follows from matrix relation (42) and read This expression looks even simpler by using the relation w 1 + w 2 + w 3 = −λ implied by the trace of the matrix W (λ) (see (41)) and finally reads The consequence of this expression (56), which is relevant to our stability analysis, is given by the following (35) is stable against perturbations δ Q whose representation (26) is given by an integral which runs only over those values of λ ∈ S x which are strictly real.

Proposition 5 Continuous wave solution
Proof If the spectral variable λ is real, then all coefficients of the characteristic polynomial of the matrix W (λ) (41), are real. Indeed, these coefficients depend on the real parameters of CW solution (35), namely the wave number mismatch q, and the parameters p and r defined by (39). Therefore, the roots w j of P W (w; λ) are either all real, or one real and two complex conjugate. In the first case, three wave numbers k j (52) are all real, and thus, the corresponding λ is in the spectrum S x , λ ∈ S x . In the second case, the three k j are all complex, i.e. with non-vanishing imaginary part, and thus λ lies in a forbidden interval of the real axis which does not belong to the spectrum S x , λ / ∈ S x . In the following, we refer to one such forbidden real interval as a gap, see Sect. 3.1. Consequently, in the first case, since λ, w 1 , w 2 , w 3 and therefore, k 1 , k 2 , k 3 , are all real, then also ω 1 , ω 2 , ω 3 , see (56), are all real, with the implication of stability. Indeed, all eigenmode matrices (49) remain small at all times if they are so at the initial time.
On the other hand, let us assume now that λ is complex, λ = μ + i ρ, with nonvanishing imaginary part, ρ = 0. Let w j = α j + i β j be the (generically complex) roots of P W (w; λ). With this notation, we have that one of the wave numbers, say k 3 , will be real only if β 1 = β 2 = β. Then, from (56), we have that ω 3 will also be real only if β 3 = ρ. Writing the polynomial P W (w; λ) as P W (w; λ) = 3 j=1 (w − α j − i β j ), and comparing the real and imaginary parts of the coefficients of same powers of w from this expression with those obtained from (57), we get a system of six polynomial equations for the six unknowns α 1 , α 2 , α 3 , β, μ and ρ, each equation being homogeneous and of degree 1, 2 or 3 in the unknowns: Then, it is immediate to show by means of elementary algebraic manipulations that the above system has real solutions for α 1 , α 2 , α 3 , β, μ and ρ, only if ρ = −β = 0 (unless the non-physical and non-generic condition p = r = 0 be met). This contradicts the original assumption ρ = 0. Therefore, for a representation (26) of a perturbation δ Q to be bounded (and thus for the corresponding CW solution to be stable), the integral in (26) must run only on those values of λ ∈ S x which are strictly real.
This Proposition 5 implies that a real part of the spectrum S x is always present, and this part may or may not have gaps (see the next Sect. 3.1). On the contrary, as it will be proved (see Sect. 3.2), a complex component of the spectrum, namely one which lies off the real axis of the λ-plane, may occur and it always leads to instabilities. This important part of the spectrum S x is made of open curves, which will be referred to as branches, and/or closed curves of the complex λ-plane which will be termed loops.
Before proceeding further, we observe that the effect of the cross-interaction is crucially influent on the stability only if the phases of two continuous waves (35) have different x-dependence, namely if q = 0 (see (35)). Indeed, if q = 0 the stability property of the continuous wave is essentially that of the NLS equation. This conclusion results from the explicit formulae w 1 = λ 2 − p, w 2 = − λ 2 − p, w 3 = −λ, and therefore which show that if p > 0 (see (39a)), the spectrum S x is the real axis with the exclusion of the gap {− √ p < λ < √ p }, while, if p < 0, the spectrum S x is the entire real axis with the addition of the imaginary branch λ = iρ, {− √ − p < ρ < √ − p }. Therefore, the stability for p > 0 follows from the reality of the frequencies ω 1 , ω 2 , ω 3 over the whole spectrum S x , see (58). The modulational instability occurs only if p < 0 since, for λ in the imaginary branch of the spectrum S x , the frequencies ω 1 , ω 2 and ω 3 have a non-vanishing imaginary part.
In the following, while computing the spectrum S x , we consider therefore only the case q = 0 and, with no loss of generality, strictly positive, q > 0. In this respect, we also note that the assumption that q be non-vanishing makes it possible, without any loss of generality, to rescale it to unit, q = 1, while keeping in mind that, by doing so, λ, w j rescale as q, while p, r , z j rescale as q 2 , i.e.
However, we deem it helpful to the reader to maintain the parameter q in our formulae, here and in the following, to have a good control of the limit as q vanishes. In addition, here and thereafter, our computations and results will be formulated for non-negative values of the parameter r according to the following Proposition 6 With no loss of generality, the relevant parameter space reduces to the half (r , p)-plane with r ≥ 0.
Indeed, the change of sign λ → −λ, r → −r takes the characteristic polynomial P W (w; λ) (57) into −P W (−w; λ) and this implies that our attention may be confined to non-negative values of r only.

Gaps
Although computing the position of the gaps of S x on the real λ-axis is not strictly relevant to the issue of stability because of Proposition 5, this information is essential to the construction of soliton solutions corresponding to real discrete eigenvalues which lie inside a gap. With this motivation in mind, we devote this subsection to this task. Using the Cardano formulae for finding workable explicit expressions of the three roots w j (λ) of P W (w; λ) for each value of p, r and q is not practical. To overcome this difficulty, we will adopt here an algebraic-geometric approach.
We begin by computing the discriminant D W = w P W (w; λ) of the polynomial P W (w; λ) (57) with respect to w, obtaining thus a polynomial in λ, with parameters q, p, r , which is positive whenever the three roots w j are real, and it is negative if instead only one root is real. Consequently, for any given value of the parameters q, p and r , the discriminant is negative, D W < 0, for those real values of λ which belong to a gap of S x , while its zeros are the end points of gaps. Here and hereafter, the notation y P(x, y, z) stands for the discriminant of the polynomial P with respect to its variable y.
In the exceptional case r = 0, computing the roots of the discriminant D W as a polynomial in the variable λ reduces to the factorization of a quadratic polynomial, allowing the formulation of the following Proposition 7 For r = 0, the gaps depend on the parameter p as follows: n og a p 0 ≤ p < q 2 two gaps, −g + ≤ λ ≤ −g − and g − ≤ λ ≤ g + q 2 ≤ p < +∞ one gap, −g + ≤ λ ≤ g + The border values g ± are The threshold values p = 0 and p = q 2 correspond to the opening and the coalescence of two gaps, respectively.
In the generic case r ≥ 0, we find that the number of gaps of the real part of the spectrum S x is either zero, or one, or two. Gaps appear at double zeros of the discriminant D W (λ; q, p, r ), thus at zeros of its own discriminant λ D W (λ; q, p, r ) = λ w P W (w; λ), namely when (61) The three polynomial factors appearing in (61) bound the regions of the (r, p)-plane characterized by different numbers of gaps.
By varying the values of the parameters q, p and r , we expect the double zeros to open and form the gaps, namely intervals where D W < 0. In order to find the number of gaps from expression (59) of the discriminant D W (λ; q, p, r ), it is convenient to plot the dependence of the parameter p = p(λ) on the spectral variable λ at the zeros of this discriminant for fixed values of the parameters r and q. This function p(λ) is therefore implicitly defined by the equation D W (λ; q, p(λ), r ) = 0. For a fixed value of q, in the (λ, p)-plane the p(λ) curve has always one cusp (a double singular point), and either two minima, or one minimum and one maximum, depending on the value of r (see "Appendix A" for details). The position of the cusp (λ S , p S ) as a function of r and q is Furthermore, p S (r ) seen as an algebraic curve in the (r, p)-plane satisfies the following implicit relation The transition between the two regimes (two minima, or one minimum and one maximum) takes place at the special threshold value r T = 4q 2 where p S (r T ) = r T (see "Appendix A"). If 0 ≤ r < r T , the cusp is also a local maximum. Equipped with these findings, we formulate Proposition 8 For real λ and r > 0, S x has the following gap structure.
The threshold values p = −r, p = r and p = p S (r ) correspond to the opening of one gap, the opening of two gaps and the coalescence of two gaps, respectively.
The threshold values p = −r, p = p S (r ) and p = r correspond to the opening of one gap, the opening of two gaps and the coalescence of two gaps, respectively.
As implied by (61) and Proposition 8, we identify three threshold curves in the (r, p)-half-plane; they are p = p ± (r ) = ±r (64) and the curve p = p S (r ) which is explicitly expressed by (62) and plotted in Fig. 1. Note that, according to Proposition 6, we will focus only on the half-plane r ≥ 0. These curves will play a role also in the next subsection where we give a complete classification of branches and loops of the spectrum.

Branches and Loops
In the previous subsection, we have described the gaps of the spectrum S x , namely those values of the spectral parameter λ which are real but do not belong to the spectrum. Here, we face instead the problem of finding the subset of the λ-plane which is off the real axis but belongs to S x . This subset is made of smooth and generically finite, open (branches) or closed (loops), continuous curves (with the possible exception of the case r = 0, see Fig. 3h and Degasperis et al. 2018). By Definition 1, λ ∈ S x if at least one of the corresponding three wave numbers k j (λ), j = 1, 2, 3 is real, see (52). Going back to the characteristic polynomial P W (w; λ) (57), we look now at its roots w j , and at the wave numbers k j , j = 1, 2, 3, see (52). If λ is not real, the roots w j cannot be all real themselves since the coefficients of P W (w; λ) are not real. As a consequence, the requirement that one of the wave numbers (52) be real implies that at least two roots of P W (w; λ), for instance w 1 and w 2 , have the same imaginary part. In Fig. 1 The three threshold curves p ± (r ) = ± r (solid black) and p S (r ) (solid red) for q = 1, as parametrically defined by (64) and (63) or explicitly by (62), are plotted. They are boundaries of regions of the (r, p)-plane, r > 0, where the number of gaps, either 0 or 1, or 2, is shown. It is also shown that p S (r ) > r for r < 4 and that p S (r ) < r for r > 4 (see Proposition 8) (Color figure  online) fact, branches and loops of S x may exist only in this case: indeed, if the three roots w j have all the same imaginary part, it can be shown that no smooth branch or loop exists. Hence, here and below, we name k 3 = w 1 − w 2 the only real wave number, while k 1 and k 2 are complex (note that k 1 + k 2 + k 3 = 0). The complex part of the spectrum S x is therefore defined as the set of the λ-plane such that Im(k 3 (λ)) = 0. In order to compute this component of the spectrum, we introduce the novel polynomial P(ζ ), defined by the requirement that its roots are the squares of the wave numbers, ζ j = k 2 j . It is plain that its coefficients d j are completely symmetric functions of the roots w j of the characteristic polynomial P W (w; λ). This is evidently so because d 1 , d 2 , d 3 are symmetric functions of ζ 1 = (w 2 − w 3 ) 2 , ζ 2 = (w 3 − w 1 ) 2 , ζ 3 = (w 1 − w 2 ) 2 , and therefore of w 1 , w 2 , w 3 , with the implication that the coefficients d 1 , d 2 , d 3 of P(ζ ) are polynomial functions of the coefficients c 1 , c 2 , c 3 of characteristic polynomial (57), P W (w) = w 3 + c 1 w 2 + c 2 w + c 3 , These relations, which read d 1 = 2(3c 2 − c 2 1 ) , combined with the expressions of the coefficients c j , see (66), yield the dependence of the coefficients d j of P(ζ ) on the parameters r , p, q and on the spectral variable λ, so that P(ζ ) ≡ P(ζ ; λ, q, p, r ). This dependence turns out to be where D W is the discriminant of P W , see its expression (59), with the implication that P(ζ ; λ, q, p, r ) is of degree 4 in the spectral variable λ. As explained at the beginning of this section, we recall that, for an arbitrary complex λ, at most one wave number (named k 3 ) is real. Therefore, our task of computing the complex part of the spectrum S x , for a given value of the parameters r , p, reduces to finding the curve in the λ-plane along which the root ζ 3 (λ) = k 2 3 (λ) of P(ζ ; λ, q, p, r ) remains real and positive, ζ 3 = k 2 3 > 0. However, it is in general much more convenient, from both the analytical and computational points of view (and indeed this is what we will be doing in the following) to reverse the perspective and regard the x-spectrum S x as the locus in the λ-plane of the λ-roots of P(ζ ; λ, q, p, r ) seen as a real polynomial in λ, for ζ spanning over the semi-line [0, +∞). In other words, for a fixed value of ζ ≥ 0, we compute the four λ-roots of the real polynomial P(ζ ; λ, q, p, r ), which are the values of λ (irrespective of being complex or purely real) such that the wave number k 3 = √ ζ is real. In this way, the two cases of Im(λ) = 0 and Im(λ) = 0 can be treated at once (allowing to retrieve and confirm all the results about gaps in the spectra presented in the previous subsection).
From the algebraic-geometric point of view, the locus of the roots of P(ζ ; λ, q, p, r ) in the λ-plane is an algebraic curve; this can be given implicitly as a system of two polynomial equations in two unknowns by setting λ = μ + i ρ and then separating the real and the imaginary parts of P(ζ ; λ, q, p, r ).
In order to explain how the spectrum changes over the parameter space, we apply Sturm's chains (e.g. see Demidovich andMaron 1981, Hook andMcAree 1990) to the polynomial Q(ζ ; q, p, r ) = λ P(ζ ; λ, q, p, r ), that is, to the discriminant of the polynomial P(ζ ; λ, q, p, r ) with respect to λ. Indeed, it turns out that the nature and the number of the components of the x-spectrum S x , seen as a curve in the λ-plane, are classified by the number of sign changes of Q(ζ ; q, p, r ), for ζ ≥ 0. This procedure requires some technical digression; thus, in order to avoid breaking the narrative here, we refer the reader to "Appendix B" for all details.
Our findings provide a complete classification of the spectra over the entire parameter space, i.e over the (r, p)-plane, r ≥ 0. First we note that the number of branches (B) can only be 0, 1, or 2, while there may occur either 0 or 1 loop (L). Moreover, in addition to the three threshold curves p ± (r ) (64) and p S (r ) (62) introduced in Sect. 3.1, one more threshold curve, p = p C (r ), is found, which, for a given non-vanishing value of q, is implicitly defined as or, explicitly, as (a) (b) Fig. 2 The structure of the spectrum S x in the (r, p) plane, q = 1 (see Proposition 9). In red, the curve p S ; in blue, the curve p C ; in grey, the curves p = ± r . a r ≤ 4q 2 . b r ≥ 4q 2 (Color figure online) It is now convenient to combine these results with those we obtained in the previous subsection on the gaps (G) on the real λ-axis of the complex λ-plane, obtaining a complete classification of the spectra. We find that only five different types of spectra exist according to different combinations of gaps, branches and loops. These are the following: where the notation nX stands for n components of the type X, with X either G, or B, or L. In Fig. 2, the four threshold curves p ± , p S , p C (with q = 1) are plotted to show the partition of the parameter half-plane according to the gap, branch and loop components of the spectrum, for both r < r T = 4 q 2 and r > r T = 4 q 2 . This overall structure of the spectrum S x in the parameter space can be summarized by the following Proposition 9 For r ≥ 0, the spectrum S x has the following structure in the parameter plane (r, p).
At the threshold values, the dynamics of the transitions between the five types of spectra can be very rich, and phenomena such as the merging of a branch and a loop to form a single branch can be observed; however, in the simplest cases (which are the majority), gaps, branches and loops open up from, or coalesce into points. An important threshold case is p = r (entailing a 2 = 0, and as such non-generic): this is the only choice of the parameters p, and r for which the spectrum is entirely real (and thus the CW solution is stable), with one gap in the interval − 2 . Examples of different spectra have been computed numerically by calculating the zeros of P(ζ ; λ, q, p, r ) as a polynomial in λ (see "Appendix C" for details) and are displayed in Fig. 3. They are representative of the five types and correspond to various points of the parameter (r, p)-space as reported in Proposition 9.
Finally, the limit case r = 0 of the parameter half-plane r > 0 deserves a special mention. In fact, four different limit spectra are found for r = 0, namely in the four intervals −∞ < p < −4 q 2 , −4 q 2 < p < 0, 0 < p < q 2 and q 2 < p < +∞. After meticulous analysis (see Degasperis et al. 2018), we find that these limit spectra are consistent with those depicted by Proposition 9. However, they have to be carefully dealt with as, in this limit, two gaps may coalesce in one, or a loop may go through the point at infinity of the λ-plane, thereby covering the whole imaginary axis. These spectra are displayed in Fig. 3f-h. We should also point out that the condition r = 0 allows for analytic computations, and explicit formulae, since, as already noted in Sect. 3.1 for the discriminant D W (λ; q, p, r ), the coefficients d j (λ) of the relevant polynomial P(ζ ; λ, q, p, r ) (65) reduce to polynomials of degree 2 in the variable λ 2 .

Modulational Instability of Two Coupled NLS Equations
In this section, we discuss the consequences of the results obtained so far in the framework of the study of the instabilities of a system of two coupled NLS equations.
In the scalar case, the focusing NLS equation has played an important role in modelling modulational instability of continuous waves (Benjamin and Feir 1967;Hasimoto and Ono 1972). This unstable behaviour is predicted for this equation by simple arguments and calculations. It is therefore rather remarkable that, on the contrary, a nonlinear coupling of two NLS equations makes the unstable dynamics of two interacting continuous waves fairly richer than that of s 1 = s 2 = − 1, a 1 = 1, a 2 = 4, as an example of a 0G 2B 0L spectrum. b S x for r = 1, p = − 1.5, i.e. s 1 = s 2 = − 1, a 1 = 0.5, a 2 = 1.118, as an example of a 0G 2B 1L spectrum. c S x for r = 1, p = 2.8, i.e. s 1 = s 2 = 1, a 1 = 1.3784, a 2 = 0.94868, as an example of a 1G 1B 0L spectrum. d S x for r = 15, p = − 13.45, i.e. s 1 = 1, s 2 = − 1, a 1 = 0.88034, a 2 = 3.7716, as an example of a 1G 1B 1L spectrum. e S x for r = 1, p = 1.025, i.e. s 1 = s 2 = 1, a 1 = 1.0062, a 2 = 0.1118, as an example of a 2G 0B 1L spectrum. f S x for r = 0, p = 0.1, i.e. s 1 = s 2 = 1, a 1 = a 2 = 0.22361, as an example of a degenerate case of a 2G 0B 1L spectrum: the imaginary axis in the spectrum is a loop passing through the point at infinity. g S x for r = 0, p = − 14, i.e. s 1 = s 2 = − 1, a 1 = a 2 = 2.6458, as an example of a degenerate case of a 0G 2B 0L spectrum: both branches are entirely contained on the imaginary axis; one branch passes through the point at infinity, whereas the other branch passes through the origin; for r = 0, two symmetrical gaps open on the imaginary axis for p ≤ − 8, as explained in Degasperis et al. (2018). h S x for r = 0, p = − 4.7, i.e. s 1 = s 2 = − 1, a 1 = a 2 = 1.533, as an example of a degenerate case of a 0G 2B 0L spectrum: this case (which also appears in Ling and Zhao 2017) when projected back onto the stereographic sphere, can be completely explained in terms of the classification scheme provided in Proposition 9 (see Degasperis et al. 2018) (Color figure online) a single wave. Since the method of investigating the stability of two coupled NLS equations that we have presented in the previous sections requires integrability, we deal here with integrable system (1), also known as generalized Manakov model, that we recall here for convenience: This system is obtained by setting v = Su in (32), where The system of two coupled NLS equations is of interest in various physical contexts and the investigation of the stability of its solutions deserves special attention. Once a solution u 1 (x, t), u 2 (x, t) has been fixed, linearized equation (33) around this solution are The two coupling constants s 1 , s 2 , if non-vanishing, are just signs, s 2 1 = s 2 2 = 1, with no loss of generality. Thus, CNLS system (1) models three different processes, according to the defocusing or focusing self-and cross-interactions that each wave experiences. These different cases are referred to as (D/D) if s 1 = s 2 = 1, as (F/F) if s 1 = s 2 = −1 and as (D/F) in the mixed case s 1 s 2 = −1.
Hereafter, our focus is on the stability of the CW solution (see (35), (36) with b j = s j a j ) where the parameter q is the relative wave number, which may be taken to be nonnegative q ≥ 0, and a 1 , a 2 are the two amplitudes, whose values, with no loss of generality, may be real and non-negative, a j ≥ 0. As for the notation, the parameter p is defined by (39a) which, in the present reduction, reads To make contact with the stability analysis presented in the two preceding sections, we introduce also the second relevant parameter r (39b) whose expression, in the present context, is r = s 1 a 2 1 − s 2 a 2 2 .
At this point, we show in Fig. 4 the (r, p)-plane as divided into octants, according to different values of amplitudes and coupling constants. We note that, as pointed out by Proposition 6, it is sufficient to confine our discussion to the half-plane r ≥ 0, where four different experimental settings are represented, one in each of these four octants. These four states of two periodic waves can be identified either by the parameters p and r or by amplitudes and coupling constants according to the relations s 1 a 2 1 = each octant being characterized as shown in Fig. 4. In particular, the octant r > p > 0 will be referred to as (D > F), whereas the octant −r < p < 0 as (D < F), to indicate the wave with larger amplitude. Moreover, we observe again that setting q = 0, which characterizes two standing waves with equal wavelength, makes the stability properties of the CW solution very similar to that of the NLS equation. Indeed, in this case solution (74) is stable if both waves propagate in a defocusing medium, s 1 = s 2 = 1, and is unstable in the opposite case s 1 = s 2 = −1. However, this solution remains stable even if the wave u 2 feels a self-focusing effect, i.e. s 2 = −1, provided the other wave u 1 in the defocusing medium, s 1 = 1, has larger amplitude a 1 > a 2 . Similarly, when q = 0, the CW solution is unstable if s 1 = 1 and s 2 = −1, and if the larger amplitude wave is the one which propagates in the focusing medium, a 2 > a 1 . These remarks follow from explicit formulae (58) and (75a), and they are clearly evident in Fig. 4 where, for r > 0, the upper two octants correspond to p > 0 and the lower two to p < 0. The limit q → 0 of the results presented in the previous sections should be considered as formally singular. Indeed, the behaviour of CW (74) against small generic perturbations becomes significantly different from that reported above if the wave number mismatch 2q is non-vanishing. For one aspect, this is apparent from our first Remark 1 If q = 0, then for no value of the amplitudes a 1 , a 2 with a 1 a 2 = 0, and for no value of the coupling constants s 1 , s 2 , the spectrum S x is entirely real.
This follows from Proposition 9 and from our classification of the spectra in the parameter space into five types, all of which have a non-empty complex (non-real) component. The condition a 1 a 2 = 0 for this statement to be true coincides with the condition that the point (r, p) neither belongs to the threshold curve p = p + (r ), nor to the threshold curve p = p − (r ), see (75). In particular, in the case p = r > 0, which is equivalent to setting a 2 = 0 and s 1 = 1, the roots of P W (w; λ), (57), can be given in closed form, with the implication that λ has to be real to be in S x . In this case, the spectrum is not immediately identifiable as one of the five types of spectra described in Proposition 9, for it corresponds to a threshold case between two of such types, and features one gap, no branches and no loops. Moreover, this spectrum coincides with that obtained for the defocusing NLS equation via a Galilei transformation. In a similar fashion, if the point (r, p) is taken on the curve p = p − (r ), namely p = −r < 0 or, equivalently, a 1 = 0, s 2 = −1, the corresponding spectrum is that of the CW solution of the focusing NLS equation modulo a Galilei transformation, namely with no gap on the real axis, one branch on the imaginary axis and no loop. Remark 1 has a straight implication on the stability of solution (74). This can be formulated as Remark 2 If q a 1 a 2 = 0, then CW solution (74) is unstable.
The relevant point here is the dependence on λ of the frequencies ω j (λ) over the spectrum S x . As we have shown in the previous section, if q a 1 a 2 = 0, the spectrum S x consists of two components: one, RS x , is the real axis Im(λ) = 0, with possibly one or two gaps (see Proposition 8), and the second one, CS x , consists of branches and/or loops where λ runs off the real axis (see Proposition 9). Therefore, spectral representation (51) of the perturbation takes the form which is meant to separate the contribution to δu due to the real part of the spectrum from that coming instead from the complex values of the integration variable λ, this being the integration over branches and loops. The 2-dim vector functions f ( j) ± (λ) do not play a role here and are not specified. On the contrary, the reality property of the frequencies ω j over the spectrum is obviously essential to stability. Since the reality of the frequencies ω 1 , ω 2 , ω 3 for λ real has been proved in Proposition 5, it remains to show that indeed ω 3 (λ) is not real if λ belongs to branches or loops. This follows from explicit expression (56), namely where k 3 is real but w 3 and λ cannot be real and cannot have same imaginary part as proved in Proposition 5, (see Sect. 3.2). Here, the imaginary part of ω 3 defines the gain function over the spectrum. This (possibly multivalued) function of k 3 plays an important role in the initial stage of the unstable dynamics. Precisely, its dependence on the wave number k 3 gives important information on the instability band and on timescales (Degasperis et al. 2018).
The physical interpretation is more transparent if these considerations are stated in terms of the amplitudes of CW solution (74) for the three choices of the coupling constants s 1 , s 2 corresponding to integrable cases. In this respect, it is convenient to translate the classification of the spectra in the (a 1 , a 2 )-plane, in particular, and with no loss of generality, in the quadrant a 1 ≥ 0, a 2 ≥ 0. The four threshold curves p ± (r ), p S (r ), p C (r ) which have been found in the (r, p)-plane, see Fig. 2, are reproduced below in the (a 1 , a 2 )-plane, see Proposition 10 For each of the three integrable cases of two coupled NLS equations (1), the spectrum is of the following types.
(D/D) s 1 = s 2 = 1 (see Fig. 5a) The lower octant (a 1 ≥ a 2 ), which corresponds to r ≥ 0, gets divided into two parts by the threshold curve p S which intersects the line a 2 = a 1 at the point (1/ √ 2, 1/ √ 2). In the lower, finite, part of this octant, the spectrum S x is of type 2G 0B 1L, while in the other part it is of type 1G 1B 0L. (F/F) s 1 = s 2 = −1 (see Fig. 5b) The upper octant (a 2 ≥ a 1 ), which corresponds to r ≥ 0, gets divided into two parts by the threshold curve p C which intersects the line a 2 = a 1 at the point ( √ 2, √ 2). In the lower, finite, part of this octant, the spectrum S x is of type 0G 2B 1L, while in the other part it is of type 0G 2B 0L. (D/F) s 1 = 1 , s 2 = −1 (see Fig. 5c) The whole quadrant corresponds to r ≥ 0. It is divided into three infinite portions by the two threshold curves p C , in the upper octant, and p S in the lower one. The spectrum is of type 1G 1B 1L in the upper part, of type 1G 1B 0L in the middle part and of type 2G 0B 1L in the lower part.
These statements are straight consequences of Proposition 9 via transformation (76). We conclude by mentioning the limit case r = 0, see the end of the previous section. This concerns only the (D/D) and (F/F) cases of Proposition 9, and it coincides with the limit a 1 = a 2 , namely with the case in which the two wave amplitudes are strictly equal. On this particular line of the (a 1 , a 2 )-plane, i.e. r = 0 in the (r, p)-plane, the spectra are of four types, according to numbers of gaps, branches and loops, as it has been shortly discussed in Sect. 3.2.

Summary and Conclusions
A sufficiently small perturbation of a solution of a (possibly multicomponent) wave equation satisfies a linear equation. If the wave equation is integrable, the solution of this linear equation is formulated in terms of a set of eigenmodes whose expression is explicitly related to the solutions of the Lax pair. We give this connection in a general N × N matrix formalism, which is local (in x) and does not require specifying the boundary condition nor the machinery of the direct and inverse spectral problem. A by-product of our approach is the definition of the spectrum S x , associated with the unperturbed solution of the wave equation. It is worth stressing that the spectrum S x for N > 2 does not coincide with the spectrum of the Lax equation x = (iλ + Q) in the complex λ-plane. This result is explicitly shown for N = 3 in the instance of continuous wave solutions of two CNLS equations. In this case, we define as functions of λ on the spectrum S x the set of wave numbers k j (λ) and frequencies ω j (λ) with the implication that the dispersion relation is given in parametric form, the parameter being the spectral variable λ which appears in the Lax pair. In general, the spectrum is a complicated piecewise continuous curve. It obviously changes in the parameter space which is the set of values of the amplitudes of the two CWs, the mismatch q of their wave numbers, and the values of the coupling constants s 1 , s 2 . Apart from particular values of the parameters, the computation of the spectrum is not achievable in analytic form and has to be done numerically. The knowledge of the spectrum is sufficient to assess the stability of the CW solution. Generally, the spectrum consists of the real λ-axis with possibly one or two forbidden bands (gaps) and few additional finite curves which may be open (branches) or closed (loops). According to these topological properties, spectra can be classified in five different types to completely cover the entire parameter space. Only few marginal cases require separate consideration. Physically relevant information comes from the λ dependence of the eigenfrequency on branches and loops. In particular, one can read out of this dependence on λ the instability band, whether at large (as in the Benjamin-Feir instability) or at small wavelengths, and the time scale of the exponential growth in time of the perturbation. This is characterized by the imaginary part of the complex frequency, namely by the gain function. However, the investigation of this interesting aspect of the stability analysis is not reported here as it requires further analysis and computations. This part of our work will be reported elsewhere (Degasperis et al. 2018).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A Structure of the Gaps of the Spectrum S x
The aim of this section is to provide additional details of the gap structure of the spectrum S x and to give a proof of Proposition 8 in Sect. 3.1, which states that the real part of the spectrum S x has either one, two or no gaps in the (r, p)-plane with r ≥ 0. The no gap case occurs if the three eigenvalues w j (λ) are real on the entire λ-axis. This happens if discriminant (59) is positive for any real value of λ. By varying the value of p and/or r a gap may open up at the double zeros of the discriminant D W (λ; q, p(λ), r ). Inside a gap, the discriminant is negative and the three wave numbers k j (λ) (52) have a non-vanishing imaginary part. We define implicitly the function p ≡ p(λ) as those values of the parameter p such that the discriminant D W (λ; q, p(λ), r ) = 0.
In order to compute the number of gaps from expression (59) of the discriminant D W (λ; q, p, r ), we plot p(λ) as a function of λ real, for few fixed values of r , and q = 1. Four such plots of p(λ) are shown in Fig. 6 for four different values of r . The intersections of this plotted curve with the straight line corresponding to a constant value of p provide the endpoints of those gaps which occur at the given value of p (and at the fixed values of r and q). For instance, the analytical results presented above for r = 0, see Proposition 7, can be clearly read out of the left plot in Fig. 6a. The two local minima at λ = ±0.5, p = 0 show the opening of the two gaps that we have analytically found for 0 < p < 1. For future reference, we also note that the local maximum of p(λ) is taken at λ = λ S = 0, namely p(λ S ) = p S = 1, and that the function p(λ) has a cusp (a double singular point), namely is not analytic, at λ S . If the parameter r increases, r > 0, this plot gets asymmetrically deformed as shown in Fig. 6b-d. Figure 6b shows that no gap exists for p < −r = −0.5, while for  λ; q, p(λ), r ) = 0; the grey area corresponds to D W < 0. a r = 0, q = 1. b r = 0.5, q = 1. c r = 4, q = 1. d r = 25, q = 1 (Color figure online) − 0.5 = − r < p < r = 0.5 only one gap occurs. Instead, for 0.5 = r < p < p S two distinct gaps appear. Finally, for p S < p < +∞ again only one gap is present. As in the previous plot, here p S denotes the local maximum and singular point of p(λ). Moreover, we observe that this graph is consistent with the following alternative expression of discriminant (59) for p = −r , which is clearly positive definite for any non-negative r , and any real λ with the exception of λ = r 4q − q 2 where p(λ) takes its minimum value. By increasing further the value of r , Fig. 6b changes and eventually looks different as the local maximum p S becomes a local minimum. The transition is illustrated in Fig. 6c computed at r = 4, with q = 1. After the transition, the curve p(λ) takes the form shown in Fig. 6d computed at r = 25, with q = 1. Again, p = −r is still the minimum of p(λ), while p = r is instead a local maximum. The singularity at λ = λ S is such that now p S = p(λ S ) is a local minimum with p S < r . Indeed, λ S is the only singularity of p(λ) and its dependence on the parameter r , λ S = λ S (r ), is implicitly expressed by its inverse relation r = r (λ S ), This yields a one-to-one correspondence between any non-negative real r and a real λ S (r ) since we prove that, for any r ≥ 0, cubic Eq. (A.1) has two complex conjugate roots. Moreover, we find also the following useful relations associated with the singular value λ S : 2) The transition between the two regimes with p S (r ) > r and p S (r ) < r takes place at the special threshold value r = r T , where p S (r T ) = r T . With self-evident notation, we find that The implicit expressions of λ S (r ) (A.1) is made explicit in (62). The function p S (r ) is plotted in Fig. 1.

B Analysis of the x-Spectrum
As pointed out in Sect. 3.2, the nature and the number of the components of the xspectrum S x , seen as a curve in the λ-plane, are classified by the number of sign changes of the discriminant of P(ζ ; λ, q, p, r ) with respect to λ, i.e. λ P(ζ ; λ, q, p, r ) = Q(ζ ; q, p, r ), for ζ ≥ 0 (see also (B.1a) in the following). The aim of this section is to give details of this result. Without loss of generality, from now on we will assume q = 1. Given the (cubic) polynomial P W (w; λ, q, p, r ), we would like to describe the geometrical locus in the λ-plane constituted by the complex values of λ such that, for each value of λ in the locus, the polynomial P W (w; λ, q, p, r ) has at least two w-roots such that their difference is real. As discussed in Sect. 3.2, since we are interested in the differences of the roots, it is convenient to introduce the polynomial P(ζ ; λ, q, p, r ), see (65), whose ζ -roots are the squares of the differences of the original w-roots of P W (w; λ, q, p, r ). The polynomial P(ζ ; λ, q, p, r ) can be conveniently rewritten as Observe that ζ = k 2 j for some j. As in Sect. 3.2, let j = 3. The polynomial P(ζ ; λ, q, p, r ) can be regarded as a 3 rd degree polynomial in ζ , as well as a 4 th degree polynomial in λ. Hence, our problem translates into finding the complex values of λ for which P(ζ ; λ, q, p, r ), seen as a polynomial in ζ , admits at least one real, positive root ζ = k 2 3 . The x-spectrum S x then coincides with the locus of the roots of P(ζ ; λ, q, p, r ) regarded as a polynomial in λ for all values of ζ > 0, namely for ζ considered as a real, positive parameter.
From the algebraic-geometric point of view, the locus of the roots of P(ζ ; λ, q, p, r ) in the λ-plane is an algebraic curve; this can be given implicitly as a system of two polynomial equations in two unknowns by setting λ = μ + i ρ and then separating the real and the imaginary parts of P(ζ ; λ, q, p, r ). In the following, we will apply Sturm's chains (e.g. see Demidovich and Maron 1981;Hook and McAree 1990) on Q(ζ ; q, p, r ), the discriminant of P(ζ ; λ, q, p, r ), and invoke Sturm's theorem to study its roots, in order to obtain a classification of the different loci S x .

B.1 Sturm Chains and Spectra Classification in the (r, p)-Plane
Let P(x) be a polynomial in x with real coefficients, and let deg(P) be its degree. A Sturm's chain (Demidovich and Maron 1981;Hook and McAree 1990) is a sequence of polynomials associated with a given polynomial P(x), and its derivative, where the notation Remainder(P 1 , P 2 ) stands for the remainder of the polynomial division between P 1 and P 2 . By Sturm's theorem, the sequence allows to find the number of distinct real roots of P(x), in a given interval, in terms of the number of changes of signs of the values of the Sturm sequence at the end points of the interval, which can even be taken to be ±∞. When applied to the whole real line, it gives the total number of real roots of P(x). Since deg(P (i+1) ) < deg(P (i) ) for 0 ≤ i < j, the algorithm terminates and it contains in general deg(P) + 1 polynomials. The final polynomial, P ( j) (x), is the greatest common divisor of P(x) and its derivative. For instance, suppose that we would like to find the number of distinct real roots of P(x) in the real interval [x i , x f ]; let ς(x i ) and ς(x f ) be the number of changes of sign of the Sturm chain at x i and x f , respectively; then, Sturm's theorem states that the number of real roots of P( For each value of ζ > 0, we have four values of λ, as P(ζ ; λ, q, p, r ) is a 4 th degree polynomial in λ. For all p, r real, there always exist four values of λ such that P W (w; λ, q, p, r ) has at least one double w-root. The nature of these four points can be classified in terms of the sign of λ w P W (w; λ, q, p, r ). Observe that P(0; λ, q, p, r ) coincides with w P W (w; λ, q, p, r ), the discriminant of P W (w; λ, q, p, r ) with respect to w, with reversed sign. Furthermore, the discriminant of P(0; λ, q, p, r ) with respect to λ, λ P(0; λ, q, p, r ), coincides with λ w P W (w; λ, q, p, r ). Thus, a first classification of S x can be achieved by imposing see Proposition 8. In order to obtain a complete classification of S x for all ζ , we observe that, by moving ζ in the interval (0, +∞), we move the four λ-roots of the polynomial P(ζ ; λ, q, p, r ) (regarded as a polynomial in λ). Two of such roots will collide if λ P(ζ ; λ, q, p, r ) = Q(ζ ; q, p, r ) ≡ Q(ζ ) vanishes. Since λ P(ζ ; λ, q, p, r ) = Q(ζ ) = 65536 Q 2 1 (ζ ) Q 2 (ζ ) , and, because Q 1 (ζ ) appears squared in the expression of Q(ζ ), the sign of Q(ζ ) depends solely on the sign of Q 2 (ζ ). Moreover, we observe that the real, positive roots of Q 1 (ζ ) are not associated with any change in the behaviour of λ, as Q(ζ ) does not change sign if ζ passes through one of those roots: indeed, if ζ passes through one of the roots of Q 1 (ζ ), then two λ-roots will collide, but after the collision they will remain on the branches that they occupied before the collision. It is easy to see, by isolating the leading term in ζ of Q(ζ ), that Q(ζ ) > 0 for large ζ . In particular, all the λ-roots of P(ζ ; λ, q, p, r ) are real if ζ is larger than the maximum of the real, positive roots of Q 2 (ζ ). After that point, there are no complex branches in the λ-plane, and this provides a simple limit for the largest value of k 3 = √ ζ for which the gain function is defined (see (79)). This also proves that the real axis or part of it is always included in the locus on the λ-plane. On the other hand, if Q(ζ ) < 0, then two of the λ roots are real and two of the λ-roots form a pair of complex conjugate values. By changing the value of ζ and keeping Q(ζ ) negative, two of the λ-roots will move on the real axis and the other two will move in the complex plane. When Q(ζ ) = 0, then two of the λ-roots collide in one point. Finally, if Q(ζ ) > 0 for ζ > 0 and smaller than the maximum of the real roots of Q 2 (ζ ), we can have either four real roots or two pairs of complex conjugate roots. Based on what we have seen so far, it is clear that the structure of the algebraic curves in the λ-plane (i.e. the x-spectrum S x ) is related to the number of changes of sign of Q(ζ ) for ζ ≥ 0, that is the number of changes of sign of Q 2 (ζ ) for ζ ≥ 0. In turn, as one varies p and r , the number of changes of sign of Q 2 (ζ ) for ζ ≥ 0 can be obtained by counting the number of positive roots of Q 2 (ζ ). In particular, as p and r vary, the number of positive roots of Q 2 (ζ ) changes in two possible ways: either two roots collide on the real ζ axis and then move onto the complex plane as a pair of complex conjugate roots (or vice versa), or a positive root passes through the origin and becomes negative (or vice versa). Therefore, in order to provide a complete classification of all the algebraic curves (i.e. the spectra) in the λ-plane, not only we need to understand the number of changes of sign of Q 2 (ζ ) for ζ ≥ 0 (which is the number of its positive zeros), but also the general structure of the zeros of Q 2 (ζ ) (namely, the nature of the its non-positive roots).
If ς(0) is the number of changes of sign for the Sturm chain {Q ( j) 2 } 4 j=0 at ζ = 0 and ς(∞) the number of changes of sign of the leading terms of the Sturm chain as ζ → +∞, then, by Sturm's theorem, the number of positive roots of Q 2 (ζ ) is simply given by ς(0) − ς(∞). Using this approach, we obtain that, for all r , p real, Q 2 (ζ ) has always at least one real, positive root. Moreover, in this way we obtain a set of ten algebraic curves in the (r, p)-plane, whose intersection defines a set of regions where either ς(0), or ς(∞), or both change value, namely where the structure of the zeros of Q 2 (ζ ) is expected to vary.

C Numerical Computation of the Spectra
In this section, we briefly illustrate how the spectra can be computed numerically. Let ζ j ≥ 0, j = 1, . . . ,n, with 1 ≤n ≤ 4, be the real, positive roots of the polynomial Q 2 (ζ ) (B.1c). We recall that Q 2 (ζ ) has always at least one real, positive root. Moreover, the non-real part of the spectrum S x is the locus on the λ-plane of the λ-roots of P(ζ ; λ, q, p, r ) for 0 ≤ ζ ≤ max jζ j . Thus, to numerically compute a spectrum for a given choice of r and p, with q = 1, it suffices to solve P(ζ ; λ, q, p, r ) = 0 for the same choice of the parameters, for a convenient set of values of ζ in the interval [0, max jζ j ]. These values will be referred to as ζ -nodes. As zeros of Q 2 (ζ ) are zeros of the discriminant of P(ζ ; λ, q, p, r ), with respect to λ, we expect that when ζ is close to one of theζ j 's the algebraic curve S x undergoes rapid changes. This implies that in order to capture all features of the spectrum, and optimize the root-extracting procedure, it is expedient to distribute the ζ -nodes according to a Chebyshev-Gauss-Lobatto distribution (Quarteroni et al. 2000) between each one of the intervals [0,ζ 1 ], [ζ 1 ,ζ 2 ],..., [ζn −1 ,ζn].
The simultaneous computation of all the roots of all polynomials has been performed via the standard technique of evaluating the eigenvalues of a companion matrix, as per implemented in MATLAB R2017a. Note also that all spectra have been verified using an objective-function technique over the imaginary part of the differences of the w j 's, computed by solving directly P W (ζ ; λ, q, p, r ) = 0.