Integrability and Linear Stability of Nonlinear Waves.

It is well known that the linear stability of solutions of 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 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 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.

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 1990, 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 is the perturbed solution of the same equation, then, at the first order of approximation, 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 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 around its continuous wave (CW) solution . 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 . 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 . Indeed, depending on boundary conditions, this technique yields a representation of the perturbation 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 . 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 (2000, 2002) for solutions which vanish sufficiently fast as . 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 , 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, (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 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 . 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 which therefore yields the long time evolution of 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 Grinevich and 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 . 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, , , which evolve according to the coupled system of NLS equationswhere is the self- and cross-interaction term. This is integrable only in three cases (Zakharov and Shulmann 1982), namely (after appropriate rescaling): (defocusing Manakov model), (focusing Manakov model) (Manakov 1973) and the mixed case . 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 , 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, 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 (), 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 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 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 matrix-valued complex functions. The existence of a fundamental (i.e. non singular) matrix solution of this overdetermined system is guaranteed by the condition that the two matrices X and T satisfy the differential equationWe recall here that, unless differently specified, a subscripted variable means partial differentiation with respect to that variable, and stands for the commutator . 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 and be polynomial in with degrees n and m, respectively. As a consequence, the matrix is as well polynomial in with degree , and therefore, compatibility Eq. (3) yields 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 , , which differs by a small change of the matrices X and T, with the implication that the pair of matrices and , at the first order in this small change, satisfies the linearized equationAgain, 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 ,and to search for its solutions and . Our main target is to find those solutions which are related to the fundamental matrix solution of Lax pair (2). To this purpose, we first note that the similarity transformationof a constant (i.e. x, t-independent) matrix yields the transformed matrix which, for any given arbitrary matrix M, satisfies the pair of linear ODEsEquation (7) are compatible with each other because of (3). Then, for future reference, we point out the following observations.

Proposition 1

If the pair solves linearized Eq. (5), then also the pairis 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 expressionsprovide a solution of linearized Eq. (9).

The validity of this statement follows from the fact that the matricesobviously 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 and , see (2), have degree one and, respectively, two, namelywhere , Q, , and 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 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 , , satisfy if , while is the identity matrix where is the (algebraic) multiplicity of the eigenvalue . Of course, . Note that this matrix induces the splitting of the set of matrices into two subspaces, namely that of block-diagonal matrices and that of block-off-diagonal matrices. Precisely, if is any matrix, then we adopt the following notationwhere is the block-diagonal part of and is its block-off-diagonal part. Consistently with this notation, the entries of an matrix are meant to be matrices themselves of dimension with the implication that the “matrix elements” 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 are generically rectangular. We also note that the product of two generic matrices and follows the rules: , while, for , the product 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 coefficientswith the following comments and definitions. The matrices , , are constant and block-diagonal, . In the following, we set , 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 propertieswhich show that also the matrix is block-off-diagonal. Also the maps , , act only on block-off-diagonal matrices according to the definitionsFinally, the matrices , , 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 , , which ultimately implies restrictions on the constant matrices and , as specified by the following

Proposition 3

The matrices and identically vanish if, and only if, the blocks of and 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 readsThis equation can be linearized around a given solution Q(x, t) by substituting Q with and by neglecting all nonlinear terms in the variable . In this way, we obtain the following linear PDEWe are now in the position to formulate our next proposition which is the main result of this section.

Proposition 4

The matrixdefined by (10), together with (12), satisfies the same linear PDE (22) satisfied by if and only if the block-diagonal matrices and satisfy the same conditions (20) which guarantee the local character of evolution Eq. (21).

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 . To this purpose, one starts by rewriting the first of ODEs (7) for the block-diagonal and block-off-diagonal components of , namelyNext, one considers the time derivative by using expression (23), along with the right-hand-side commutator of the second relation in (7), namely , and one replaces the matrix with its expression as in (12) and (16). Then, all terms of the form , 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 do indeed cancel out, provided conditions (20) are satisfied. The remaining terms of the expression of can be rearranged as to coincide with those which appear in the right-hand-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 , 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 . 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 , and 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 , and . By setting instead and , one obtains the three-wave resonant interaction system for (Kaup 1976b) and many-wave interaction type equations for Ablowitz and Segur (1981) (e.g. see Calogero and Degasperis 2004; Degasperis and Lombardo 2007, 2009).

x-Spectrum of the Solution Q(x, t)

The result stated by Proposition 4 in the previous section implies that any sum and/or integral of over the spectral variable is a solution of linear Eq. (22). Here and in the following, we assume that such linear combination of matrices , which formally looks as the Fourier-like integral representationis such to yield a solution of (22) which is both bounded and localized in the x variable at any fixed time t (periodic perturbations are not considered here). The boundedness of implies that the matrix 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 defines a special subset of the complex -plane, over which the integration runs, which will be referred to as the x-spectrum 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 (like, for instance, when its entries are in ), then coincides with the spectrum of the differential operator defined by the ODE (i.e. the x-part of Lax pair (2)). However, if instead Q(x, t) goes to a non-vanishing and finite value as , this being the case for continuous waves, the spectrum associated with the solution Q(x, t) does not coincide generically with the spectrum of the differential operator . In the present matrix formalism, this happens for , and it is due to the fact that the spectral analysis applies to the ODE , see (7), rather than to the Lax equation itself. We illustrate this feature for in the next section, where we consider the stability of CW solutions. The spectrum consists of a piecewise continuous curve and possibly of a finite number of isolated points. Note that at any point of the spectrum the matrices , 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 and we rather consider solutions of (22) which have integral representation (26) and vanish sufficiently fast as , this being a requirement on the matrix 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 remains small as time grows, say . Thus, the basic ingredient of our stability analysis is the time dependence of the matrices for any . As we show below, getting this information requires knowing the spectrum whose computation therefore is our main goal.

Wave Coupling and Spectra: An Example

In order to illustrate the results of the previous section and to capture, at the same time, a class of nonlinear wave phenomena of physical relevance (including the system of two coupled NLS equations), the rest of the paper will focus on the matrix case , , withandHere and below, the asterisk denotes complex conjugation and the four field variables are considered as independent functions of x and t and are conveniently arranged as two two-dimensional vectors, that isWhile 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 , and T-matrix (16), with , , specializes to the expressionThen, matrix PDE (21) becomeswhich is equivalent to the two vector PDEsHere, 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. (31), linearized Eq. (22) for a small change readsMoreover, Proposition 4 guarantees that the matrix satisfies this same linear PDE, namelyand, for , 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 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 , 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 () 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 self-interaction 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 , , , and, in a crucial manner (see below), on the real parameter q which measures the wave number mismatch between the two wave components and (or and ).

The main focus of this section is understanding how the spectrum changes by varying the parameters , , , 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) readswhere the matrix has expression (27), while the matrix isand we conveniently introduce the real parameters which will be handy in the following. Next we observe that a fundamental matrix solution of the Lax equations has the expressionwhere the x, t-independent matrices W and Z are found to be with the property that they commute, , consistently with compatibility condition (3). In order to proceed with plain arguments, we consider here the eigenvalues and , , of and, respectively, of as simple, as indeed they are for generic values of . In this case, both and are diagonalized by the same matrix , namelyNext we construct the matrix via its definition, see (23) and (6),which, because of explicit expression (40), readsAs for the matrix , it lies in a nine-dimensional linear space whose standard basis is given by the matrices , whose entries arewhere is the Kronecker symbol ( if and otherwise). However, the alternative basis , which is obtained via the similarity transformationwhere diagonalizes W and Z (see (43)), is more convenient to our purpose. Indeed, expanding the generic matrix in this basis asthe scalar functions being its components, and inserting this decomposition into expression (45), leads to the following representation of Fwhere we have introduced the x, t-independent matricesThe advantage of expression (49) is to explicitly show the dependence of the matrix F on the six exponentials .

