Versatile rogue waves in scalar, vector, and multidimensional nonlinear systems

This review is dedicated to recent progress in the active field of rogue waves, with an emphasis on the analytical prediction of versatile rogue wave structures in scalar, vector, and multidimensional integrable nonlinear systems. We first give a brief outline of the historical background of the rogue wave research, including referring to relevant up-to-date experimental results. Then we present an in-depth discussion of the scalar rogue waves within two different integrable frameworks—the infinite nonlinear Schrödinger (NLS) hierarchy and the general cubic-quintic NLS equation, considering both the self-focusing and self-defocusing Kerr nonlinearities. We highlight the concept of chirped Peregrine solitons, the baseband modulation instability as an origin of rogue waves, and the relation between integrable turbulence and rogue waves, each with illuminating examples confirmed by numerical simulations. Later, we recur to the vector rogue waves in diverse coupled multicomponent systems such as the long-wave short-wave equations, the three-wave resonant interaction equations, and the vector NLS equations (alias Manakov system). In addition to their intriguing bright–dark dynamics, a series of other peculiar structures, such as coexisting rogue waves, watch-hand-like rogue waves, complementary rogue waves, and vector dark three sisters, are reviewed. Finally, for practical considerations, we also remark on higher-dimensional rogue waves occurring in three closely-related (2  +  1)D nonlinear systems, namely, the Davey–Stewartson equation, the composite (2  +  1)D NLS equation, and the Kadomtsev–Petviashvili I equation. As an interesting contrast to the peculiar X-shaped light bullets, a concept of rogue wave bullets intended for high-dimensional systems is particularly put forward by combining contexts in nonlinear optics.

Manakov system). In addition to their intriguing bright-dark dynamics, a series of other peculiar structures, such as coexisting rogue waves, watchhand-like rogue waves, complementary rogue waves, and vector dark three sisters, are reviewed. Finally, for practical considerations, we also remark on higher-dimensional rogue waves occurring in three closely-related (2 + 1)D nonlinear systems, namely, the Davey-Stewartson equation, the composite (2 + 1)D NLS equation, and the Kadomtsev-Petviashvili I equation. As an interesting contrast to the peculiar X-shaped light bullets, a concept of rogue wave bullets intended for high-dimensional systems is particularly put forward by combining contexts in nonlinear optics.
Keywords: rogue wave, modulation instability, integrable turbulence (Some figures may appear in colour only in the online journal)

Motivation and historical background
The field of rogue waves is currently one of the most active multidisciplinary areas of research encompassing oceanography, hydrodynamics, optics and photonics, plasma physics, and Bose-Einstein condensation [1,2]. Scientists from a wide range of disciplines acknowledge that this general rogue wave concept can describe in a relevant way extreme wave events of different origins, from the killer waves in the ocean [3][4][5][6] to the extreme pulses in optics [7,8], and even the financial crises in economy [9].
Considering the ubiquity of extreme wave dynamics, and the multiplicity of their contextual denominations, it is important to begin our review article by providing an accessible rogue wave definition that most scientists would be inclined to approve. Indeed, this issue of shaping up a unifying rogue wave concept has been debated in past years [1]. It is therefore of concern to recall the observational origin of this concept arising from oceanography. Whereas it had been nourished long ago by sailor folklore of monster or killer waves, it was scientifically confirmed by live accounts of maritime disasters involving tankers and liners as well as anomalous water elevation recordings on offshore platforms, such as the henceforth notorious New Year's wave at the Draupner platform in the North Sea, recorded on January 1st, 1995 [3].
From scientific accounts in oceanography, three main features have arisen to qualify a rogue wave phenomenon [7,10]. The dynamics should manifest rare transient events of giant amplitude. This is first quantified by comparing the wave amplitude to the significant wave height (SWH), defined as the average amplitude of the highest third of the wave events. To qualify, a wave transient should have a maximum amplitude more than twice (or larger than 2.0-2.5, according to different authors) the SWH. Secondly, at the observational level, the wave transient should appear and disappear unexpectedly. Even though this second criterion appears to be merely qualitative, it remains essential to rule out high-amplitude solitary waves from the picture. For instance, the second criterion excludes the tsunamis in hydrodynamics, since once triggered by a primary catastrophic event, usually of geological origin, their propagation can be followed over large distances. For a similar reason, in photonics, the second criterion rules out stable mode-locked laser pulses, whatever gigantic their peak power may be.
It is always conceivable that transient events of very large amplitude might be resulting from the linear superposition of random waves [11]. However, in that situation, the likelihood of waves transients exceeding twice the SWH would be extremely scarce. Such likelihood of wave events is assessed by recording the probability distribution function (PDF) of the wave amplitudes. Since linear superposition models generally entail Gaussian distributions, the recurrent observation of rogue wave phenomena, beginning from the field of oceanography and flourishing in photonics, has given nonlinearity an essential role in most rogue wave contexts. Therefore, comes the third criterion, stating that rogue wave events should arise more frequently than they would if they were obeying a classical Gaussian distribution. This is an indispensable technical aspect in rogue wave experimentation, which is the scientific translation of the wicked maritime experience of witnessing a catastrophic event that should not have reasonably occurred during the lapse of time spent on the sea, would the forecast statistics be based on a classical probability distribution. Hence, the typical PDF of rogue wave dynamics will display an L-shaped distribution of amplitudes [12][13][14][15][16].
Rogue wave observation and experimentation in various propagation media has naturally been a strong incentive for fundamental rogue wave modelling, mathematically and numerically. As an integrable scalar equation of nonlinear wave propagation in one dimension, with versatile applications from deep water hydrodynamics to nonlinear fibre optics, the nonlinear Schrödinger (NLS) equation has played a central role in the systematic investigation of rogue wave phenomena. In 1983, Howell Peregrine proposed a first-order rational solution of the NLS equation, endowed with localization in both space and time (1 + 1) dimensions [17]. This solution, although usually termed as 'Peregrine soliton', is not a soliton in the ordinary sense, since its profile never stays stationary. Instead, this doubly-localized pulse gradually develops from a background excitation, reaches a climax worth three times the background amplitude, then vanishes back into the background in an otherwise symmetrical fashion. Its high peak amplitude as well as its asymptotic connection with a constant background, which supports the aesthetic idea that this extreme wave comes from nowhere and disappears without leaving a trace [18], has made the Peregrine soliton a central prototype of rogue wave manifestations, even though the notion of statistical distributions, which is meant to confront the experimental appearance of rogue waves among chaotic wave fields, namely, our third criterion, has been mostly absent to date from mathematical investigations. Nevertheless, the Peregrine soliton has until now continued to play a pivotal role for modelling rogue waves in scalar NLS systems, as well as a gauge to assess other rogue wave solutions proposed for systems other than the NLS equation. Analytically, it can be considered as the limiting case of breather solutions on a finite background, when their recurrence period tends to infinity [19]. Therefore, it is anticipated that Peregrine-type solitons will bear a special relationship with modulation instability (MI), which we shall specify and illustrate at several places along our review.
In nonlinear science, the diversity of the nonlinear propagation media that can be considered in realistic settings entails a wide range of modelling equation candidates [20]. Above all, in some media, the dominating nonlinearity will be quadratic, instead of cubic as in the NLS equation. In addition, considering the bandwidth and intensity of the medium excitation, the inclusion of higher-order dispersive and/or nonlinear terms can become a necessity to reach an acceptable modelling accuracy. Finally, the vector nature of interacting wave fields can also manifest, as well as the higher dimensionality of the propagation, when we go beyond one-dimensional guided propagation. Thus, the possible existence of rogue wave phenomena in various experimental fields of nonlinear science justifies their mathematical and numerical investigations within a wide range of propagation equation models. Note that in the following, we have privileged the mathematical exploration of integrable equations, in order to explicitly obtain exact rational rogue wave solutions. Therefore, dissipative nonlinear systems [21] are not discussed here, albeit their practical importance in several areas of rogue wave experimentation [14], from current-driven hydrodynamics to ultrafast lasers [3,8,22]. Also, the constraint of integrability often dictates a specific relationship between the model parameters, weighing in a fixed fashion several physical effects that otherwise can feature a large variability among realistic propagation media. Nevertheless, the mathematical exploration of exact rational solutions in various integrable models can be considered as the safest ground to discover new possibilities of rogue wave existence, endowed with novel features. Sometimes, the latter can even be extrapolated numerically to closely related models whose parameter sets break the integrability, but still possess qualitatively similar solutions [23].

Outline of the review
The structure of the review is organized in the following. In section 2, a short peregrination around soliton solutions on a background of the scalar NLS equation will bring us to the Peregrine soliton, the premier rogue wave mathematical prototype (section 2.1.1). We then present a series of exact rogue wave solutions recently obtained from scalar nonlinear equations that can be considered as extensions of the NLS equation. This extension begins with the Hirota equation (section 2.1.2), and the infinite NLS hierarchy (section 2.1.3), which generalizes the construction of higher-order integrable equations. These extended equations do also accept the rational Peregrine soliton as a valid solution, exclusively in the self-focusing case. The Sasa-Satsuma equation is then presented as an interesting situation where the structure of the rogue wave solution can significantly differ from the Peregrine soliton, with the possibility to involve two peaks instead of a single one, and four holes (section 2.1.4).
Whereas the ubiquity of rational rogue wave solutions among a wide range of nonlinear propagation equations may look surprising, it provides a hint for a common generating mechanism. As mentioned above, the Peregrine soliton can be related to a limiting case of MI, when the modulation period tends to infinity. Therefore, the existence of such rational solution would be precisely contingent on the existence of a net positive MI gain for any arbitrarily small modulation frequency, namely, the existence of baseband MI. This fundamental relationship is developed in section 2.1.5, and illustrated for other propagation equations subsequently investigated.
Analytical rational solutions provide us with prototypical rogue wave profiles that are essential in experimental confrontation. Nevertheless, to assess whether rogue waves can be excited in realistic conditions, it is also crucial to undertake systematic numerical simulations that incorporate noise into the continuous-wave (cw) background. Noise and nonlinearity will produce turbulent wave fields, over which rogue wave events can be statistically analysed. This approach is presented in section 2.1.6, in the framework of the integrable NLS equation, revealing, besides Peregrine solitons, the impact of breather and soliton collisions on the generation of localized extreme events.
Adding higher-order terms with more parameters into the propagation equation is also a way to produce new possibilities of balancing several physical effects, which is known to impact the domains of existence of soliton solutions. This is also true for rational solutions: we show in section 2.2 that a cubic-quintic extended NLS equation allows for rational rogue waves solutions in both focusing and defocusing nonlinearity cases, a feature absent from the NLS hierarchy. These Peregrine-like rational solutions feature a strongly localized frequency chirping effect, hence their denomination of chirped Peregrine solitons. Rogue waves formed with either focusing or defocusing nonlinearity are also illustrated within the Fokas-Lenells equation (section 2.2.3), revealing the subtle role of the self-steepening term that plays in rogue wave formation.
Coupled nonlinear systems also constitute an abundant class of experimentally relevant dynamical systems. They consist of interacting wave components of different modes, frequencies or polarizations, commonly termed vector systems. Similar to the inclusion of additional physical effects and free parameters, it is known that the model extension from a scalar to a vector one usually expands the domains of existence for MI and solitons. We also expect the same with respect to rogue waves. Moreover, the vector character will provide additional possibilities for rogue wave profiles, such as coupled bright-dark and dark-dark rogue waves.
These possibilities are explored in section 3, starting in section 3.1 with the long-wave short-wave (LWSW) resonant system, which is a general parametric process occurring when the group velocity of the short wave matches the phase velocity of the long wave. The (1 + 1)-component LWSW system demonstrates dark-bright rogue waves, whereas the study of the (2 + 1)-component LWSW equation reveals the unusual coexistence of two different rogue waves solutions on a common cw background. The three-wave resonant interaction system (TWRI) is another popular vector system that is relevant in many areas, such as in the nonlinear optics of quadratic propagation media. Its investigation in section 3.2 unveils bright-dark-bright and dark-dark-bright rogue wave triplets, along with some peculiar rogue wave patterns such as watch-hand-like (WHL) super rogue waves and complementary rogue waves. In section 3.3, the Manakov system, which models for instance the propagation in optical fibres with randomly varying birefringence, is shown to exhibit vector dark rogue wave solutions, as well as higher-order solutions that take the form of vector dark triplets, which we dub dark three sisters, the dark counterpart of the three-sister rogue waves reported to occur in the Great Lakes in North America. All the models investigated in section 3 provide the occasion to verify the direct relationship between the existence of rogue wave solutions and that of baseband MI. Particularly, the vector dark rogue waves predicted for the Manakov system have just recently observed in telecommunication fibres.
As soon as the fields are no longer confined into unidimensional waveguides, the extension to higher-dimensional models represents a physically natural albeit mathematically challenging task, as accounted for in section 4. Conserving integrability as a guiding rule, several (2 + 1)D models are presented: Finally, the last section 5 concludes this review and proposes an outlook for future possible directions of theoretical research in this fast growing field.  [24] plays a central role in nonlinear science, and is well known for its relevance to capturing any weakly nonlinear evolution of narrow-band processes [19]. For instance, it can describe well the small-amplitude gravity waves in deep water [25], the Langmuir waves in hot plasmas [26], the beam/pulse propagation in optical waveguides or fibres [20], and even the Bose-Einstein condensates confined to highly anisotropic cigar-shaped traps in the meanfield regime [27,28]. In dimensionless variables, this equation can be cast into the following standard form