Proposition 4 stated in the previous section guarantees that, for any choice of the functions , expression (49) be a solution of linearized Eq. (33), see (34). It is plain (see (26)) that the requirement that such solution be localized in the variable x implies the necessary condition that the functions be vanishing for , , . The further condition that the solution be bounded in x at any fixed time t results in integrating expression (49) with respect to the variable over the spectral curve of the complex -plane (see (26)):Here, according to Sect. 2.1, can be geometrically defined as follows:

Definition 1

The x-spectrum , 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 , , or explicitlyis real.

Observe that the ’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 of (33). Indeed, it is sufficient to compute the eigenfrequenciesas suggested by the exponentials which appear in (49). Their expression follows from matrix relation (42)and readThis expression looks even simpler by using the relation implied by the trace of the matrix (see (41)) and finally readsThe consequence of this expression (56), which is relevant to our stability analysis, is given by the following

Proposition 5

Continuous wave solution (35) is stable against perturbations whose representation (26) is given by an integral which runs only over those values of which are strictly real.

Proof

If the spectral variable is real, then all coefficients of the characteristic polynomial of the matrix (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 of are either all real, or one real and two complex conjugate. In the first case, three wave numbers (52) are all real, and thus, the corresponding is in the spectrum , . In the second case, the three 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 , . In the following, we refer to one such forbidden real interval as a gap, see Sect. 3.1. Consequently, in the first case, since , , , and therefore, , , , are all real, then also , , , 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, , with non-vanishing imaginary part, . Let be the (generically complex) roots of . With this notation, we have that one of the wave numbers, say , will be real only if . Then, from (56), we have that will also be real only if . Writing the polynomial as , 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 , , , , 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 , , , , and , only if (unless the non-physical and non-generic condition be met). This contradicts the original assumption . Therefore, for a representation (26) of a perturbation to be bounded (and thus for the corresponding CW solution to be stable), the integral in (26) must run only on those values of which are strictly real.

This Proposition 5 implies that a real part of the spectrum 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 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 (see (35)). Indeed, if the stability property of the continuous wave is essentially that of the NLS equation. This conclusion results from the explicit formulae , , , and thereforewhich show that if (see (39a)), the spectrum is the real axis with the exclusion of the gap , while, if , the spectrum is the entire real axis with the addition of the imaginary branch , . Therefore, the stability for follows from the reality of the frequencies , , over the whole spectrum , see (58). The modulational instability occurs only if since, for in the imaginary branch of the spectrum , the frequencies , and have a non- vanishing imaginary part.

In the following, while computing the spectrum , we consider therefore only the case and, with no loss of generality, strictly positive, . 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, , while keeping in mind that, by doing so, , rescale as q, while p, r, rescale as , 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 -plane with .

Indeed, the change of sign , takes the characteristic polynomial (57) into 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 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 of 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 of the polynomial (57) with respect to w, obtaining thus a polynomial in , with parameters q, p, r,which is positive whenever the three roots 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, , for those real values of which belong to a gap of , while its zeros are the end points of gaps. Here and hereafter, the notation stands for the discriminant of the polynomial P with respect to its variable y.

In the exceptional case , computing the roots of the discriminant as a polynomial in the variable reduces to the factorization of a quadratic polynomial, allowing the formulation of the following

Proposition 7

For , the gaps depend on the parameter p as follows:

The border values areThe threshold values and correspond to the opening and the coalescence of two gaps, respectively.

In the generic case , we find that the number of gaps of the real part of the spectrum is either zero, or one, or two. Gaps appear at double zeros of the discriminant , thus at zeros of its own discriminant , namely when The three polynomial factors appearing in (61) bound the regions of the -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 . In order to find the number of gaps from expression (59) of the discriminant , it is convenient to plot the dependence of the parameter on the spectral variable at the zeros of this discriminant for fixed values of the parameters r and q. This function is therefore implicitly defined by the equation . For a fixed value of q, in the -plane the 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 as a function of r and q is Furthermore, seen as an algebraic curve in the -plane satisfies the following implicit relationThe transition between the two regimes (two minima, or one minimum and one maximum) takes place at the special threshold value where (see “Appendix A”). If , the cusp is also a local maximum. Equipped with these findings, we formulate

Proposition 8

For real and , has the following gap structure.

If

The threshold values , and correspond to the opening of one gap, the opening of two gaps and the coalescence of two gaps, respectively.

If

The threshold values , and 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 areand the curve 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 . These curves will play a role also in the next subsection where we give a complete classification of branches and loops of the spectrum.

Fig. 1

The three threshold curves (solid black) and (solid red) for , as parametrically defined by (64) and (63) or explicitly by (62), are plotted. They are boundaries of regions of the (r, p)-plane, , where the number of gaps, either 0 or 1, or 2, is shown. It is also shown that for and that for (see Proposition 8) (Color figure online)

Branches and Loops

In the previous subsection, we have described the gaps of the spectrum , 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 . This subset is made of smooth and generically finite, open (branches) or closed (loops), continuous curves (with the possible exception of the case , see Fig. 3h and Degasperis et al. 2018). By Definition 1, if at least one of the corresponding three wave numbers , is real, see (52). Going back to the characteristic polynomial (57), we look now at its roots , and at the wave numbers , , see (52). If is not real, the roots cannot be all real themselves since the coefficients of are not real. As a consequence, the requirement that one of the wave numbers (52) be real implies that at least two roots of , for instance and , have the same imaginary part. In fact, branches and loops of may exist only in this case: indeed, if the three roots have all the same imaginary part, it can be shown that no smooth branch or loop exists. Hence, here and below, we name the only real wave number, while and are complex (note that ). The complex part of the spectrum is therefore defined as the set of the -plane such that . In order to compute this component of the spectrum, we introduce the novel polynomial ,defined by the requirement that its roots are the squares of the wave numbers, . It is plain that its coefficients are completely symmetric functions of the roots of the characteristic polynomial . This is evidently so because , , are symmetric functions of , , , and therefore of , , , with the implication that the coefficients , , of are polynomial functions of the coefficients , , of characteristic polynomial (57), These relations, which readcombined with the expressions of the coefficients , see (66), yield the dependence of the coefficients of on the parameters r, p, q and on the spectral variable , so that . This dependence turns out to bewhere is the discriminant of , see its expression (59), with the implication that 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 ) is real. Therefore, our task of computing the complex part of the spectrum , for a given value of the parameters r, p, reduces to finding the curve in the -plane along which the root of remains real and positive, . 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 as the locus in the -plane of the -roots of seen as a real polynomial in , for spanning over the semi-line . In other words, for a fixed value of , we compute the four -roots of the real polynomial , which are the values of (irrespective of being complex or purely real) such that the wave number is real. In this way, the two cases of and can be treated at once (allowing to retrieve and confirm all the results about gaps in the spectra presented in the previous subsection).