Scalar rogue waves
where u(ξ, τ ) is the complex wave envelope traveling along ξ, with τ being the transverse variable. Along this review, the subscripts stand for partial derivatives, unless otherwise stated. In optical contexts, while the second term on the lefthand side of equation (2.1) denotes the group-velocity dispersion (GVD) or diffraction responsible for the pulse broadening or beam divergence, the third term accounts for the self-focusing (σ = 1) or self-defocusing (σ = −1) Kerr nonlinearity that acts to possibly counterbalance the former effect. Here for brevity, we fix the GVD in (2.1) to be anomalous while letting the nonlinearity sign be changeable, which can be equivalently transformed to the case with a fixed nonlinearity sign but a variable dispersion [20].
Since the pioneering work of Zakharov and Shabat in 1971 [29], there has been an intense and continued investigation into the analytic solutions of the NLS equation, both for their intrinsic scientific interest and for their potential to provide new insights into many important applications. Thanks to integrability, the NLS equation is now found to admit several exact elementary analytical solutions (solitons, cnoidal waves, breathers), but one in particular describing a localized soliton on a finite background is of much concern as the latter can help understand the physics of extreme waves. The first solution of this type is the Kuznetsov-Ma (KM) soliton, which was derived by Kuznetsov in 1977 [30] and separately by Ma in 1979 [31]. Considering the general plane wave background given by the KM soliton can be expressed as where m = σa 2 sinh(2φ), n = 2 √ σa sinh(φ) and φ ∈ R. The second breather-type solution on a finite background was originally presented by Akhmediev in 1986 [32], starting with the MI, hence the term Akhmediev breather (AB). On the same background (2.2), we can express AB as where now m and n are defined by m = σa 2 sin(2φ) and n = 2 √ σa sin(φ), with φ ∈ R. We should point out that, compared to their original forms shown in [30][31][32][33], both solutions (2.3) and (2.4) have been generalized by including two extra free parameters ω and σ. Particularly, as σ = −1, they are still, although singular, solutions of equation (2.1) and are easy to verify by recalling that cos(ix) = cosh(x). Despite the similarity between solutions (2.3) and (2.4) (basically with trigonometric cosines and hyperbolic cosines interchanged), the former indeed exhibits distinctly different dynamics than the latter. As seen in figures 1(a) and (b), the KM soliton periodically breathes along the propagation direction ωξ − τ = 0, whereas the AB always oscillates in the transversal direction, but with its localized peaks being oriented along the line ωξ − τ = 0. The period of both breathers is determined by the free real parameter φ; the larger the value φ, the more intensely the breathers oscillate. Recently, both ABs and KM solitons with ω = 0 were successively realized in optical fibres [34,35] and then observed at the same time in a water wave experiment [36]. However, both of them are only observable for σ = 1, as predicted by (2.3) and (2.4).
Interestingly, as the frequency parameter φ approaches zero, either the solution (2.3) or (2.4) can boil down to the much simpler rational form which is also an exact solution of the NLS equation (2.1). This rational solution was first proposed in 1983 by Howell Peregrine in the domain of hydrodynamics [17], and hence the name Peregrine soliton. It is now evident that the Peregrine soliton corresponds to the limiting case of either KM solitons or ABs when the period tends to infinity [33][34][35][36], and only exists in the self-focusing regime (as σ < 0, the solution (2.5) will be singular and thus physically irrelevant). To date, Peregrine solitons have been observed in a water wave tank [33], in optical fibres [34], and in plasmas [37]. One may notice that the solutions (2.3)-(2.5) presented above indeed hold true for an arbitrary value of σ, although the latter is frequently normalized to ±1 in the model (2.1). We let σ be an arbitrary constant in the solutions because this can help inspect the transitionary behavior of solutions, if it exists, in the vicinity of σ = 0. In this review, we shall follow such a convention to present solutions, unless otherwise stated.
Among the above three types of solitons built on a finite background, the Peregrine soliton has the peculiarity of being localized in both space and time, with a peak amplitude which can reach three times the surrounding background height [38], as seen in figure 1(c). It thus depicts a wave that appears from nowhere and disappears without a trace [18], matching the fleeting nature of oceanic rogue waves frequently witnessed by seafarers [3,10,12]. For this reason, the Peregrine soliton has often been thought of as the prototype of rogue waves [39], as mentioned in section 1.1.
Later, with the great success of the Peregrine soliton in modelling realistic rogue waves, researchers began to seek higher-order (multiple) rogue wave solutions to the NLS equation [40][41][42][43][44], using either the Darboux transformation (DT) [45] or Hirota bilinear method [46]. Almost in the same breath, super rogue waves with order up to the fifth have been observed in a water wave tank [47][48][49]. Here we do not present the analytical results of higher-order rogue waves for scalar systems but leave this interesting issue in section 3 for the analysis of vector rogue waves. On the other hand, there was also an intensive study on controllable rogue waves in dispersion-and nonlinearity-managed systems whose model coefficients can be made variable [50,51]. Basically, such types of rogue waves, usually constrained by stringent parameter conditions, can often be constructed from the known solutions of the constant-coefficient equations by means of self-similar symmetry reduction. Hence, in this review, we are primarily concerned with the unmanaged rogue wave systems, with scope beyond the NLS framework.
In practice, the NLS equation is a little bit poetical to realistic problems, as there are also other significant physical effects that need to be considered, aside from the GVD and Kerr nonlinearity. For example, one needs to take into account the third-order dispersion (TOD), self-steepening, and/or delayed nonlinear response when studying the ultrashort pulse propagation in highly dispersive optical fibres [52,53]. Under these circumstances, the integrable Hirota equation [46] and Sasa-Satsuma equation [54] can gain more ground for application as compared to the simple NLS equation, although the former two usually involve higher-order terms with fixed proportions.

The Hirota equation.
Let us first consider the Hirota equation [46], which, in dimensionless form, can be written as The last two terms that enter with the real coefficients ε and γ are responsible for the TOD and a time-delay correction to the cubic term, respectively [55,56]. For the sake of discussion, here and in what follows we will adopt the denomination proposed in [57,58], namely selfsteepening, for the last term |u| 2 u τ , although in some works, it is the term (|u| 2 u) τ that was referred to as self-steepening [52,53]. Meanwhile, we will assume γ 0 without loss of generality. The parameter relation = γ/(6σ) is chosen so that equation (2.6) is integrable. It is easy to see that if γ = 0, equation (2.6) can reduce to the NLS equation (2.1). Noteworthily, in (2.6), we have taken into account both the self-focusing and self-defocusing scenarios which are specified by σ. As shown in [55,56,59], the Hirota equation is equivalent to the compatibility, R τ ξ = R ξτ , of the two Lax linear operators where R = [r, s] T (T means a matrix transpose), and with λ being the spectral parameter, and Here the asterisk denotes complex conjugation. Hence, using the Darboux dressing technique [55,56,60] based on the above Lax pair, one can readily derive the rogue wave solutions of the Hirota equation (2.6). We begin with the plane wave solution that now reads Seeding it into the Darboux formalism yields the fundamental rogue wave solution (2.14) A comparison between solutions (2.5) and (2.13) shows that the latter is none other than the Peregrine soliton, but with its peak being oriented along the line χξ − τ = 0. Meanwhile, it is clear that this Peregrine solution is only physically valid for the self-focusing nonlinearity, as in the NLS case.
15) where the coefficients γ, 4 , 5 etc are arbitrarily real, and with p 0 = 0 and p 1 = σu * . Using the functional derivative defined by equation (2.16) and the recursion relation (2.17), the first few K terms can be found to be Obviously, for these specified K terms, one can identify the first two terms in (2.15) as the NLS equation (2.1) and the first three terms as the Hirota equation (2.6). The first four terms are known as the Lakshmanan-Porsezian-Daniel equation [62] which has also been intensively investigated in recent years [63]. It is worth noting that, as in (2.1) and (2.6), both kinds of nonlinearity, labeled by σ, have been considered in (2.15). As regards this infinite NLS hierarchy, no matter how long it is, it strikingly admits the Peregrine soliton solution [61] which is similar in form to solutions (2.5) and (2.13), but with the defining parameters Now, the Peregrine soliton given by solution (2.22) propagates along the line χ ξ − τ = 0. It includes those defined by (2.5) and (2.13) as special cases. As pointed out in section 2.1.2, the whole infinite NLS hierarchy does not admit the Peregrine structure in the self-defocusing regime, as its solution becomes now singular. It is not surprising that these solutions are so similar in structure, as the whole infinite NLS hierarchy shares the same linear spectral problem constructed from a loop algebra of sl(2). However, not all scalar integrable wave equations of one field variable can be squeezed into this infinite NLS hierarchy. For example, the Sasa-Satsuma equation [54], although comprising the higher-order dispersion and nonlinearity terms as well, involves a 3 × 3 linear spectral problem resulting from the loop algebra of sl(3) [64], and hence admits more complicated rogue wave dynamics [65][66][67].

A special case: the Sasa-Satsuma equation.
The Sasa-Satsuma equation, named after its pioneers Sasa and Satsuma [54], has been an active area of research for the past two decades. In dimensionless form, it can be written as where = γ 6σ and ν = γ/2 so as to satisfy the integrability condition. Compared to the Hirota equation (2.6), an additional term proportional to the coefficient ν, sometimes referred to as self-frequency shift in literature [68,69], is included, which has its origin in the delayed Raman response [52,53]. Due to the introduction of this term, equation (2.27) possesses a 3 × 3 Lax pair and hence an unusual DT [64]. Proceeding with this DT, one can obtain the fundamental rogue wave solution [66]: where u 0 is the plane wave solution defined by (2.11) and (2.12), and In like manner, we have generalized the solutions to include the self-defocusing situation as well. We clearly see that this solution involves polynomials of fourth order, more complicated than the Peregrine soliton structure allowed by the Hirota equation, although both integrable equations differ by only one term.
To avoid the singularity, we have obtained from equation (2.31) the parameter condition β > 0 for existence of rogue waves, which can be simplified to give It follows that the rogue waves within the Sasa-Satsuma framework exist only in the focusing (σ > 0) situation, yielding the same conclusion as above, but however in a limited regime of ω [66], different from that discussed for the NLS hierarchy where the Peregrine rogue wave tends to exist in the whole regime of ω. Moreover, different from the Peregrine soliton structure in the NLS hierarchy, the rogue wave solution (2.28) involves one or two peaks, depending on the value of ω [66]. For illustration, we present in figure 2 a twisted rogue wave structure which features two peaks and four holes, as ω = −1/4 (the other two parameters have been fixed: a = 1, γ = 3). It can therefore not be identified as Peregrine soliton because the latter involves one peak and two side holes [17]. In [67], it was numerically demonstrated that this peculiar elongated two-peak structure can be excited from a chaotic wave field, hence endowed with a strong stability property, see figure 7 therein. In the latter simulation, another set of parameters a = 1, γ = −3, and ω = 1/4 was used, which could result in an identical structure to that shown in figure 2 by a time inversion τ → −τ.
2.1.5. Baseband MI as origin of rogue waves. Although the existence of rogue waves has been confirmed by multiple observations, uncertainty still remains on their fundamental origins. It is now generally recognized that the MI is among the several mechanisms that may lead to rogue wave excitation [70][71][72][73][74]. We know that MI is a fundamental property of many nonlinear dispersive systems, which is associated with the growth of periodic perturbations on an unstable cw background [7,24]. However, as revealed in recent works [75,76], not every kind of MI necessarily leads to rogue wave generation; it is generally the baseband MI that plays such a pivotal role. Simulations showed that in the passband MI region, only nonlinear small oscillations were observed. Here, by baseband MI we mean that the cw background undergoes instability in a region where the perturbations can have infinitesimally small frequencies. Conversely, the passband MI is referred to as the situation where the perturbation grows up in a spectral region that does not include Ω = 0 as a limiting case [75,76]. Below, we take the Hirota and Sasa-Satsuma equations as enlightening examples.
Suppose the plane wave solution (2.11) experiences small perturbations, given by where p and q are small amplitudes of the Fourier modes, Ω signifies the modulation frequency (Ω 0), and κ denotes the complex propagation parameter of perturbations. Substituting (2.37) into (2.27) and linearizing the resulting equation, we obtain a system of two coupled linear equations for p and q. This system has a nontrivial solution only when κ and Ω satisfy the following dispersion relation: where ν = 0 is given for the Hirota equation and ν = γ/2 for the Sasa-Satsuma one, with the same = γ 6σ . In principle, MI occurs whenever the root κ of the quadratic equation (2.38) has an imaginary part. Therefore, in the limit of Ω = 0 (that is, at a sufficiently low modulation frequency), the necessity of a complex κ requires that η 2 σ > a 2 ν 2 , under which there occurs the baseband MI that gives rise to rogue waves [75,76]. It is obvious that within either the Hirota or the Sasa-Satsuma framework, rogue waves exist only for self-focusing nonlinearity (σ > 0). Specially, when ν = γ/2, the above condition can exactly reduce to (2.36) obtained with the Sasa-Satsuma equation.
Solving equation (2.38) directly for κ and defining the MI growth rate as γ h = Ω|Im(κ)|, one can obtain the MI map versus Ω and ω. Figures 3(a) and (b) display such MI maps associated with the Hirota equation (2.6) and with the Sasa-Satsuma equation (2.27), respectively, both in the self-focusing (σ = 1) situation. It is clear that for the Hirota equation, the MI is of baseband type in the whole ω regime, while for the Sasa-Satsuma equation, it features baseband type as well, but exists in the regime ω < 1/γ − a/2 or ω > 1/γ + a/2. No passband MI was exhibited for both situations. These results are completely consistent with the parameter conditions for rogue wave existence that have been predicted in our analytical solutions (2.13) and (2.28).
In this topical review, we will further illustrate the equivalence between the rogue wave existence and the baseband MI by other convincing examples. We will show in particular that some vector nonlinear systems will also give birth to passband MI which, when acting alone, contributes little to the rogue wave excitation; only the baseband MI has such a pivotal role. Based on our observations, we believe that the baseband MI approach has a merit in that it enables us to predict extreme wave events not only for integrable systems, but also for nonintegrable physical models of great interest, whose analytical rogue wave solutions might usually be unknown [23].
2.1.6. Integrable turbulence and rogue waves. Nonlinear dynamical systems often exhibit behaviors characterized by erratic changes of their characteristic local parameters. These irregular behaviors are usually denoted as being turbulent [2]. The interaction of a multiplicity of waves can cause such turbulence. In integrable systems, these waves are basically solitons and breathers and their presence is fully determined by the initial conditions. Radiation waves are also modes of the system but their role is less important because they have low amplitudes.
Integrable turbulence has attracted much attention in recent years [77][78][79][80]. Among the most interesting phenomena observed in the turbulent wave fields is the appearance of rogue waves, those that have peak amplitude more than twice the SWH. The highest waves appear as a result of the multiple collisions of breathers and/or solitons [81] rather than as a single Peregrine soliton. When solitons dominate such collisions, we have 'soliton turbulence', while when breathers do the term 'breather turbulence' is then used [82]. The number of rogue waves, the PDF of the chaotic wave fields, and their physical spectra are all specific for either of these two situations. In general, a mixed regime will participate of the characteristics of both.
Solitons, breathers, and radiation waves are indeed deterministic in evolutions once the initial conditions are given, and they can be found using the inverse scattering transform (IST) technique. Let us take the focusing NLS equation, i.e. equation (2.1) with σ = 1, as a paradigm of integrable systems, which originates from the compatibility condition of the following set of linear matrix equations [29]: , as defined in (2.7), and the 2 × 2 matrices A and B are given by We should point out that the Lax pair given above is completely consistent with that defined by equations (2.7)-(2.10), of course in the NLS equation limit = 0. It is seen that the spectrum of eigenvalues of equation (2.39) does not depend on ξ. This means that the dynamics of the field u(ξ, τ ) at any ξ is defined by the spectrum of eigenvalues of (2.39) at ξ = 0, that is, can be fixed by the initial condition.
Recently, Akhmediev, Soto-Crespo, and Devine gave a clear evidence of soliton turbulence that plays an otherwise unperceived role in generating rogue waves, by solving the eigenvalue problem (2.39) numerically with given initial chaotic field inputs u(0, τ ) [82,83]. In their numerical calculations, a constant background field of unit amplitude perturbed by a random function was taken as initial conditions. This random function is bound to have zero mean value and be both Gaussian distributed and correlated. The resulting initial field intensity, normalized to mean-value unity, is then defined by two parameters: its standard deviation (μ) and its correlation length ( L c ) that provide good estimates of the mean height and width of the input waves, respectively. Besides, a large number of different realizations with the same initial statistical parameters have been considered so as to obtain the statistical characteristics of the resulting field after a certain propagation distance, at which some kind of stationary turbulence is reached [83]. It should be pointed out that, for un unperturbed constant background input, while the eigenvalues that are purely real correspond to radiation waves, the eigenvalues that are purely imaginary correspond to breathers (e.g. Peregrine solitons, KM solitons, or ABs) and those involving both real and imaginary parts may correspond to solitons. We will see that, when the initial noise level is increased, there is a strong excitation of solitons whatever happens in this competing process.
The sets of complex eigenvalues λ solved for three different values of μ are shown in figure 4. Figure 4(a) shows the eigenvalues found for µ = 0.1, which gives an initial u function close to the unperturbed cw. The upper limiting point λ = i in the spectrum corresponds to the Peregrine breather. The latter has a peak amplitude of 3 which is the maximum possible for breathers. The eigenvalues located below the point λ = i on the imaginary axis correspond to ABs that have lower amplitudes. Some eigenvalues have the imaginary part slightly exceeding 1, which may correspond to KM solitons. Clearly, the predominant presence of the cw component in the initial conditions leads to the excitation of different types of breathers through MI. Figure 4(b) shows the set of eigenvalues obtained when the chaotic perturbation is larger, µ = 0.5. Many of the eigenvalues correspond to ABs and some, those with both real and imaginary parts, to solitons. It is suggested that there appear a greater number of solitons as compared to breathers, giving rise to a higher probability of generating rogue waves. Larger deviations from the cw in the initial conditions may result in that practically no eigenvalues are left exactly on the imaginary axis, as seen in figure 4(c) where µ = 0.9. All eigenvalues now correspond to solitons rather than breathers. Their amplitudes are twice the imaginary part of the eigenvalue and they will acquire a velocity which is defined by the real part of the complex eigenvalue.
The increasing presence of solitons of higher amplitude and velocity in the chaotic wave field, which is true as μ grows, prophesies a higher probability of appearance of rogue waves. That is also fully confirmed by numerical simulations. Equation (2.1) can be easily solved numerically using a standard split-step Fourier method. By performing such a task for a high number of different realizations, we can construct the corresponding PDF. Figure 5 shows the PDF of the field intensity obtained after propagating a distance of 100 units along the ξ axis. The PDF curves are not settled at the initial stages of propagation, but converge to a nearly stationary regime and practically do not change after ξ = 20. These curves show clearly that the tails of the PDF increase considerably with the value of μ, in relationship with the relative content of solitons and breathers in the initial conditions.  [82]. In all cases L c = 0.79. As the eigenvalues appear in complex conjugate pairs, only the upper half of the complex plane is shown. The above results are universal and can be extended to other integrable models that have both solitons and breathers as essential parts of their complex dynamics, such as the TWRI equations (section 3.2), the vector NLS equations (section 3.3), and so on.

The general cubic-quintic NLS equation
At this stage, one may be curious whether the self-defocusing media could support the Peregrine soliton structure. To answer this question, let us consider another scalar equation framework [84,85] different from the NLS hierarchy, namely, In this cubic-quintic (CQ) type NLS equation, we still retain, for the sake of comparison, the same coefficients σ and γ for the Kerr nonlinearity and the pulse self-steepening effect, respectively, as defined above. However, to attain integrability, the combined coefficients (µ, γ) related to the last two terms-the nonlinearity dispersion and quintic nonlinearityare specially selected as shown. It is worth noting that in this integrable model, we have excluded the higher-order dispersion terms beyond the GVD. This will enable us to catch the significance of the self-steepening effect in the absence of higher-order dispersion. In optical contexts, such CQ NLS equation can model the propagation of ultrashort pulses in a singlemode optical fibre [52,53] or in quadratic nonlinear media considering the group-velocity mismatch [57]. The general scalar model (2.42) has several reductions under specific conditions. In addition to the standard NLS equation (2.1) corresponding to the case µ = γ = 0, this model can also reduce to many known derivative or modified NLS equations, e.g. the Chen-Lee-Liu (CLL) type NLS equation (µ = γ = 0) [86], the Kaup-Newell type NLS equation (µ = 2γ, γ = 0) [87], the Gerdjikov-Ivanov equation (µ = 0, γ = 0) [88], and the Kundu-Eckhaus (KE) equation (γ = 0, µ = 0) [84,89], all of which were recently explored for understanding rogue waves [90][91][92][93]. No doubt, once the rogue wave solutions of the CQ-NLS equation (2.42) were found, they should include those of the above equations as special cases.
Starting from the known solutions of the CLL equation followed by a gauge transformation, we obtain the general fundamental rogue wave solution of (2.42) [94] where the initial plane-wave seed u 0 and the extra phase Φ are defined as follows It is easy to see from (2.46) that if and only if η = σ − γω > 0 can this rational solution behave like a rogue wave. This condition implies that it becomes possible for rogue waves to exist in a defocusing regime as well, when the self-steepening effect comes into play, i.e. γ = 0. Besides, we find that, apart from the phase Φ, the complex polynomial function, given in the square brackets in equation (2.43), which characterizes the Peregrine structure, has an inherent phase denoted by Θ, (2.48) In contrast to the acquired phase Φ given by (2.45), however, this inherent phase Θ does not depend on the parameter μ and thus will be the same across different integrable systems for given γ and given initial parameters. As a special case, when γ approaches zero, equation (2.43) can boil down to which is none other than the first-order rational solution of the KE equation [93]. Correspondingly, the parameter condition for existence of rogue waves becomes now η = σ > 0, which suggests that in the absence of the self-steepening effect, the solution (2.49) is only valid for the focusing nonlinearity, as in the NLS hierarchy.

MI and the chirped Peregrine soliton.
Actually, the general existence condition of rogue waves, η = σ − γω > 0, can also be obtained via a baseband MI analysis [76]. We assume that the perturbed plane wave is given by (2.37) but with u 0 defined by (2.44). Substituting it into (2.42) and linearizing the resulting equation, we obtain a system of two coupled linear equations for p and q. This system has a nontrivial solution only when κ and Ω satisfy the following dispersion relation: which is a quadratic equation of κ. Obviously, in the limit of Ω = 0, one can obtain from (2.51) the parameter condition η > 0 for the rogue wave existence, a sufficient condition exactly consistent with our analytical prediction above. By use of the quadratic equation (2.51), we plot in figure 6 the MI maps of the background field, defined by for both the focusing and defocusing cases. It is seen that the MI is of baseband type in both cases and the Peregrine soliton states only exist in the ω < σ/γ regime for γ = 0, independently of the value of μ. Besides, for the same initial parameters, the MI in the defocusing case (see figure 6(b)) is weaker than in the focusing case (see figure 6(a)). Here we should point out that if η 0, equation (2.43) may also serve as a valid solution to (2.42), but generally manifests itself in the form of algebraic solitons [95].
It is now evident that the solution (2.43) exhibits a Peregrine structure, but will be relieved from any singularity because of the self-steepening effect, even in the defocusing nonlinearity situation. As such, it generalizes the original Peregrine soliton concept of the NLS hierarchy that only exists in the focusing situation [17]. Considering that this solution acquires an extra instantaneous frequency shift [53] we will dub it the chirped Peregrine soliton. It is interesting to note that the chirp δω is also doubly localized, as a rogue wave does, but on a zero background, markedly different from that of the traveling dissipative solitons which is usually of tanh-shape in the transversal dimension [96,97]. This chirp will have a central value of 8(µ − γ)a 2 , hence displaying a peak or a dip at the center, depending on whether µ > γ or < γ. Obviously, such kind of strong localization entails potentially important frequency chirping, for large (µ − γ) values, in both temporal and spatial domains. Figure 7 illustrates a typical chirped Peregrine soliton as well as its chirp characteristic in a self-defocusing medium. It is clear that for such a defocusing Kerr nonlinearity, there occurs a deterministic Peregrine soliton structure, usually involving an extra nonplanar phase (see figures 7(a) and (b)). More intriguingly, the chirp distribution will take a dark doublylocalized structure on a zero background for µ < γ, as seen in figure 7(c). In fact, as one can verify, even if the Kerr nonlinearity is zero (i.e. σ = 0), the Peregrine soliton state remains alive as long as η > 0.
We need to emphasize that for the current CQ-NLS equation, the chirped Peregrine soliton form found for the focusing case has no bearing upon the one in the defocusing case. In a physical sense, that a defocusing (or equivalently, normally dispersive) scalar system could allow a bright solution is not surprising. One may recall that in the CQ-type dissipative system, there also occur the bright dissipative solitons in both the anomalous and normal dispersion regimes [98,99].

Numerical simulations and experimental proposals.
Extensive numerical simulations were also performed to study the stability of chirped Peregrine solitons using the standard split-step Fourier method [82]. Particularly, in order to see whether the chirped Peregrine solitons would be easily generated in realistic conditions, we intend to excite them numerically by using the plane-wave solution (2.44) at ξ = 0 as the initial condition, perturbed by white noise of strength ε = 0.001. Typical simulation results are presented in figure 8, using otherwise identical parameters as in figure 7. In the amplitude evolution plot (see figure 8(a)), the first 5 distance units have been removed, as there are hardly any visible changes on the chosen scale. It is seen that after 10 distance units, the MI has developed completely and there would appear rogue waves, ABs, and KM solitons in the sea of waves. The part selected by a black rectangle is enlarged and plotted in figure 8(b). It is clear that the wave in the fore part of the 3D plot is exactly the Peregrine soliton shown in figure 7, as its profile agrees very well with the analytic solution (2.43) (see figure 8(c)). This result implies that the chirped Peregrine solitons, for either the focusing or defocusing nonlinearity, can be observed in a laboratorial environment as long as the model equation (2.42) applies.
In [94], several possible experimental settings for observation of chirped Peregrine solitons were put forward, in the context of nonlinear optics. Usually, one may consider the temporal dynamics of ultrashort pulses in soft glasses or organic polymers, which show observable quintic nonlinearity besides the cubic response [20,53]. In a bulk polydiacetylene paratoluene sulfonate (PTS) crystal, a stable propagation of 2D spatial solitons was reported, thanks to the stabilizing effect of a saturating quintic nonlinearity [100]. Hence, as one might expect, the generic chirped Peregrine solitons may also be observed in a 1D PTS waveguide, of course taking into account the self-steepening effect that always accompanies the ultrashort pulse propagation [53]. Additionally, one may consider the propagation of ultrashort optical pulses in quadratic crystals (e.g. periodically-poled lithium niobate or tantalate crystals) in the high phase-mismatch cascading regime, which may mimic effective Kerr and self-steepening effects [57,58,101]. Further, one might also consider mode-locked fibre lasers, since their distributed modelling involves quintic nonlinearity [21,22]. However, the latter is essentially of dissipative type, which breaks down the integrability required by the current analytic approach. Although beyond analyticity, dissipative rogue waves [14] should in general exhibit frequency chirping effects as discussed above. All these interesting issues are topics for future investigation among the soliton physics community.

2.2.
3. An extended case: the Fokas-Lenells equation. The significant role played by the self-steepening term in rogue wave evolution dynamics can be further understood within the following scalar wave equation Here ∂ τ is a partial derivative with respect to τ and again the coefficient γ accounts for the selfsteepening effect [57,58]. This equation is completely integrable [102], first derived by Fokas using the bi-Hamiltonian method [103] and then obtained by Lenells from the Maxwell's equations under appropriate envelop approximations [104], hence the name Fokas-Lenells equation [105]. It differs from the CLL equation [86] by containing an extra space-time correction term, −(γ/σ)u ξτ , and thus can not be classified into the CQ-NLS framework (2.42). When γ = 0, the NLS equation (2.1) arises again. Using the standard DT method [95,106], one can obtain readily the fundamental rogue wave solution of (2.54): Once again due to the presence of the self-steepening term, the rogue wave solution (2.55) contains a complex denominator and can therefore exist in both the focusing and defocusing situations without singularity. It is of typical Peregrine structure, as seen in figure 2 in [95], but, however, there is no chirp being introduced to it, differently from the chirped Peregrine solitons discussed above.
As the function D in (2.56) should be positive definite everywhere, it follows that which is the parameter condition for rogue wave existence. Equivalently, one can solve the inequality (2.58) further and find that ω should lie between [95] σ γ and (2.59) Obviously, this condition is much tighter than the one obtained for existence of the rogue wave solution (2.43) to the CQ-NLS equation. Let us check the parameter condition (2.58) using the baseband MI analysis [76]. Substitution of the perturbed background field equation (2.37) into the governing equation (2.54) followed by linearization, we obtain the dispersion relation Apparently, in the limit of Ω = 0, the necessity of a complex κ requires that η 3 (ηγ 2 a 2 − σ 2 ) < 0, a condition exactly the same as (2.58). The MI maps determined by (2.60) are shown in figure 9, for both the focusing and defocusing situations. For better view, here we plotted the physical quantity ln(γ h ) versus Ω (> 0) and ω, where γ h is the growth rate as defined before. It is seen that the MI exhibits a baseband type in both cases, and for the same parameters, the MI in the defocusing situation is weaker than that in the focusing situation, as exhibited in figure 6. In particular, as Ω approaches zero, the baseband MI lies in a region exactly defined by (2.59). We point out that the results shown in figure 9 are totally consistent with those presented in [76], see for example figure 3 therein.