Fig. 3

Examples of spectra for different values of r and p, with . a for , , i.e. , , , as an example of a 0G 2B 0L spectrum. b for , , i.e. , , , as an example of a 0G 2B 1L spectrum. c for , , i.e. , , , as an example of a 1G 1B 0L spectrum. d for , , i.e. , , , , as an example of a 1G 1B 1L spectrum. e for , , i.e. , , , as an example of a 2G 0B 1L spectrum. f for , , i.e. , , 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 for , , i.e. , , 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 , two symmetrical gaps open on the imaginary axis for , as explained in Degasperis et al. (2018). h for , , i.e. , , 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)

From the algebraic-geometric point of view, the locus of the roots of in the -plane is an algebraic curve; this can be given implicitly as a system of two polynomial equations in two unknowns by setting and then separating the real and the imaginary parts of .

In order to explain how the spectrum changes over the parameter space, we apply Sturm’s chains (e.g. see Demidovich and Maron 1981, Hook and McAree 1990) to the polynomial , that is, to the discriminant of the polynomial with respect to . Indeed, it turns out that the nature and the number of the components of the x-spectrum , seen as a curve in the -plane, are classified by the number of sign changes of , for . 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, . 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 (64) and (62) introduced in Sect. 3.1, one more threshold curve, , is found, which, for a given non-vanishing value of q, is implicitly defined asor, explicitly, asIt 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 , , (with ) are plotted to show the partition of the parameter half-plane according to the gap, branch and loop components of the spectrum, for both and . This overall structure of the spectrum in the parameter space can be summarized by the following