Vector rogue waves
While rogue wave investigation is flourishing in several fields of science, there is a necessity to go beyond the scalar NLS framework in order to model important classes of physical systems in a relevant way. One development of major importance consists in the study of coupled wave systems, as numerous physical systems comprise interacting wave components of different modes, frequencies, or polarizations. When compared to scalar dynamical systems, vector systems generally allow energy transfer between their additional degrees of freedom, which potentially yields families of intricate vector rogue wave solutions. In recent years, vector rogue wave dynamics were explored in diverse coupled integrable systems such as the LWSW resonant equations [107][108][109], the TWRI equations [110][111][112][113][114], the vector NLS (VNLS, sometimes alias Manakov system [115]) equations [116][117][118][119][120] or the Hirota extension [60], and others [121][122][123].
In this section, we will take the LWSW, TWRI, and VNLS equations as typical examples, to reveal the unique vector rogue wave dynamics. First of all, we will show that for coupled systems, there may appear a passband gain spectrum in the MI map, but it does not contribute to the rogue wave generation when acting alone, consistently with the baseband MI conjecture observed for scalar wave systems (see figures 3, 6 and 9). However, in coupled systems that involve more components, the baseband spectrum may extend significantly the domain such that the latter overlaps partially that of the passband spectrum. In this case, the passband MI may also be possible to dominate over the baseband one and excite the rogue waves. Secondly, we will show that there exist new rational solutions with dark structures that are also relevant to practical physical systems. In particular, in a defocusing VNLS system, we demonstrate vector dark rogue waves as well as their dark three-sister counterparts [75,119]. One may recall that in scalar NLS equations, no single dark rogue wave solutions take place, even in the case of a defocusing nonlinearity. Lastly, we wish to present other novel rogue wave phenomena such as coexisting rogue waves on the same background [109], WHL super rogue waves arising from three-wave interaction [113], complementary rogue waves originating in equal group-velocity propagation [114], and some other higher-order structures (e.g. rogue wave doublets, triplets, quartets, and sextets) [110,112,120].

Long-wave short-wave resonant system
Among coupled wave dynamics, the LWSW resonance is a general parametric process that manifests when the group velocity of the short (high-frequency) wave matches the phase velocity of the long (low-frequency) wave. It has been predicted in plasma physics [124] and nonlinear optics [125]. In hydrodynamics, the LWSW resonance can result from the interaction between capillary and gravity waves [126,127]. The fact that capillary rogue waves have been discovered experimentally [128] may stimulate further the investigation of rogue wave solutions for media manifesting LWSW resonance.
The nonlinear interaction of the complex short wave, u, and the real long wave, φ, is modeled by the following two coupled equations, expressed in a normalized form [125] (3.1) This equation is integrable, with its fundamental rogue wave solutions given by [107] is an arbitrary constant defining the background of the long-wave field). The parameters m and n in (3.2) are real, defined by Here ω n ≡ (2a 2 ) 1/3 , and the parameter has been given in piecewise form to avoid any ambiguity. It follows from (3.4) and (3.5) that only when ω 3ω n /2 do the above analytical solutions exhibit rogue wave structures.

MI and vector dark-bright rogue wave dynamics.
Let us first examine the MI of the background fields being perturbed in accordance with where p, q, and s are small amplitudes of the Fourier modes. As usual, a substitution of (3.7) into (3.1) followed by linearization yields the dispersion relation which is a cubic equation of κ. As the MI arises from a non-real value of κ, it is required that the discriminant of the cubic equation (3.8) should satisfy Hence, according to the baseband MI theory, we let Ω = 0 in the inequality (3.9) and obtain 4ω 3 < 27a 2 , which is just the parameter condition for rogue wave existence, exactly consistent with our analytical result obtained above. Figure 10 displays the MI map related to (3.8), where γ h = Ω|Im(κ)| denotes the growth rate. The dash-dotted lines denote the marginal instability given by ∆ = 0. It is seen that the MI map consists of two parts separated by a green dashed line-the baseband MI, which determines the domain of rogue waves as Ω approaches zero, and the passband MI, which does not have zero-frequency component and contributes little to the rogue wave formation. Clearly, we can see that not every MI is responsible for the rogue wave generation [75,76], which has nonetheless been overlooked in previous studies for years. Within the domain ω < 3ω n /2, the rogue wave structure will change with ω. As revealed in [107,129], the short-wave field exhibits a dark rogue wave state when ω (8a 2 /3) 1/3 ≈ 1.1ω n . Particularly, at ω = (8a 2 /3) 1/3 , a black rogue wave structure occurs, by which we mean the two side holes coalesce and the field amplitude can fall to zero at the dip center. We remind that in the scalar NLS equation, a self-defocusing Kerr nonlinearity results in a rational solution becoming singular, thus precluding a dark rogue wave solution [94,95]. Here, owing to the complex interplay between anomalous dispersion and the nonlinearity resulting from the coupling between u and φ fields, we anticipate that the contribution part of φ integrated along with propagation would wipe out the singularity, giving rise to an unusual dark rogue wave solution.
Field coupling provides a drastic change of the effective nonlinearity characterized by φ. It also favors strongly asymmetric solutions. If we consider a symmetric bright pulse formed at a given time, it has both temporally increasing and decreasing intensity parts, initially symmetric. One can infer from the governing equation (3.1) that, due to the coupling with φ, the increasing intensity part of the pulse will induce an increase of φ with propagation, as in a self-focusing process. Whereas the decreasing intensity part of the pulse will make a decrease of φ as in a self-defocusing process. This asymmetrizes the pulse in a peculiar way, giving rise to bright and dark localized pulsed solutions that depend on the magnitude of the normalized frequency ω/ω n .
Numerical results of the dark-bright rogue wave state at ω = (8a 2 /3) 1/3 are demonstrated in figure 11, with or without exerting white noise perturbation [107]. The initial pulse inputs are given by analytical solutions (3.2) at ξ = −5. We mimic the initial noise conditions by multiplying the real and imaginary parts of u and φ by [ where u and v are two complex short-wave field envelopes and φ is again the real long-wave field. Using the bilinear Hirota method [46], one obtains readily the exact fundamental rogue wave solutions [109] where the initial plane waves u 0 and v 0 are defined by their respective amplitude (a j ), wavenumber (k j ), and frequency (ω j ) ( j = 1, 2) according to The dispersion relations for the above seeding plane waves are given by k j = 1 2 ω 2 j − b. The real parameters m and n in (3.11) must satisfy (3.14) As a result, for given initial parameters, one can solve the algebraic equations (3.13) and (3.14) for the values of m and n, hence determining the solutions (3.11). It is easy to show that when one field vanishes (e.g. let a 2 = 0), these solutions can reduce to the ones (3.2), with m and n being explicitly given by (3.4) and (3.5).
For illustration, we plot in figure 12(a) the evolutions of m and n versus ω 2 that are allowed by (3.13) and (3.14), for given ω 1 = 0 and a 1 = a 2 = 1. It is interesting to note that, first, in such a multicomponent case, the rogue waves can exist in the whole domain, as m and n are allowed to have valid values everywhere, differently from the LWSW case that requires ω 3ω n /2; second, when ω 2 < 2.4600 (truncated to four decimal places), each nonzero ω 2 value yields two sets of valid (m, n) values which correspond to two different families of rogue waves. The former observation can be confirmed further using the baseband MI theory [75,76], which states that rogue waves could occur whenever the dispersion relation has a non-real root κ in the limit Ω = 0. This is always the case, as seen in figure 12(b) which shows the MI growth rate γ h = Ω|Im(κ)| versus Ω and ω 2 , obtained from the quintic equation (3.15) under the same parameter condition as in figure 12(a). Compared to that shown in figure 10, the baseband MI in figure 12(b), although quite involved, now extends to the whole domain, as expected from figure 12(a). Noteworthily, the MI map shown in figure 12(b) is made up of several baseband and passband parts which overlap partially. In the positive frequency region, the passband MI will also act a part in the rogue wave evolution. As regards the latter observation, it means that, for given initial background parameters, there may appear two coexisting rogue wave structures. For example, as ω 1 = 0 and ω 2 = −1.2469, equations (3.13) and (3.14) will give two valid values (m, n) = (−1.3514, 0.7803) and (−0.4287, 0.6442), specified by the red crosses and solid circles in figure 12(a), respectively. Substituting each pair of (m, n) into the solutions (3.11) gives rise to two different families of rogue wave structures as shown in figure 2 in [109]. This remarkable coexisting behaviors can be reminiscent of the bistable states occurred in soliton evolutions [130][131][132][133], although typically rogue waves are transients while solitons are stationary states. We then solved the underlying model equations (3.10) numerically to inspect the stability of these coexisting rogue wave families. A small amount of white noises is initially added to the analytical solutions (3.11) at ξ = −5, by multiplying the real and imaginary parts of both the u and v fields and the real φ field with the factor [1 + εr i (τ )] (i = 1, 2, · · · 5), respectively, where r i (τ ) are uncorrelated random functions uniformly distributed in the interval [−0.5, 0.5], and ε is a constant defining the noise level (here we used ε = 10 −4 ). The simulation results are illustrated in figure 13 (left and middle), whose (m, n) values are indicated in the caption, corresponding to the red crosses and solid circles in figure 12(a). It is revealed that these two coexisting rogue wave families are stable on an unstable background, despite the fact that the MI tends to interfere strongly with the trailing edge of the localized solutions after some propagation distance.
In order to see whether these two rogue wave solutions would coexist in realistic conditions, we intended to excite them numerically by using initial conditions significantly different from the exact solution profiles. To reduce the number of variables to play with, we used the plane-wave solutions (3.12) at ξ = 0 as initial conditions for u and v, and used φ(ξ = 0, τ ) = 0.4cos(2πτ /40)sech[(τ − 2)/8]. Surprisingly, for this set of initial conditions, we found it possible to excite both types of fundamental rogue waves on a background (see the right column in figure 13). It is clearly seen that, at around ξ = 10, there appear simultaneously two markedly different rogue wave types, well separated and corresponding to those shown in figures 13 (left and middle). This demonstrates further, from the numerical perspective, that the coexistence of diverse fundamental rogue waves is possible. The evolution that follows features additional multiple rogue wave dynamics, all with a combination of both rogue wave types. These subsequent multiple rogue waves could also be triggered by the onset of MI, which promotes quasi-periodic structures by patterning the cw background.

Three-wave resonant interaction system
As is known, TWRI enjoys a prominent status in nonlinear science (e.g. plasma physics, optics, fluid dynamics, and acoustics) [134]. In optical contexts, TWRI describes different processes such as parametric amplification, frequency conversion, transient stimulated Raman scattering (SRS) and backward or forward stimulated Brillouin scattering (SBS). As such, TWRI provides the basis for our understanding of diverse pattern-forming systems [134][135][136][137][138]. Other important domains of application of TWRI in nonlinear optics are group-velocity pulse control [139,140], ultrashort pulse train generation [141], laser-plasma interaction [142], and so on.
As early as 1970s, the integrability of the governing equations was established, and soliton solutions were identified [134]. These solitons are coherent localized structures that result from a dynamic balance between the energy exchanges due to the nonlinear interaction and the convection due to the group velocity mismatch [143]. This is in contrast to the case of quadratic solitons, where the energy flow among the waves is counterbalanced by GVD (or diffraction) [144]. Interestingly, TWRI solitons propagate with a common (or locked) velocity, despite the fact that the three waves travel with different linear group velocities before being mutually trapped [145,146]. This property makes such solitons very alluring in applications, since the walk-off caused by group-velocity mismatch, which usually limits the parametric frequency conversion efficiency, can be circumvented by nonlinear coupling. Moreover, when two optical waves are coupled to an acoustic wave via the SBS process, TWRI solitons may permit to considerably slow down the speed of light [139].
TWRIs can be conveniently classified according to the signs of the nonlinear coupling coefficients, as well as the ordering of the linear group velocities of each component wave. As discussed in [134], depending on these parameters, TWRIs may feature either soliton exchange (SE) dynamics (usually termed parametric three-wave mixing in the context of nonlinear optics) or stimulated backscattering (SB). TWRI may also exhibit an explosive behavior in that the coupled waves may develop a collapse in a finite time. Obviously, interactions of different types display very different behaviors. For instance, in the SE situation, velocitylocked solitons possess bright structures [145], whereas in the SB situation their dark counterparts would appear [147].
Recent works show that besides velocity-locked traveling solitons, TWRI equations also admit families of spatiotemporally localized solutions known as rogue waves [110][111][112][113][114]. In [110], fundamental rogue wave solutions were presented for SE-type interactions, modelling the sudden appearance of amplitude peaks in a basic multicomponent nonlinear wave system. In [113], the intriguing dynamics of watch-hand-like super-rogue waves were demonstrated to occur in such a three-wave mixing process. Once the group velocity of the two waves is equal, a pair of complementary rogue wave structures would appear [114]. Besides, it turns out that there exist rogue wave structures in the SB case as well, which behave differently from that in the SE case [112]. All these results will be systematically reviewed in this subsection.
3.2.1. Rogue waves in parametric three-wave mixing and coherent stimulated scattering. The TWRI equation that governs the propagation of three coupled waves in a weakly dispersive nonlinear medium can be written in dimensionless form [134] where u n (ξ, τ ) (n = 1, 2, 3) are the slowly varying complex envelopes of the three fields. From a physical standpoint, these fields may denote pump, signal and idler optical waves in the parametric mixing process that occurs in a quadratic medium [134][135][136], or describe the SBS (SRS) process where the optical pump wave scatters off a material acoustic (optical) phonon wave to form the Stokes wave [137,138]. The coefficients V n correspond to the relative group velocities of the three waves and we suppose V 1,2 > V 3 . Below we consider a reference frame comoving with u 3 , hence V 3 = 0. The above choice of signs before the quadratic terms is indicative of the nonexplosive character of the interaction [134]. Basically, as V 1 > V 2 , the interaction features the SE property, whereas the condition V 1 < V 2 corresponds to the SB process. In either situation, equation (3.16) admits exact rogue wave solutions because of its integrability and nonexplosive property [112]. Once V 1 = V 2 , one then comes to the degenerate TWRI equation, whose rogue wave solutions exist as well [114]. We will discuss the degenerate case in section 3.2.3.
Considering the resonant conditions for the frequencies and momenta, the initial planewave seeds that satisfy (3.16) can be expressed as where with a n (> 0) being the respective background heights. For convenience, we use with Γ j = Vj V1−V2 ( j = 1, 2 and the same below). Using the standard DT method [119,120], we obtain the fundamental (first-order) rogue wave solutions [110] u [1] where Here and for later use, we define the parameters β n as Once the parameters λ 0 and µ 0 are known, the general second-order rogue wave solutions can be cast in the following compact form [112] u [2] 1 = u 10 where with γ n (n = 1, 2, 3, 4) being four arbitrary complex constants. As the rogue wave structures strongly depend on γ n , we will term these parameters structural parameters in this paper so as to distinguish them from the plane-wave parameters. The specific spectral parameter λ 0 in the above formulas is the complex root (with a nonzero imaginary part) of the discriminant condition under which the cubic equation will have a double root µ 1 = µ 2 = µ 0 . Here (3.28) In general, equations (3.25) and (3.26) can be solved analytically, even with B = 0, although in most cases the root expressions are quite lengthy. For convenience of reference, we have provided these analytical solutions in appendix A. Usually, the assumption of B = 0 or a 2 = a 1 V 1 /V 2 may simplify the results greatly [112]. A close inspection of the cubic equation (3.26) as well as its discriminant (3.25) reveals that the rogue wave solutions for V 1 > V 2 could exist in the whole regime |δ| < +∞, while those for V 1 < V 2 only exist in the limited regime defined by These parameter conditions can also be obtained by use of the baseband MI analysis. To check this, we add small-amplitude Fourier modes to the plane-wave solutions and express where p n and q n are small amplitudes of the Fourier modes, and the parameters Ω and κ are assumed to be positive and complex, respectively. A substitution of these perturbed solutions into (3.16) followed by linearization yields the dispersion relation As a result of the baseband limit Ω = 0, equation (3.30) reduces to a real-coefficient quartic equation (B = 0) whose solution κ can be solved algebraically. It is easy to prove that as V 1 > V 2 (the SE case), this quartic equation always permits a pair of complex conjugate roots in the whole regime of δ, whereas for V 1 < V 2 (the SB case), it only supports complex roots when the parameter condition (3.29) is fulfilled, as predicted by the analytical solutions. On the other hand, one can solve the sextic equation (3.30) numerically for κ and get the MI gain map by defining it as γ h = Ω|Im(κ)|. It follows from (3.30) that the MI map will be same if the frequency difference δ has the same absolute value, for given amplitudes and relative velocities. We provide in figure 14 the MI map for the SB case. It is shown that the MI map in this case consists of baseband and passband parts, but only the baseband one determines the domain of rogue waves. We need to point out that the fundamental solutions (3.20) and the second-order solutions (3.23) can apply to both the SE and SB processes. To give some specific examples, we illustrate in figure 15 the dark-dark-bright (DDB) triplets for the SB case (left column) obtained with a 1 = 2, a 2 = 1, V 1 = 1, V 2 = 4, and δ = 2, and the bright-dark-bright (BDB) triplets for the SE case (middle column), obtained with a 1 = 1, a 2 = 2, V 1 = 4, V 2 = 1, and δ = 1.949. In both situations, we set κ = 0, as the rogue wave structures do not depend on κ , and use the same structural parameters γ 1 = 5, γ 2 = 1, and γ 3 = γ 4 = 0. It is evident that there is a distinctly different evolution dynamics between the SE and SB rogue waves. On the other side, in the right column of figure 15, we provided the numerical results of the DDB rogue  (3.30) in the SB case with a 1 = 2, a 2 = 1, V 1 = 1, and V 2 = 4 [112]. It has been separated into baseband and passband parts by a dashed white line. The green cross in the map indicates the maximum growth rate for |δ| = 2, which corresponds to a modulation frequency of Ω max = 2.41.
waves shown in the left column, under the same parameter conditions but imposing small amounts of white noise on initial profiles. As seen, in the presence of a tiny perturbation, the rogue wave triplets can still propagate very neatly for a rather long time, till eventually the spontaneous MI of the background fields grows up. Moreover, one can further compare the numerical simulations with the MI analysis, when the same initial plane-wave parameters are used. We note that in figure 14, when δ = 2, the maximum gain (green cross) corresponds to a modulation frequency of 2.41. Strikingly, in the right column in figure 15, it is shown that the period of the MI-induced waves is around 32/12, corresponding to a modulation frequency of 3π/4 2.36, almost the same as the former value. Once again a good consistency between numerics and theory is exhibited.

(3.34)
We should point out that, in addition to the six structural parameters γ n , the rogue wave solutions (3.32) only depend on the free parameters a 1 and on the relative group velocities V 1 and V 2 . Intriguingly, by choosing appropriate sets of structural parameters, the analytical solutions can exhibit highly asymmetric super rogue waves located around the origin, like the three hands of a watch. We now fall to discussing these novel watch-hand-like (WHL) super rogue waves, while leaving the discussion of other complex patterns such as doublets, quartets, and sextets in section 3.3.3.
Here we simply choose γ 3 = 1 (with all the other structural parameters set to zero), and let a 1 = 1, a 2 = 3, and a 3 = 2 √ 2, which implies that V 1 = 9 and V 2 = 1. Typical WHL rogue wave states are illustrated in figure 16, clearly showing that the three rogue wave components have a watch-hand-like distribution (see the contour plots (b), (d), and (f)). The watch trait will be more obvious if these three components are superimposed together in one image (see panel (g) where for comparison we have normalized the wave fields to have the same background height). We understand that the main spatiotemporal orientations of the WHL rogue wave components are determined by their relative group velocities. For instance, V 3 = 0 implies that, in the reference frame of u 3 , the hump orientation of the component u 3 is frozen at around ξ = 0 (see figure 16(f)). Also, as one can verify, increasing the value of V 2 towards V 1 would change the orientation of the other two field components, but still would maintain a watch trait, apart from its hands being extended significantly along the ξ dimension. Nevertheless, each rogue wave hand, although narrow in shape, always involves a highly asymmetric but still smooth profile. This unusual distribution is markedly different from the super rogue wave distributions that have been found to date in other coupled nonlinear systems [118][119][120], where the waves always feature a major overlap of their spatiotemporal distributions. In this regard, as they do not significantly overlap in time and space, the three components of such super rogue wave states could be efficiently separated using an appropriate filtering technique: this represents an important feature that would facilitate the experimental diagnostics and observation.
In addition to the above mentioned watch trait, all the three rogue wave components feature a giant wall-like hump, with a peak amplitude more than five times the respective background height. For illustration, we define an amplification factor g as the ratio of the peak amplitude to the average background. It is clearly seen in figures 16(a), (c) and (e) that, for the three rogue wave components, the corresponding g values are 5.43, 5.43 and 6.09, respectively. All of them are larger than 5, a characteristic value that the super (second-order) rogue waves in the scalar NLS equation can reach [18,47]. Here we point out that, for other parameter scenarios where γ 3 = 0, the g value is found to be smaller than 5, although usually a WHL distribution remains.
It should be mentioned that, under the parameter condition (3.31), the dispersion relation (3.30) resulting from MI analysis can be factorized into two cubic equations: These cubic equations can be exactly solved using the Cardano's formulas and hence the growth rate γ h = Ω|Im(κ)| can be readily calculated. Figure 16(h) shows the combined map of the MI gains in the plane (V 2 , Ω), calculated from (3.35) with given parameters a 1 = 1 and V 1 = 9. It is clear that the MI is of baseband type and as a result the WHL super rogue waves could exist for 0 < V 2 < V 1 . We finally simulated the governing TWRI equation and showed that these WHL rogue waves are themselves stable despite the onset of MI activated by small amounts of white noise (see figures 4 and 6 in [113]). It is even found that such WHL rogue waves can indeed be numerically excited from chaotic background fields. A typical result has been provided in figure 17, where we used otherwise identical parameters as in figure 16 except V 2 = 4. The noise strength we imposed on the initial background is now ε = 5 × 10 −4 , much larger than that used in figure 15. It is clearly seen that with this amount of noise perturbation, millions of waves were generated from MI after several propagation units, among them there appear some highest amplitudes that certainly bear a resemblance to (although due to the surrounding interference not so high as) the analytical ones, see those selected by black ellipses in figure 17. This no doubt reinforces the stability of these WHL rogue waves and our view on their potential applications.

The degenerate case: complementary rogue waves.
Since the rogue waves in the SE and SB cases behave differently, as seen in the left and middle columns in figure 15, one may naturally ask whether the TWRI equation (3.16) admits exact rogue wave solutions in the degenerate case V 2 = V 1 = V . The answer is of course affirmative, as revealed in [114]. By taking the limit V 2 → V 1 ≡ V in (3.20), one can obtain the exact fundamental rogue wave solutions u [1] 1 = u 10 where u n0 (n = 1, 2, 3) are specified as before by (3.17) with identical definitions for k 1 , k 2 , κ , δ and a 3 , but now A and B are defined by In the same way, the second-order rogue wave solutions can be found from (3.23), namely, where , (here and below j = 1, 2), Here γ n (n = 1, 2, 3, 4) are again four arbitrary complex constants. In what follows, we will assume ω 1 > ω 2 , i.e. δ > 0, without loss of generality. It is of interest to note that the sum of the intensities of the components u 1 and u 2 will always be conserved and can be given by independently of what parameters γ n we choose for deterministic rogue wave structures and whether or not their background heights are equal. As a matter of fact, the relation (3.40) is a natural consequence of the governing equation (3.16) which suggests that Therefore, if the field u 1 takes a bright Peregrine soliton [17,34] or a bright triplet [148] structure, the field u 2 will take the corresponding dark counterpart so that the intensity conservation (3.40) can be fulfilled. In other words, they are spatiotemporally complementary, hence the name complementary rogue waves. For illustration, we demonstrate in figure 18 the evolution dynamics of the second-order rogue waves with the same background height (see captions for specific parameter values).
Apparently, the former two butterfly-type patterns are spatiotemporally complementary (see figures 18(a) and (b)), while the third one, however, is always bright (see figure 18(c)).
It is easily concluded that the complementary rogue waves could exist in the whole parametric space. In fact, one can verify this using the baseband MI theory intended for rogue waves [75,76], which equates the existence regime of rogue waves with that of nonzero MI gain at arbitrarily low modulation frequency. We note that, in this degenerate case, the dispersion relation (3.30) can now be simplified to a quartic equation It is easy to find that, in the limit of Ω = 0, equation (3.41) always permits a pair of complex conjugate roots in the whole regime of δ (of course excluding δ = 0), as analytically predicted above. Figure 18(d) shows the MI gain map related to (3.41) in the plane (Ω, δ) for specific parameters a 1 = 1, a 2 = √ 3 and V = 4. It is clear that the baseband MI could extend over the whole range of δ, since the passband MI that occurs in coherent SB situation (refer to figure 14) is absent now. In particular, at δ = 1, the maximum of the growth rate γ h will occur at the modulation frequency Ω 1.23 (see the green cross in panel (d) or (e)).
Practically, of most concern to general soliton community is the stability of these rogue waves with respect to background broadband noise sources (e.g. quantum noise) [149,150]. Recent work showed that the Fermi-Pasta-Ulam recurrence of Akhmediev breather solutions of the scalar NLS equation may eventually break down in the presence of competing spontaneous noise-activated MI [149]. For this reason, we perturbed the initial deterministic rogue wave profiles by small amounts of white noise, and inspected whether the complementary rogue waves are still observed in the presence of the MI activated by such quantum noise. As in [112,113], we multiplied the real and imaginary parts of all three field components u n at sufficient negative times by a factor [1 + εr i (ξ)] (i = 1, . . . , 6), respectively, where r i are six uncorrelated random functions uniformly distributed in the interval [−1, 1] and ε is a small parameter defining the noise level. Figure 19 shows the numerical results intended for the fundamental rogue waves, either unperturbed (i.e. letting ε = 0) or perturbed by a small amount of white noise for which we choose ε = 10 −8 . It is clear that without any perturbations (see left column), numerical simulations produce almost identical results as predicted by solutions (3.36), hence giving evidence that the rogue wave solution is robust against numerical integration noise accumulation. On the other hand, in the right column, we showed that all three rogue wave components, under tiny perturbations, can still propagate very neatly for a rather long time, till eventually the MI of background fields grows up. Moreover, one can infer from the first figure in the right column that the period of the MI-induced periodic wave is around 18/3.5, corresponding to a modulation frequency of 7π/18 1.22, almost the same as calculated from the MI analysis, which has been indicated by the green crosses in figures 18(d) and (e).
Experimentally, we expect that there is a possibility to realize complementary rogue waves in a dual-mode optical fibre where a forward stimulated Brillouin scattering can occur [151]. In this case, the group velocity of the pump and Stokes optical waves can be nearly identical, while in comparison, the velocity of the acoustic wave is close to zero, resulting in an intermodal coupling governed by the degenerate TWRI equation. Under these circumstances, the complementary rogue wave dynamics would occur and one could observe in optical fibres two-color optical rogue waves of bright-dark type, thanks to the coupling with the acoustic wave and the relatively long interaction length.

Manakov system
Another coupled integrable model of great interest is the VNLS equation (also known as Manakov system [115]), which usually takes the following dimensionless form where u(ξ, τ ) and v(ξ, τ ) are the complex envelopes of the two field components and the sign symbol s indicates the anomalous dispersion (s = 1) or the normal dispersion (s = −1) regime. It is easy to show that by changing the variable ξ → −ξ, the above equation can be equivalently transformed to the focusing or defocusing VNLS equation [20]. We remark that for integrability condition, the ratio of the self-phase modulation to the cross-phase modulation in (3.42) should be equal to unity [111,115,152]. Albeit very special, such Manakov system has constituted a universal essential vector model for the fundamental exploration of complex coupled soliton dynamics (e.g. dark-dark and dark-bright soliton pairs) [153][154][155]. From a practical perspective, it can model pulse propagation in elliptically birefringent optical fibres [156] or the crossing sea waves in open ocean [157]. Particularly, there has been an experimental report of observation of Manakov spatial solitons in AlGaAs planar waveguides [158].
In recent years, the fundamental and higher-order rogue wave solutions of different kinds of VNLS equations, either focusing or defocusing, have been found [75,111,[116][117][118][119][120], with one or more field components exhibiting peculiar dark structures that are generally unattainable in the scalar NLS systems. Here, as in the TWRI case, we wish to present universal rogue wave solutions up to the second order applicable for both the anomalous and normal dispersion situations. We will highlight the dynamics of vector dark rogue waves and particularly the dynamics of vector dark three sisters that occur in the normal dispersion regime. We find that there is in fact an intimate relation between the solutions in the anomalous (normal) VNLS equation and those in the SE (SB) TWRI equation, if the field variables u 1 and u 2 of the latter are compared. and again let κ = ω 1 + ω 2 , and δ = ω 1 − ω 2 , as in section 3.2. Then, by virtue of DT, the fundamental rogue wave solutions of the VNLS equation (3.42) can be obtained as (3.21), but with β 2 = s(µ 0 + λ 0 − κ/2). Here j = 1, 2 and the same below. u 0 and v 0 represent the seeding plane-wave solutions (3.45) whose amplitudes (a j ), frequencies (ω j ), and wavenumbers (k j ) are connected by k j = s(A − ω 2 j /2). As one can verify, the complex parameters λ 0 and µ 0 in above formulas can be exactly determined by equations (3.25) and (3.26). Hence, for brevity, we do not present them here again and one can refer to appendix A for solutions.
Also, for given values of λ 0 and µ 0 , the second-order rogue wave solutions can be expressed in a compact form 48) and γ 1 , γ 2 , γ 3 , and γ 4 being four arbitrary complex constants. We point out that the fundamental and second-order rogue wave solutions given above are universal and can apply to both the anomalous and normal dispersion situations. Moreover, as one can see, these solutions possess rogue wave structures very similar to those defined by TWRI solutions [111]; this is evident if one compare the solutions (3.44) or (3.46) with the TWRI solutions (3.20) or (3.23). To be more specific, the rogue wave hierarchy allowed by the anomalous VNLS equation is almost the same as occurred in the SE-type TWRI equation, whereas that allowed by the normal VNLS equation bears a striking resemblance to the one occurred in the SB-type TWRI equation. A slight difference is that in the TWRI situation, the rogue wave structures do not depend on the parameter κ , but in the VNLS situation they definitely depend on it, which has been included in β 2 and hence in ϑ. Besides, there always exists a rogue wave hierarchy in the degenerate TWRI case, i.e. V 1 = V 2 = V [114], but, however, no rogue wave solutions exist in the zero dispersion (s = 0) limit of the VNLS equation.
Apart from the above observations, we find further that, due to the two-wave coupling, the normally dispersive (or equivalently, defocusing) VNLS equation can admit nonsingular rational solutions [75], as opposed to its scalar counterpart whose solutions are singular [94,95]. However, in the normal dispersion regime, the rogue waves only develop under the parameter condition [75,119] differently from those in the anomalous dispersion regime which can occur in the whole parameter space. This parameter condition is identical to (3.29), as it is a natural result of equations (3.25) and (3.26). One can also check it using the baseband MI theory [75,76], which requires, as usual, for given plane wave solutions (3.45) being perturbed according to should have a non-real root μ even in the limit Ω = 0. As usual, p n and q n (n = 1, 2) are small amplitudes of the Fourier modes. This is of course true, and for the sake of convenience, we give a detailed proof in appendix B. On the other hand, for a visualized comparison, we plot in figure 20 the MI maps related to (3.50), which is defined by γ h = Ω|slm(µ)|/2, for both the anomalous and normal dispersion situations (For simplicity, we only consider the B = 0 case so that the MI maps depend only on the absolute value of δ, as seen in (3.50)). It is shown that in the anomalous dispersion regime (see figure 20(a)), the MI map consists of two segments, both of which have a baseband gain spectrum that can occupy altogether the whole domain of δ. Particularly, one can find further that, although the passband gain spectrum does not determine the domain of existence of rogue waves, it may indeed contribute to the formation of rogue waves as long as it dominates over the baseband spectrum of other MI segments for given δ. By comparison, in the normal dispersion regime (see figure 20(b)), the MI map is composed of only one segment that comprises baseband and passband gain spectra. It is obvious that in this case, the baseband spectrum only occupies the domain |δ| δ m , which is exactly the same as the parameter condition (3.49) for rogue-wave existence. However, contrarily to the former case, the passband spectrum contributes little to the formation of rogue waves.
More interestingly, intriguing vector dark rogue wave dynamics can be observed in the normal dispersion regime, namely, the two polarization components can be kept being dark simultaneously. This behavior can be reminiscent of another class of closely-related nonlinear entities, vector dark solitons, which form as a result of a strong nonlinear coupling balanced by the normal dispersion [153,155]. As one might know, such vector dark rogue wave states are absent in most focusing coupled systems [107,109,116,118]. In the following, we are merely concerned about the normal dispersion (defocusing) regime where the vector dark rogue wave dynamics can occur.

Vector dark rogue waves and recent experimental progress.
For simplicity, we consider the case of equal background height, a 1 = a 2 = a (i.e. B = 0), and let s = −1 (i.e. normal dispersion). In this case, equations (3.25) and (3.26) can be exactly solved, with results given by (see appendix A) where now A = −2a 2 (< 0) and the parameter η is defined by For this special case, it follows easily that as −3A/2 |δ| < δ m , both field components could feature a dark structure. Particularly, as |δ| = −3A/2, the rogue-wave intensity can exactly fall to zero at the dip center. We will dub such special dark rogue waves as being black, similar to the definition of the black soliton terminology in contemporary soliton science [159]. This is nicely illustrated in figure 21, where we demonstrate both a vector dark Peregrine soliton given by (3.44) (see surface plots (a) and (b)) and a vector dark triplet given by (3.46) (see surface plots (c) and (d)), using the same set of initial plane-wave parameters a 1 = a 2 = a = 1, δ = − √ 3, and κ = 1/2 (the values of γ n for the triplet pattern are specified in the caption).
Concerning the intriguing triplet pattern shown in figures 21(c) and (d), one might be reminded of that there is a similar three-wave phenomenon occurring in Great Lakes known as the 'three sisters' [160], wherein the first wave will hit the ship and before its water drains away the second wave strikes, followed by the third wave. This sequential impact adds together suddenly overloading the deck with tons of water. These types of waves are now believed responsible for the sinking of the SS Edmund Fitzgerald on Lake Superior in 1975. Our vector dark triplets shown here can be thought of as the dark counterparts of these three-sister rogue waves. For this reason, we vividly term them 'dark three sisters'. Numerical results of the vector dark triplets without white noise perturbation were demonstrated in figures 22(a) and (b), using the analytical solutions (3.46) at ξ = −20 as initial conditions. It is clearly seen that the numerical solutions exactly reproduce the analytical solutions at least until ξ = 20, anticipating the stability of these dark triplets. We further numerically inspect the stability by perturbing the above initial condition. Specifically, we multiplied the real and imaginary parts of both the u and v components at ξ = −20 by a factor [1 + εr i (τ )] (i = 1, · · · , 4), respectively, where r i are four uncorrelated random functions uniformly distributed in the interval [−1, 1] and ε is a small parameter defining the noise level (here we used ε = 10 −4 ). Numerical results are provided in figures 22(c) and (d), showing that the former two sisters in either component can evolve as before, but the third one seems to be strongly disturbed by a complicated set of wave structures arising from the MI. This is not surprising because the MI can grow exponentially with the propagation distance and eventually it will develop to a large value to form waves interfering with the third sister. In other words, the dark three sisters themselves are stable, although the backgrounds where they are built are always unstable. Moreover, we note from figure 22(d) that the period of the MI wave is around 5, corresponding to a modulation frequency 2π/5 1.2566, a value very close to that given by the MI analysis above, see red cross in figure 20(b).
The experimental observation of an optical spatiotemporal dark rogue wave has been recently reported by using standard telecommunication equipment [161]. In this experiment, two orthogonally polarized cw pumps, provided by diode laser pumps, were injected in a 3 km long standard telecommunication fibre in normal chromatic dispersion and low polarization mode dispersion regime. Under these conditions, the nonlinear spatiotemporal propagation effects in the telecommunication fibre can be well described by the Manakov system [162]. When the total pump power injected in the fibre grows larger than a critical value, the band of unstable sideband frequency shifts extends down to zero, reaching the baseband MI condition [162]. It turns out that it is precisely such condition that admits the emergence of vector rogue waves. In addition, to excite a spatiotemporal dark rogue wave an input modulation is necessary to break the stability of the orthogonal pumps. When the modulating perturbation signals are present, the nonlinear reshaping of the polarized pumps leads to the generation of a periodic train of temporal black intensity notches, as shown in the numerical intensity plots in figure 23. Each of the generated notches matches well the shape of a single isolated optical dark rogue wave defined by the exact solutions (3.44), expressed in dimensional units.
The numerical and theoretical dark rogue wave predictions were confirmed by the experimental results. Figure 24 compares input and output (after 3 km of optical fibre) intensities from the experiments with an input periodic intensity modulation, and as well as their corresponding numerical and analytical dark rogue wave solutions. As can be seen, an overall excellent quantitative agreement is obtained between theory, numerics and experiments. Only slight discrepancies appear due to the non-ideal initial conditions used for the generation of the optical rogue waves.

The rogue wave doublets, quartets, and sextets.
Let us now give a final remark on the rogue wave solutions of Manakov system. We find that there exists a class of more general solutions when the parameter conditions B = 0 and A = 2δ 2 are fulfilled, as in the TWRI case. It is easy to show that, under such conditions, which mean that s > 0 (i.e. anomalous dispersion) and √ s|δ| = a 1 = a 2 ≡ a, (3.53) the cubic equation (3.26) will have a triple root µ 1 = µ 2 = µ 3 = µ 0 = 0, and the spectral parameter satisfying the discriminant (3.25) is found to be λ 0 = ±i √ 3δ/2. As a result, one can generalize the second-order rogue wave solutions (3.46) as [120] u [2] where u 0 and v 0 are still given by (3.45), with γ n (n = 1, 2, · · · , 6) being arbitrary complex constants. The other functions are defined (see caption). It is clearly shown that the constituents in these doublets, quartets or sextets are similar, but need not be identical; some have a peak amplitude more than twice the background height, but some do not, a few even having a double-peak structure. Also the position of the constituents is not fixed, which depends strongly on the choice of the structural parameters involved.

Multidimensional rogue waves
Echoing the multidisciplinary diffusion of soliton concepts a few decades ago, rogue-wave research has continued its expansion, while providing novel perspectives for the manifestation of extreme waves in a variety of nonlinear media. This requires the study of the propagation models that go beyond the scalar wave equations (see section 2) or vector multicomponent systems (see section 3). From the physical point of view, the extension of rogue wave models to higher dimensions is essential [15,16,163,164]. Oceanic rogue waves are manifestly (2 + 1)D phenomena, whereas in the context of ultrafast optics, the propagation of intense light pulses in a nonlinear slab or bulk medium entails a complex multidimensional dynamics, where the spatial and temporal degrees of freedom can not be treated separately [165][166][167]. Further, the extension of nonlinear dynamics to a higher spatiotemporal dimensionality can never be considered as a trivial problem. For instance, a pure χ (3) Kerr medium, which is appropriate for transverse optical confinement in the 1D case, does not provide stable confinement for higher dimensions due to wave collapse [165]. This situation explains why the practical quest of spatiotemporal solitons in higher dimensions, also termed light bullets, has become a holy grail for nonlinear optics [168][169][170][171][172][173][174][175][176][177][178][179]. Concerning spatiotemporal dimensionality, analytical rogue wave investigations have been mostly confined to (1 + 1)D models so far, due to the difficulty in finding integrable models in higher dimensions. Typically, a natural extension of the scalar NLS equation (2.1) is the (2 + 1)D NLS equation: which has great potential in applications, e.g. providing a standard description for the propagation of optical pulses in planar waveguides [167,180]. However, this model equation is non-integrable and has no exact rogue wave solutions with the help of IST or DT. But this is not the end, just the beginning of the full story. To seek for an integrable fungible version of equation (4.1), several significant efforts were taken. The best known integrable example is the DS equation [6,181,182] iu z + 1 2 2) which models the evolution of a 2D wave packet on water of finite depth in the shallow water limit [183,184], or describes the propagation of acoustic-ion waves in a magnetized plasma [185]. The DS equation (4.2) can be divided into the DSI and DSII equations [6,181]. In the context of fluid dynamics, DSI corresponds to β = 1, which occurs for water depths where surface tension dominates, while DS II corresponds to β = −1, which occurs for water depths where surface tension can be neglected. In the infinite depth limit, however, the DS equation can be reduced back to the (2 + 1)D NLS equation (4.1) that is non-integrable. A major development in the understanding of DS equations came in 1988, with the discovery of a class of exponentially localized 2D solutions-'dromions' [186][187][188]. Lately, rogue wave solutions were discovered for both the DSI and DSII equations [189,190], revealing novel spatiotemporal rogue wave dynamics. We need to mention here that there were also research activities on the third-type DS equation [191][192][193] whose rogue wave solutions exist as well [194].
Recently, another integrable extension of equation (4.1) was derived, first from the basic hydrodynamic equations [195] and then from the dynamics of ion acoustic waves in a magnetized plasma [196], namely, In this extended (2 + 1)D NLS equation, while the second term accounts for the diffraction and dispersion effects, the third term has been justified to serve as an effective Kerr nonlinearity [195,196]. Owing to integrability, the corresponding rogue wave solutions of (4.3) have been explored in a recent work [197].
On the other side, by using the multiscale method, equation (4.1) can also be transformed to the KP equation [198], up to the first-order approximation [199][200][201]: (4.4) As an integrable (2 + 1)D extension of the KdV equation, the KP equation plays a fundamental role in nonlinear wave theory, allowing the formation of stable solitons or rational localized solutions pertinent to systems involving quadratic nonlinearity, weak dispersion, and slow transverse variations [202]. As in the DS equation (4.2), the KP equation has two distinct versions (KP-I for β = 1 [203][204][205][206] and KP-II for β = −1 [207,208]), which bear a similarity in form but differ significantly in their underlying mathematical structure and the corresponding solution dynamics. Particularly, intriguing original nonsingular lump solutions of the KP-I equation were highlighted in the literatures [203][204][205][206].
Besides the above typical circumstances, there are also other possible interesting integrable extensions of (4.1), for instance, taking the (2 + 1)D differential-integral form as in [209,210] to account for the nonlocal nonlinearity. For our present purposes, in this section, we wish to only overview those interesting high-dimensional rational solutions allowed by equations (4.2)-(4.4). For convenience, we will normalize the dispersion parameter β and the nonlinearity parameter σ in the above equations to be ±1, without loss of generality. Specifically, we begin with the rogue wave dynamics in the DS equations, and then make an attack on the higher-dimensional rogue waves defined by a composite (2 + 1)D model. The latter model is a combination of the Hirota equation and the complex mKdV equation expressed in different dimensions [59], which can yield the extended (2 + 1)D NLS equation (4.3) and the KP-I equation (4.4), respectively. Naturally, once the rogue wave solutions of such a composite model are obtained, they should also satisfy both (4.3) and (4.4). Finally, we will show further that the KP-I equation in fact admits a general hierarchy of nonsingular lump solutions, but only a limited number of solutions can be identified as rogue waves.