Fig. 2

The structure of the spectrum in the (r, p) plane, (see Proposition 9). In red, the curve ; in blue, the curve ; in grey, the curves . a . b (Color figure online)

Proposition 9

For , the spectrum has the following structure in the parameter plane .

If

If

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 (entailing , 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 .

Examples of different spectra have been computed numerically by calculating the zeros of 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 -space as reported in Proposition 9.

Finally, the limit case of the parameter half-plane deserves a special mention. In fact, four different limit spectra are found for , namely in the four intervals , , and . 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 allows for analytic computations, and explicit formulae, since, as already noted in Sect. 3.1 for the discriminant , the coefficients of the relevant polynomial (65) reduce to polynomials of degree 2 in the variable .

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 equationhas 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 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 in (32), whereThe 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 , has been fixed, linearized equation (33) around this solution areThe two coupling constants , , if non-vanishing, are just signs, , 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 , as (F/F) if and as (D/F) in the mixed case .

Hereafter, our focus is on the stability of the CW solution (see (35), (36) with )where the parameter q is the relative wave number, which may be taken to be non-negative , and , are the two amplitudes, whose values, with no loss of generality, may be real and non-negative, . 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 At this point, we show in Fig. 4 the -plane as divided into octants, according to different values of amplitudes and coupling constants.

Fig. 4

The -plane divided according to amplitudes and coupling constants , (Color figure online)

We note that, as pointed out by Proposition 6, it is sufficient to confine our discussion to the half-plane , 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 relationseach octant being characterized as shown in Fig. 4. In particular, the octant will be referred to as , whereas the octant as , to indicate the wave with larger amplitude. Moreover, we observe again that setting , 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, , and is unstable in the opposite case . However, this solution remains stable even if the wave feels a self-focusing effect, i.e. , provided the other wave in the defocusing medium, , has larger amplitude . Similarly, when , the CW solution is unstable if and , and if the larger amplitude wave is the one which propagates in the focusing medium, . These remarks follow from explicit formulae (58) and (75a), and they are clearly evident in Fig. 4 where, for , the upper two octants correspond to and the lower two to .

The limit 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 , then for no value of the amplitudes , with , and for no value of the coupling constants , , the spectrum 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 for this statement to be true coincides with the condition that the point (r, p) neither belongs to the threshold curve , nor to the threshold curve , see (75). In particular, in the case , which is equivalent to setting and , the roots of , (57), can be given in closed form,so thatwith the implication that has to be real to be in . 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 , namely or, equivalently, , 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 , then CW solution (74) is unstable.

The relevant point here is the dependence on of the frequencies over the spectrum . As we have shown in the previous section, if , the spectrum consists of two components: one, , is the real axis , with possibly one or two gaps (see Proposition 8), and the second one, , consists of branches and/or loops where runs off the real axis (see Proposition 9). Therefore, spectral representation (51) of the perturbationtakes the formwhich is meant to separate the contribution to 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 do not play a role here and are not specified. On the contrary, the reality property of the frequencies over the spectrum is obviously essential to stability. Since the reality of the frequencies , , for real has been proved in Proposition 5, it remains to show that indeed is not real if belongs to branches or loops. This follows from explicit expression (56), namelywhere is real but and cannot be real and cannot have same imaginary part as proved in Proposition 5, (see Sect. 3.2). Here, the imaginary part of defines the gain function over the spectrum. This (possibly multivalued) function of plays an important role in the initial stage of the unstable dynamics. Precisely, its dependence on the wave number 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 , corresponding to integrable cases. In this respect, it is convenient to translate the classification of the spectra in the -plane, in particular, and with no loss of generality, in the quadrant , . The four threshold curves , , which have been found in the (r, p)-plane, see Fig. 2, are reproduced below in the -plane, see Fig. 5, according to coordinate transformations (75) and (76). The outcome of this analysis is summarized by the following

Fig. 5

The -plane (see Proposition 10). In Fig. 5a and 5b, grey portions correspond to . a (D/D) . b (F/F) . c (D/F) (Color figure online)

Proposition 10

For each of the three integrable cases of two coupled NLS equations (1), the spectrum is of the following types.

(see Fig. 5a)

The lower octant , which corresponds to , gets divided into two parts by the threshold curve which intersects the line at the point . In the lower, finite, part of this octant, the spectrum is of type , while in the other part it is of type .

(see Fig. 5b)

The upper octant , which corresponds to , gets divided into two parts by the threshold curve which intersects the line at the point . In the lower, finite, part of this octant, the spectrum is of type , while in the other part it is of type .

(see Fig. 5c)

The whole quadrant corresponds to . It is divided into three infinite portions by the two threshold curves , in the upper octant, and in the lower one. The spectrum is of type in the upper part, of type in the middle part and of type in the lower part.

These statements are straight consequences of Proposition 9 via transformation (76).

We conclude by mentioning the limit case , 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 , namely with the case in which the two wave amplitudes are strictly equal. On this particular line of the -plane, i.e. 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 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 , associated with the unperturbed solution of the wave equation. It is worth stressing that the spectrum for does not coincide with the spectrum of the Lax equation in the complex -plane. This result is explicitly shown for in the instance of continuous wave solutions of two CNLS equations. In this case, we define as functions of on the spectrum the set of wave numbers and frequencies 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 , . 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).