The Davey-Stewartson equation
Among the deterministic and stochastic models describing ocean waves as well as their destructive influence, the DS equation [183] plays an important role in the areas of fluid dynamics because of its integrability [181]. It results from a multiple-scale analysis of modulated nonlinear surface gravity waves propagating over a horizontal sea bed. In normalized forms, it can be expressed by the model equation (4.2). Here we should stress that, in oceanography, z represents the evolution time, x and t stand for the two-dimensional space coordinates, and correspondingly, u and φ give the meanings of the surface wave elevation and the particle velocity along the dominant wave direction, respectively [6]. It is now known that this model equation has many identifiable coherent structures and waves, including solitons, dromions, Stokes waves, and velocity-dependent vortices [186][187][188]. Of most concern for our present purposes are those coherent structures referred to as rogue waves, which exist as well for DS equations [189,190]. To proceed, let us first inspect the existence condition of rogue waves following the baseband MI analysis.
The plane-wave solution of the DS equation (4.2) can be expressed as where k = a 2 σ − βω 2 /2 − κ 2 /2 − 2b, with b being the background of the φ field. Adding small-amplitude Fourier modes to the plane-wave solution (4.5) gives where p, q, and s are small complex amplitudes, and κ is a complex propagating constant while Ω and μ are real. Substitution of (4.6) into the DS equation (4.2) followed by linearization yields the dispersion relation: which is a quadratic equation of κ. In the baseband limit Ω = 0, this quadratic equation admits a pair of complex roots if and only if the discriminant satisfies According to the baseband MI conjecture repeatedly confirmed in previous sections, equation (4.8) will give the existence condition of rogue waves associated with DS equations. It is evident that in the DSI equation (β = 1), the inequality condition (4.8) can be always met for specific μ values independently of whether σ = 1 or −1, suggesting that there always exist rogue wave solutions in either case. However, in the DSII equation (β = −1), this condition holds only when σ = −1, which means that the rogue wave solutions exist solely for σ = −1.
As will be shown below, the σ = 1 case gives rise to singular rational solutions, which are physically irrelevant.
For further illustration, we solve (4.7) directly for the growth rate, which is defined as before by γ h = Ω|Im(κ)| (here we assume Ω > 0 without loss of generality). The MI maps for either β = 1 (DSI) or β = −1 (DSII) situation were shown in figure 26, under given parameter conditions κ = ω = 0 and a = 1. For better view, we plotted the logarithmic gain, i.e. ln(γ h ) in the DSI equation situation (see figures 26(a) and (b)), while in the DSII equation situation we used γ h for the MI map (see figure 26(c)). It is clear that for given constant background (κ = ω = 0) the MI that is responsible for the rogue wave generation can occur in the DSI equation, whatever sign the σ takes, but, however, only exists for the DSII equation with σ = −1. This conclusion can indeed apply to all background field scenarios, as confirmed by calculations.
On the other side, to seek for the rogue wave solutions, one can convert the DS equation (4.2) into the following Hirota bilinear forms by using the substitutions where F and G are functions of z, x and t (G is complex, but F should be real), and D is the Hirota bilinear differential operator defined by Then, in terms of the new variables, the Hirota bilinear forms (4.9) and (4.10) can be expressed as As the bilinear equation (4.14) is hard to solve, we can narrow the range of solutions by separating it into two parts Now equations (4.15)-(4.17) can be readily solved, with solutions given in terms of the determinants τ n that should satisfy τ * n = τ −n , for details see [41,189,190]. We do not repeat the derivations here but write the rational solutions as below.
The DSI equation case (β = 1). It is clear that both θ 1 and θ −1 are now real. As shown in [189], the nonsingular rational solutions can be constructed as with the matrix elements m (n) jk given by Here, n j are arbitrary non-negative integers and j and c jι are arbitrary complex constants. The shorthand notation ∂ j ( * j ) signifies the differentiation with respect to j (or * j ). As m (n) * jk = m (−n) kj , the condition τ * n = τ −n can be satisfied by the property that a matrix has the same determinant as its transpose. In appendix C of [190], the proof of the nonsingularity of the solutions, F = τ 0 > 0, was also given.
The DSII equation case (β = −1). In this case, θ −1 = σθ * 1 , and thus the nonsingular rational solutions need to be carefully constructed. As shown in [190], the solutions can still be given by (4.18), but the determinants τ n are made to be 2N × 2N order, namely, As in the DSI case, the constant parameters c jι , d kl , j , and k are all arbitrarily complex. With the help of the property of block matrices one can easily show that the determinants (4.22) fulfil τ * n = τ −n entirely. Particularly, as n = 0, one can find that only in the case of σ = −1 does the denominator F = τ 0 become positive definite everywhere, see appendix B in [190]. This suggests that the DSII equation admits the rational solutions only for σ = −1, differently from the DSI equation that admits rogue wave solutions for σ = ±1. This finding is completely consistent with our baseband MI analysis presented above, see figure 26.
Hence, the rogue wave solutions of the DS equation (4.2) are given by (4.11) and (4.18), with τ n given by (4.19) for the DSI equation and that given by (4.22) for the DSII equation, each followed by substitution of (4.13) into the determinants. A special note is that all these rational solutions are built on constant background fields (4.5), with k = κ = ω = 0. As demonstrated in [189,190], there are plenty of interesting rogue wave structures related to these solutions, for example, line rogue waves and even exploding rogue waves. We do not reproduce these results here and one can refer to [189,190] for more information.
Lastly, we wish to conclude this section with explicit fundamental rogue wave solutions obtained from the above general solutions. Let N = 1, n 1 = 1, and m 1 = 0 in (4.20), (4.23) and (4.24), one can get the fundamental rational solutions

(4.30)
Here ε is an arbitrary complex constant with nonzero real part. It is easily seen from these fundamental solutions that only when α > 0 do the rogue waves exist, which is sufficiently met in the DSI equation, but holds true only for σ = −1 in the DSII equation, once again confirming our baseband MI conjecture.

The composite (2 + 1)D wave equation and the rogue wave bullet concept
Let us now consider an otherwise naive combination of two 1D wave equations: which are expressed in different dimensional spaces, but with a common time variable t. It is seen that the first equation of (4.31) is the complex mKdV equation [211] and the second one the Hirota equation [46], both of which are integrable [56,212]. Here, ε is a constant parameter characterizing the perturbation to the NLS equation (we usually let > 0 below, unless otherwise stated). We recall that, if considered separately, either the Hirota or the complex mKdV equation has significant implications in modelling the propagation of ultrashort optical pulses in a general medium [213,214]. The physical significance of the composite model (4.31) can be viewed in a different light. On one side, by changing variables to Z = z − x, X = x, and T = t, and with the chain rule (4.32) which is none other than the extended (2 + 1)D NLS equation (4.3). It was reported that this single (2 + 1)D equation is suitable for modelling the capillary fluid waves [128,195], optical cavity waves [215], and ion acoustic waves [196]. On the other side, defining an intensity quantity E = |u| 2 , one can also obtain another (2 + 1)D integrable equation: which is nothing but the KP-I equation defined by (4.4), in addition to having different equation coefficients which can be made identical by rescaling variables. Since originally proposed for modelling the ion-acoustic waves in plasmas [198], this equation has now been obtained as a reduced model in hydrodynamics, ferromagnetics, nonlinear optics, Bose-Einstein condensates, and string theory [202,[216][217][218][219][220][221]. Therefore, once the rogue wave solutions of the composite equation (4.31) are obtained, they should also satisfy the above single (2 + 1)D wave equations (4.32) and (4.33), as the latter two can be directly derived from the former. The reverse may not be true, but this does not degrade the intrinsic interest of the combined model (4.31). We note that the integrable system (4.31) can be cast into a 2 × 2 linear eigenvalue problem, with the following Lax triad [59]: where R = [r, s] T (T means a matrix transpose), and with λ being the spectral parameter. It is easy to check that (4.31) can be exactly reproduced from the compatibility conditions based on the Lax triad (4.34).
In terms of the above Lax triad, a DT can be constructed, which is similar to that for the 1D scalar NLS equation. Hence, for the initial plane-wave seeds u 0 (z, x, t) = a exp (ikz + iκx + iωt) , (4.36) where κ = a 2 K − ω 2 ( ω + 1 2 ) and k = ω 3 − 6a 2 ω, with K = 6 ω + 1, one can obtain readily the fundamental (first-order) rogue wave solutions of (4.31): where Similarly, the general second-order rogue wave solutions u [2] can be given by where Here γ n (n = 1, 2, 3, 4) are arbitrary complex constants. We note that, apart from its propagation factor, the fundamental solution (4.37) can be expressed as a Peregrine soliton form [17,33,34] in terms of the combined variables ξ and τ, while the second-order solution (4.39) depends also on x besides ξ and τ. However, it should be stressed that these rogue wave solutions can not be factorized into a trivial product of two rogue waves obtained from the Hirota and the complex mKdV equations separately; they are totally inseparable. Moreover, both rogue wave solutions (4.37) and (4.39) can exactly satisfy (4.32) and (4.33) after a change of variables, as pointed out above. Let us take a close look at the MI of the background fields for such a composite system. By adding small-amplitude Fourier modes to the plane-wave solution (4.36), i.e.
]} (where p and q are small amplitudes of perturbations, Ω is the modulation frequency (Ω 0), and μ and κ are complex propagation parameters), one can obtain, after linearizing operations, two dispersion relations that μ, κ, and Ω should satisfy, It is obvious that in the baseband limit Ω = 0, both equations always involve non-real roots μ or κ for arbitrary ω (of course excluding ω = − 1 6 or 0), suggesting that the MI will be of baseband type. This will be seen by MI maps in figure 27, where we define the MI growth rates by γ x = Ω|Im(µ)| and γ z = Ω|Im(κ)|, Therefore, it is easily concluded that the rogue wave states would occur in the whole domain of ω. Specially, as ω = −(6 ) −1 or 0, they can be manifested as the line rogue wave forms in the plane (t, x) or (t, z) [189,190]. The typical 3D intensity distributions of the fundamental and second-order rogue wave solutions are demonstrated in figure 28, each plotted in (t, z), (t, x), and ( x, z) planes, respectively. It is shown that the fundamental rogue wave in either plane takes the form of Peregrine soliton (see intensity plots (a)-(c)), while the second-order rogue wave could take a triplet structure [44,49,148] consisting of three Peregrine solitons as γ 2 = 0 (see intensity plots (e)-(g)). The isosurface plots presented in figures 28(d) and (h) reveal that both the fundamental and second-order rogue waves have a directional preference or a bullet nature by which we mean they can propagate in a certain direction ζ with transverse double localization (see figure 28(d)). For this reason, we prefer to term such a special 3D rogue wave behavior a rogue-wave bullet. Basically, the propagation direction ζ of Peregrine rogue-wave bullet can be determined by the spatial line which passes through the origin, as indicated by ξ = 0 and τ = 0 in (4.38). So is the case with the rogue wave triplet, for which its three rogue wave components can propagate almost along the direction vector ( K 6 , ω, a 2 + ω 2 2 ) of the spatial line (4.44), as seen in figure 28(h). However, it may become involved for other second-order rogue wave states, as their directional preference depends also on the choice of γ n . The robustness of these rogue-wave bullets has been confirmed numerically, in spite of the onset of spontaneous MI, see for example figure 4 in [59]. One may draw an analogy between these unusual rogue-wave-shaped bullets and the peculiar X-shaped light bullets that were discovered in 2-component quadratic media [222]. Besides departing from a conventional bell-shaped pulse, in both cases, the localization of the light bullet requires feeding from a substantial background field. Moreover, the combination of a rogue wave and a bullet-like propagation may also echo the question arose whether one could launch a rogue wave onto a target [223]. On the other hand, there is also a possibility to realize such kind of rogue wave bullets in laboratory, owing to their good stability against the spontaneous MI [59]. In higherdimensional media, the experimental complexity multiplies in both material fabrication and pulse characterization aspects. However, experimental advances have recently provided a strong incentive in the area of rogue wave investigation in complex media. As examples, the deterministic dark counterparts of optical rogue waves predicted for the (1 + 1)D VNLS equation [116] have just been observed in a telecommunication fibre [161]; the 3D optical rogue waves arising from inhomogeneity have also just recently been generated in a linear spatial system [224]. In the wake of these developments and particularly in view of the close link of the composite model (4.31) to the widely-used KP-I equation [199,217], we anticipate that the prediction of rogue wave bullets may spark significant experimental research on the generation of 3D rogue waves in nonlinear optical or hydrodynamic systems.

The Kadomtsev-Petviashvili I equation
We mentioned in section 4.2 that the rational solutions of the composite (2 + 1)D model (4.31) definitely satisfy the KP-I equation, but the reverse may not be true. Let us now clarify this last yet interesting issue. As a matter of fact, the search for the nonsingular rational solutions to the KP-I equation has been a subject of much interest and was extensively studied over the past decades [203][204][205][206]. However, these solutions are mostly built on the zero background. In a recent work [225], a family of nonsingular lump solutions that propagate on a finite background, with the order extended to infinity, were presented. It is shown that these compact solutions enable us to gain an insight into the complicated multi-lump spatiotemporal dynamics, as well as the rogue-wave bullet dynamics that they may exhibit under a certain parameter condition.
For this end, we write the KP-I equation as 45) where u(z, x, t) is a scalar function describing a nonlinear long wave of small amplitude propagating along the distance z, but with an asymmetric dependence on the coordinates t and x. It should be pointed out that the above interpretation of dependent variable u and of the coordinate variables can be adapted to the particular applicative contexts [199,202,217,218]. In addition, by rescaling variables, equation (4.45) can be transformed to other equivalent forms with different coefficients in front of each term.
Before proceeding, we should point out that equation (4.45) is invariant with respect to the Galilean transformation [226] t → χ ≡ t − 2ωx + 12ω 2 z, (4.46) x → ρ ≡ x − 12ωz, (4.47) z → z, and u(z, x, t) → u(z, ρ, χ), (4.48) where ω is an arbitrary real parameter. In other words, if u(z, x, t) is a solution of equation (4.45), so is u(z, ρ, χ), which can be obtained by simply replacing the original variables t and x with χ and ρ, respectively. As we will see, this property is very helpful for generalizing the special solutions. Besides, one can recall that the KP-I equation (4.45) can be cast into a linear eigenvalue problem, with the Lax pair defined through [227] iψ x + ψ tt + uψ = 0, ψ z + 4ψ ttt + 6uψ t + wψ = 0, (4.49) where u is the scattering potential, ψ is the eigenfunction, and w is defined by w t = 3u tt − 3iu x . Obviously, the compatibility of equations (4.49) requires that the potential u must satisfy the KP-I equation (4.45).
In terms of the above Lax pair, the binary DT [206,227] can be constructed to give the solutions. To be specific, letting w = 0 and the initial solution be u = a, one can find from (4.49) the corresponding eigenfunction ψ, where φ is a complex function of the spectral parameter κ. Here without loss of generality, we assume κ to be real, as any imaginary part of κ can be removed by the Galilean transformation defined above. Consequently, a repeated use of the binary DT gives the nth-order solution u (n) of (4.45) and the corresponding ψ (n) : Usually, there are two ways to define M, and hence M . One convenient choice is to use the differential form called Wronskian representation [228] det M = W(ψ 1 , ψ 2 , · · · ψ n ), (4.53) det M = W(ψ 1 , ψ 2 , · · · ψ n , ψ), (4.54) where ψ m (m = 1, 2, · · · n) are the linearly independent solutions of equations (4.49) for given seed u at specific κ = κ 0 , and W is the Wronskian defined by W = det A, with the matrix [227]. However, this approach may result in complex and singular solutions that are physically irrelevant [229], unless the functions ψ j are carefully constructed so that M is positive definite, for example, separating each ψ j into two according to equations (13) and (14) in [230]. The other efficient way is to use the integration form proposed in [206] for the matrix M that will be intrinsically positive definite. Usually, M has many possible constructions in terms of integrals [206]. Here we adopt the simplest integration form where with m = 1, 2, · · · , n and ψ given by (4.50). The dagger in (4.56) indicates the complex conjugate transpose and for the scalar function, ψ † m = ψ * m . Obviously, m will be a polynomial of degree m in t, x, z and satisfies the recursion relation  Here and in what follows, we suppressed the 0 from κ 0 for the sake of brevity. For convenience, the polynomials m up to m = 6 are provided below: 59) 4 = ξ 4 + 6ξ 2 ϑ + 4ξη + 3ϑ 2 + γ 4 , (4.60) 5 = ξ 5 + 10ξ 3 ϑ + 10ξ 2 η + 5ξ(3ϑ 2 + γ 4 ) + 10ϑη + γ 5 , (4.61) where Then, substituting (4.55) into the first equation of (4.52) followed by some manipulations, we obtain where F n is a positive polynomial of 2n degree, given by a sum of n + 1 absolute squares, (4.68) Here n k = n! k!(n−k)! is the binomial coefficient and = 2κ. The resultant general solution u (n) (z, x, t) follows after trivial replacement operations t → χ and x → ρ, with χ and ρ defined by (4.46) and (4.47). Note here that the polynomials defined by (4.68) will not be the sole rational expressions satisfying the KP-I equation (4.45) [206,231], as they are obtained from the specific M form (4.55).