-∞<p<0	no gap
0≤p<q^2	two gaps, -g_+≤λ≤-g_- and g_-≤λ≤g_+
q^2≤p<+∞	one gap, -g_+≤λ≤g_+

-∞<p<-r	no gaps
-r≤p<r	one gap
r≤p<p_S(r)	two gaps
p_S(r)≤p<+∞	one gap

-∞<p<-r	no gaps
-r≤p<p_S(r)	one gap
p_S(r)≤p<r	two gaps
r≤p<+∞	one gap

-∞<p<p_C(r)	0G2B0L
p_C(r)<p<-r	0G2B1L
-r<p<r	1G1B0L
r<p<p_S(r)	2G0B1L
p_S(r)<p<+∞	1G1B0L

-∞<p<-r	0G2B0L
-r<p<p_C(r)	1G1B1L
p_C(r)<p<p_S(r)	1G1B0L
p_S(r)<p<r	2G0B1L
r<p<+∞	1G1B0L

-∞<p<p_C(r)	2 positive, 0 negative and 2 complex conjugate roots
p_C(r)<p<-r	4 positive, 0 negative and 0 complex conjugate roots
-r<p<r	1 positive, 1 negative and 2 complex conjugate roots
r<p<p_S(r)	2 positive, 2 negative and 0 complex conjugate roots
p_S(r)<p<-p_C(r)	1 positive, 3 negative and 0 complex conjugate roots
-p_C(r)<p<+∞	1 positive, 1 negative and 2 complex conjugate roots

-∞<p<-r	2 positive, 0 negative and 2 complex conjugate roots
-r<p<p_C(r)	3 positive, 1 negative and 0 complex conjugate roots
p_C(r)<p<-p_C(r)	1 positive, 1 negative and 2 complex conjugate roots
-p_C(r)<p<p_S(r)	1 positive, 3 negative and 0 complex conjugate roots
p_S(r)<p<r	2 positive, 2 negative and 0 complex conjugate roots
r<p<+∞	1 positive, 1 negative and 2 complex conjugate roots