Abnormal lump evolution dynamics.
Inspection on the evolution of the lump structures with z reveals that the lump components always tend to locate parallel to the χ = 0 direction at z = −∞ but along the ρ = 0 direction at z = +∞, independently of the choice of structural parameters. This abnormal transition is contrary to what might be normally expected in soliton collisions [20] where the horizontally propagating solitons could still propagate horizontally after a collision. We note also that the two-lump interaction demonstrated in [231] is a normal process where a t (or x)-aligned lump distribution is still t (or x)-aligned after the collision, distinctly different from the case under discussion. Moreover, we find that, despite the complexity of the lump patterns, the central location of each lump component can be accurately given by the real roots of the n-degree polynomial equation:  Usually, for the polynomial degree higher than three, one needs to resort to some routines to solve such an algebraic equation. However, compared with the lengthy polynomial (4.68), equation (4.69) provides a more efficient pathway to insight into the complicated multi-lump spatiotemporal dynamics.
To exemplify the spatiotemporal dynamics, we consider the sixth-order lump solutions, whose polynomial F 6 is given by [225]  One can refer to appendix C for other low-order polynomials. Figure 29(a) shows the isosurface plot of the sixth-order lump for arbitrary set of given parameters a = 1, = 1, ω = 0, γ 1 = 1, and γ 2 = −60, clearly indicating that as z increases, its six components can change from a vertical distribution (see figure 29(c)) into a horizontal distribution (see figure 29(d)), both with respect to the t axis. This intriguing evolution can be seen by the t-z view in figure 29(b), where the lump distribution tends to be perpendicular to the (t, z) plane at large negative z, but will be parallel to the (t, z) plane at large positive z, as suggested by the green arrows in figure 29(a). We should point out that for a nonzero ω value, the lump distribution at z = ±∞ is still horizontally or vertically aligned, but in the (χ, ρ) plane. If converting back to the original (t, x) plane, the lump distribution will be obliquely aligned parallel to the straight line t = 2ωx at z = −∞, but still horizontally aligned at z = +∞. This is a trivial result of the Galilean transformation invariance of the KP-I equation.

4.3.2.
Existence of higher-dimensional rogue waves. We note that each lump component exhibits a Peregrine-like structure [17,94], namely, a hump flanked by two dips on a nonzero background. However, contrarily to the genuine Peregrine soliton, the two dips of the lump are not necessarily to be of zero amplitude, and its peak amplitude does not have a fixed factor of the background height, both of which depend on the value of the free parameter ε [199]. One may now wonder if these rational nonsingular solitons could be identified as rogue waves. Let us now address this interesting issue below. We take the first-order lump solution as an example. Inserting the polynomial F 1 (refer to (C.1) in appendix C) into (4.67) followed by the replacement t → χ and x → ρ, we obtain the general explicit fundamental solution: where χ and ρ are again defined by (4.46) and (4.47). We note that in (4.71) the parameter γ 1 has been removed by the coordinate translation so that the lump peak is located at the center. It is clear that this lump solution involves a free parameter ε besides a and ω, and therefore involves a variable peak amplitude and two variable minimums, which are given respectively by u max = a + 4 2 , u min = a − 2 /2. Apparently, this feature does not match with the strict definition of the Peregrine soliton solution of the NLS equation, as the latter has a peak amplitude fixedly three times the background height and two side holes that can fall to zero in the dip center [17,94]. In this regard, the solution (4.71) resembles a soliton in dynamical structure more than a rogue wave, hence often termed the lump in the early papers [204,205].
However, a recent work [59] showed that the KP-I equation nevertheless allows rogue wave solutions as well. It suggests that the solution (4.71) may be closely related to the Peregrine soliton profile. In fact, by letting = √ 2a, we find that this fundamental solution can be expressed as the intensity of a typical Peregrine soliton, namely, u = |ϕ(z, x, t)| 2 , (4.73) where ϕ(z, x, t) is the (2 + 1)D Peregrine soliton given by with the dispersion relations κ = a − ω 2 and k = 3aω − ω 3 . It is not surprising because the latter is the exact solution of the following complex mKdV and NLS equations: which can give birth to the KP-I equation (4.45) through the relation (4.73) [59]. Figure 30 illustrates this special fundamental lump structure ( = √ 2a) at z = 0, with either ω = 0 or ω = 1/2. It is clear that the lump structure can take the shape of the typical Peregrine soliton intensity (see figures 30(a) and (b)), but exhibit an extra traveling wave behavior, as seen in figures 30(c) and (d). In [59], we termed such a special structure a roguewave bullet in the context of nonlinear optics, by which we mean a nonlinear wave packet that has a characteristic rogue-wave profile in two dimensions, while it propagates without distortions in a third dimension, analogous to the X-shaped light bullet in a normally dispersive nonlinear medium [222].

Conclusion
In conclusion, we presented a timely and comprehensive review of rogue wave dynamics occurring in diverse physical systems ranging from scalar and vector integrable models to high-dimensional ones. We mainly focused on the analytical and numerical predictions of these rogue wave structures, as well as on their versatile evolution dynamics, while including reference to relevant up-to-date experimental results. The contents are multifarious, covering almost all the latest developments on this horizon. In this review, we attempted to organize these diverse results while clarifying several fundamental questions posed on rogue waves. They are summarized as below.
First, we addressed the controversial topic concerning the rogue wave origin and inquired whether a self-defocusing nonlinearity could allow for rogue wave solutions as well. For this end, we began with the simple scalar integrable models, which can be classified as two different frameworks, namely, the infinite NLS hierarchy (see section 2.1) and the general CQ NLS equation (see section 2.2), with both the focusing and defocusing nonlinearities being considered. In the former NLS hierarchy (e.g. NLS equation [33], Hirota equation [56], infinite NLS equation [61]), we showed by analytical examples that there occurs a typical Peregrine soliton structure with focusing nonlinearity, but no rogue wave structures exist in the defocusing regime. So is the case with the Sasa-Satsuma equation [65][66][67], although its rogue wave structures may be more complicated. However, in the latter CQ NLS framework, due to the presence of the self-steepening effect (which is not balanced out by the TOD), a special type of fundamental rogue wave of Peregrine-like structure could exist in both the focusing and defocusing regimes. As it usually acquires an instantaneous frequency shift (i.e. chirping) during evolution, we dubbed this type of rogue wave the chirped Peregrine soliton [94], analogous to the chirped soliton concept which was coined with respect to the soliton in nonlinear photonics. Why are there so different behaviors for two scalar models of seemingly minor difference (for instance, the Hirota equation as compared to the CLL equation)? Guided by this curiosity, we intended to explore the fundamental origin of rogue waves, and introduced two important topics recently developed. One is the analysis of the MI gain spectrum [75,76] that equates the rogue wave existence with that of the baseband MI as its frequency approaches zero (see section 2.1.5), the other is the integrable soliton turbulence [82,83] which has also been thought as an origin of rogue wave generation (see section 2.1.6). Both mechanisms are deeply addressed and confirmed by numerical simulations.
In the above discussion, there is an unheeded fact that in scalar models, even those with defocusing nonlinearity, no dark counterparts of Peregrine solitons exist. However we had clues to suppose the existence of dark rogue waves in coupled or vector integrable systems. To address this basic question, we explored the vector rogue wave dynamics within three multicomponent systems, i.e. the LWSW, TWRI, and VNLS equations. We showed both analytically and numerically that in vector systems, there definitely appear dark structures in field components. The reason is that in vector systems, the energy is able to transfer among additional degrees of freedom, hence leading to intricate vector rogue wave dynamics. Among those vector dynamics, the most interesting are: rogue wave coexistence [109], by which we mean different rogue waves structures could coexist simultaneously on the same background (see section 3.1.2); the super WHL rogue waves [113] that are endowed with three rogue wave components extending in different directions (see section 3.2.2); the complementary rogue waves [114] that occur as a result of optical wave interactions of equal group velocity (see section 3.2.3); and the vector dark three sisters [119], which were coined to reflect the dark characteristic of the realistic 'three sisters' phenomenon occurring in Great Lakes in North America (see section 3.3.2). All these dynamics are discussed in detail, and their existence is confirmed by numerical simulations as well as by the baseband MI analysis. In particular, we verified that in vector systems, the rogue wave existence condition is still well consistent with that of the baseband MI, although the MI maps become in these cases more complex due to the appearance of passband MI. Additionally, we also reviewed the recent experimental progress on observing vector dark fundamental rogue waves using standard telecommunication equipment [161,162].
Finally, in view of practical implications (e.g. considering that oceanic rogue waves are manifestly (2 + 1)D phenomena), one may wonder if the higher dimensional integrable systems support rogue wave solutions. Therefore, it is justified to consider the high-dimensional rogue waves, as a key ingredient of this review. For illustrative purposes, we merely discussed the three closely-related (2 + 1)D nonlinear systems, i.e. the DS equation (see section 4.1), the composite (2 + 1)D NLS equation (see section 4.2), and the KP-I equation (see section 4.3), by presenting the exact explicit nonsingular rational solutions. Novel spatiotemporal dynamics, such as linear rogue waves [189,190,194], lumps [203][204][205][206]225], and dark-lump waves [199], are unveiled. Particularly, in the study of a composite model, we highlighted the concept of rogue wave bullets [59] by combining concepts in nonlinear optics. We anticipate that the prediction of rogue wave bullets may spark experimental research on the generation of 3D rogue waves that has already been a topic among experimental physicists [224].
We should point out that this review is solely centered on the rogue wave dynamics in integrable systems, whether they are scalar, vector, or high-dimensional. However, from the experimental perspective, integrable models are too ideal to be real, particularly when additional factors must be considered in a laboratorial environment, which usually destroy the exact integrability [20,232]. In recent years, by awaking to the relevance and importance of nonintegrable models, there has been an increasing trend of rogue wave research on non-or near-integrable systems [23, 199-201, 233, 234]. It might be expected that this will become one of the future possible directions of theoretical research in the field of rogue waves. Actually, as we have highlighted in the review, the baseband MI analysis furnishes a route to predicting rogue waves in such nonintegrable systems. Special attention should be paid on that in order for g to be positive, (G + F) 1/3 + (G − F) 1/3 should be always kept real. Hence, this gives the parameter condition for rogue wave existence. One can verify that, while G is always real, F will be purely imaginary for A > 0, but real for A < 0 and for a limited domain of δ defined by It is worth noting that equation (A.8) is a natural result of the inequality H − 4K > 0 which enables F to be real when A < 0. Accordingly, as A > 0 (e.g. SE-type TWRI, focusing VNLS), one can get certainly that g > 0 in the whole domain of δ. However, in the case of A < 0 (e.g. SB-type TWRI, defocusing VNLS), g can only be positive under the above condition (A.8). Of course, in the latter case, one should take the real of the cubic roots of G + F and G − F so as to get a real g [119].
On the other hand, once the roots λ ± are known, one can substitute each of them into the cubic equation (A.2) and find the double root µ 0 . Let us define µ ± = λ 2 ± + δ 2 /12 + A/3. (A.9) Then, as one can check, the solutions (λ 0 , µ 0 ) that satisfy both (A.1) and (A.2) can be given by two of the four combinations (λ+, ±µ +) and (λ − , ±µ − ). However, in the special case B = 0, the above solutions can be greatly simplified. Specifically, these solutions can be expressed as below: It should be noted that, as A > 0, the value of η will be real for |δ| A/2, but it becomes complex when |δ| > A/2. However, as A < 0, η will always be real in the regime |δ| δ m ≡ 2 √ −A. Further, if B = 0 and A = 2δ 2 are met, one can readily get λ 0 = ±i √ 3δ/2, and the cubic equation (3.26) or (A.2) will have a triple root µ 0 = 0.

Appendix B. Baseband MI analysis of equation (3.50)
In the baseband limit Ω = 0, the quartic equation (3.50) can be rewritten as which is just the parameter condition (3.49) for rogue waves to develop. However, in the anomalous dispersion regime (s > 0), one can find that the discriminant (B.2) is always positive because of A > 0 and A 2 > B 2 . Meanwhile, among the two conditions r − p 2 /4 > 0 and p 0, at least one of them holds true, since, if A < 2δ 2 , the former condition is true, otherwise the latter condition is met. That is to say, in the anomalous dispersion regime, the rogue waves can exist in the whole parameter space.