Mathematical Modeling of Rogue Waves: A Survey of Recent and Emerging Mathematical Methods and Solutions

Abstract

Anomalous waves and rogue events are closely associated with irregularities and unexpected events occurring at various levels of physics, such as in optics, in oceans and in the atmosphere. Mathematical modeling of rogue waves is a highly active field of research, which has evolved over the last few decades into a specialized part of mathematical physics. The applications of the mathematical models for rogue events is directly relevant to technology development for the prediction of rogue ocean waves and for signal processing in quantum units. In this survey, a comprehensive perspective of the most recent developments of methods for representing rogue waves is given, along with discussion of the devised forms and solutions. The standard nonlinear Schrödinger equation, the Hirota equation, the MMT equation and other models are discussed and their properties highlighted. This survey shows that the most recent advancement in modeling rogue waves give models that can be used to establish methods for the prediction of rogue waves in open seas, which is important for the safety and activity of marine vessels and installations. The study further puts emphasis on the difference between the methods and how the resulting models form the basis for representing rogue waves in various forms, solitary or with a wave background. This review has also a pedagogic component directed towards students and interested non-experts and forms a complete survey of the most conventional and emerging methods published until recently.

1. Introduction

Anomalous waves, or “rogue waves”, represent a rare phenomenon at sea, which occurs on multiple occasions yearly [1,2] and causes millions of dollars of loss of cargo and loss of lives annually [3]. Rogue waves are abnormally elevated waves, with a height 2–3× that of the normal average wave and unusually having steep shapes [4,5]. Rogue waves were recorded for the first time in 1995 during a winter storm in the North Sea, when the “New Years Wave” hit the Draupner platform with a wave height of 27 m and 2.25× the average wave height [4]. The laser installation on the deck, which regularly records the elevation of the platform over the sea bed, registered the solitary giant wave with its elevation 15.4 m above and 11.6 m below the zero-level [4]. The shape of the wave was symmetrical (Figure 1) with a Gaussian bell shape and with a particularly narrow wavelength. This shape and behavior of anomalous waves are conserved across several observations made in the last 25 years, including the rogue wave that hit the North Alwyn platform in November 1997 [6], the Gorm platform in 1984 [6] and from Storm 172 on the North Alwyn field 100 miles east of Shetland [5]. The latter was particularly unusual, with a height 3.19× the average (Figure 1).

Rogue waves are known to have sunk over 20 supercarriers since 1970 [6] and carry a force of 16–20-times (100 metric tons/m${}^2$) that of a 12-m wave, and they can easily break ship structures, which are designed to withstand far lower impact forces (6 metric tons/m${}^2$) [3]. Rogue waves are an eminent threat to shipping and naval activities and increase in prevalence with climate-change weather patterns [7]. In this context, the insurance sector has searched for new models for predicting rogue waves and for fortifying naval structures [3], as both off-shore installations, shipping and also cruise-ships have been increasingly exposed to rogue waves in the last few decades [3,6]. This development has also sparked the project “Max Wave” [2], which has contributed with new models and algorithms for predicting rogue waves by the use of satellite observation data. Rogue waves occur also in optical systems [8], in the atmosphere [9], in plasma [10], as well as in molecular systems during chemical reactions [11].

Earlier mathematical models and derived algorithms that were used to predict wave patterns were originally developed by using the linear Gaussian random model, and rogue phenomena at sea were largely disregarded as superstition. The linear Gaussian model is essentially a superposition of elementary waves and predicts the occurrence of a rogue event at a very low probability. This low probability is however incorrect accounting for the laser readings made in the last two decades at off-shore installations. Non-linear models, which show a better agreement with the frequency of rogue events at sea, are therefore gradually replacing the Gaussian model used in the insurance industry. Non-linear models have been studied by several groups and include the modified non-linear Schrödinger equation (NLSE) [6], the Peregrine soliton model [12] the Levi–Civita and Nekrasov models [13,14], the Davey–Stewartson model [15], the fourth order partial differential equation of Kadomtsev–Petviashvili, the one-dimensional Korteweg–de Vries equation for shallow water surfaces, the second-order Zakharov partial differential equation [16] and the fully-nonlinear potential equations. Other systems have recently been developed and are herein reviewed in detail given their relevance to rogue wave ocean phenomena, including the inhomogenous non-linear Schrödinger equation [17], the Akhmediev model [18,19,20,21] and the recent models developed by Cousins and Sapsis [22,23,24].

2. The Non-Linear Schrödinger Equation in the Prediction of Rogue Waves

Rogue waves occur both in oceans, as well as in optical systems [12], as well as in other wave systems (see above). For fiber optic systems, rogue waves are normally entirely one-dimensional; however, two-dimensional rogue waves have been recently documented, in the form of the two-dimensional dissipative rogue waves [25]. These recently discovered optical rogue waves occur when a delayed feedback is generated in the transverse plane of the the cavity, forming an overlap of counter-directional fiber optic signals, which leads to a rogue amplitude [25]. These two-dimensional signals in optical systems are described by there own form of PDE, the Lugiato–Lefever equation [26], which allows for 2D rogue wave solutions to be modeled without collapse dynamics. This model is also used for describing a large spectrum of nonlinear phenomena in optical systems, such as bistability, localized structures, self-pulsating localized structures and also complex spatiotemporal behavior through an extended quasi-periodicity [27]. Rogue waves formed in fiber optic systems have also been recently considered as a new field of research in optics, given their anharmonic and nonlinear properties, which can be a future application of optical technologies [28]. In particular, a class of rogue waves which have the potential for application in optical technologies are the self-similar pulses [29,30,31,32]. Self-similar pulses are wave amplitudes measured in fiber amplifiers [33], which experience an optical gain together with a Kerr nonlinearity [34]. During the induction of the self-similar impulse in the solid, a fluid or any wave-carrying medium, the shape of the resulting rogue wave no longer depends on the shape or duration of the seed pulses, but depends only on the seed pulse energy (chirping). This creates a large effect on the amplitude, which is largely independent of the initial conditions of the wave pattern. This event, or formation of a rogue component in the wave train, has also been observed in ocean wave systems [35] and has attracted various groups to develop prediction methods using the variations of the non-linear Schrödinger equation (NLSE) [17,29,33,36,37]. One group in particular developed the variable coefficient inhomogenous nonlinear Schrödinger equation (vci-NLSE) for optical signals [17]:

$$i{\psi }_x+\frac{1}{2}\beta \left(x\right){\psi }_{tt}+\chi \left(x\right){\left|\psi \right|}^2\psi +\alpha \left(x\right)t^2\psi =i\gamma \left(x\right)\psi ,$$

(1)

which derives from the Zakharov equation [38]. Here, $\psi (t,x)$ is the complex function for the electrical (wave) field and x, and $\mathit{t}$ are respectively the propagation distance function and retarded time function. The parameter $\alpha \left(x\right)$ defines the normalized loss rate, and the function $\alpha \left(x\right)t^2$ accounts for the chirping effects (which indicates that the initial chirping parameter is the square of the normalized growth rate). The parameter $\beta \left(x\right)$ defines the group-velocity dispersion (i.e., for an entire wave train), while $\chi \left(x\right)$ defines non-linearity parameters, and $\gamma \left(x\right)$ defines the loss or gain effects of the wave signal. This equation is adaptable both for oceanic waves, as well as for optical non-linear wave guides. Equation (1) is essentially the same as the generalized Gross–Pitaevskii equation with the harmonic oscillator potentials in the Bose–Einstein condensates [39] and can be solved by applying the similarity transformation [40] by replacing $\psi (t,x)$ in Equation (1) with:

$$\psi (t,x)=\rho \left(x\right)\mathsf{\Psi }(T,X)e^{i\varphi (t,x)},$$

(2)

where $\rho \left(x\right)$ is the amplitude and T and X represent the differential functions describing the original propagation distance and the similarity variable, while $\varphi (t,x)$ is the linear variable function of the exponential term, which all must be considered well to avoid singularity of the system $\psi (t,x)$ [17]. T and X are given as:

$$T=\frac{t-t_c\left(x\right)}{W\left(x\right)}$$

(3)

$$X={\int }_0^z\frac{\beta \left(s\right)}{W^2\left(s\right)ds}$$

(4)

and hence, the similarity transformation gives:

$$i{\mathsf{\Psi }}_X+\frac{1}{2}{\mathsf{\Psi }}_{TT}+{\left|\mathsf{\Psi }\right|}^2\mathsf{\Psi }=0,$$

(5)

which is the standard non-linear Schrödinger equation.

The transformation and integrability conditions derived by [17] show that the factors of the wave system, such as effective wave propagation, distance, central position amplitude, the width and phase of the pulse, are ultimately dependent on the group velocity dispersion and on the non-linearity parameters of the system ($\alpha ,\beta ,\gamma ,\chi$). The “self-similar” solution found in the process of the transformation of the variable coefficient inhomogenous nonlinear Schrödinger equation into the standard nonlinear Schrödinger equation can ultimately be controlled under dispersion and non-linearity management [17]. Once transformed from the iNLSE, the solutions to the NLSE are derived by the derivation of polynomial conjugates to the root exponential function. This process is reviewed in detail here from the studies by [19].

The Solutions to the NLSE

The NLSE equation has been solved by various groups, including [16,19,40,41,42]. Following one of the most recent works by [18,19] in particular, the steps for deriving exact solutions to the NLSE are defined by identifying rational solutions [18] for the homogenous nonlinear system in Equation (5) by using the Darboux transformation [43]. This method is often used to derive rational solutions for non-linear systems and is adaptable to specific optical rogue waves, as well as ocean rogue waves, when represented by the NLSE. The main definition of a rogue event is that the wave “appears from nowhere and vanishes without a trace”, which is feature partly related to the behavior of solitons, which are independent waves that self-propagate and exit a collision unchanged. The origin of solitons arises from the first observation of a single solitary wave in the North Sea, made in 1834 by J. S. Russell, who later reproduced the solitary wave in a tank. Since then, solitons have been mainly studied in optical systems and are represented as solutions to several types of nonlinear PDEs, including the NLSE, the Korteweg–de Vries equation and the Sine–Gordon equation. This type of rogue behavior is described by the rational solutions derived from the NLSE [19], which describe an induction of a system instability to the top of a plane wave amplitude, which is transferred to the highest amplitude and then decays exponentially towards zero [18]. This behavior is represented by Ma-solitons and by Akhmediev breathers or “Akhmediev solitons” [18,19,44,45,46]. The difference between these two soliton models lies in the initial conditions, where the Ma-solitons originate from the initial conditions, while the Akhmediev solitons arise during evolution of the system given by modulation instability [44,47,48].

When solving the NLSE according to the Akhmediev scheme [19], their method describes the modeled envelope function ($\psi$) as a solution ranked into an order of hierarchy, starting from first and progressing to the second, third or fourth order [19]. The difference between each order is the increasing amplitude of the rogue wave (first order, lowest amplitude, fourth order, sharpest peak and highest amplitude). The envelope function $\psi$ is expressed as a ratio of polynomials multiplied to the complex exponential root function, $e^{ix}$. The polynomials, which are given by functions of variable x and t, are identified by performing the Darboux transformation on the NLSE system [19]. Akhmediev and colleagues furthermore applied a compatibility-check between the root function $e^{ix}$ and the reference-state for two specified column matrix elements, which define initial conditions for the NLSE. These matrix elements (vectors) are given specifically by Akhmediev and colleagues [19] as two differential equations:

$$r_x=i{l}^{2}r+il{\psi }^{\ast }s-il2\left|\psi \right|2r+\frac{1}{2}{\psi }^{\ast }s$$

(6)

$$s_x=i{l}^{2}s+il\psi r-\frac{1}{2}\psi tr+\frac{i}{2}\left|\psi \right|2s$$

(7)

which are split into real and imaginary parts, before being simplified and solved to fit into the modified Darboux scheme [19,43] to give the two linear differential forms:

$$r_l(x,t)=\sqrt{2}\left[xt-\frac{1}{2}+ix\right]e^{-ix/2}$$

(8)

$$s_l(x,t)=\sqrt{2}\left[x-i(t+\frac{1}{2})\right]e^{-ix/2}$$

(9)

where the two vectors (8), (9) are used in the Darboux scheme to find $\psi j$, where $\mathit{j}$ is the order of hierarchy. The general solution to the NLSE, derived from this scheme [19], is given by the following general form (for any order in the hierarchy):

$${\psi }_j(x,t)=\left[{(-1)}^j+\frac{G_j(x,t)+ix{H}_j(x,t)}{D_j(x,t)}\right]e^{ix}$$

(10)

where G, H and Dare the polynomials of the two variables x and t (mentioned above). The first order-solution [19] has the following polynomials: G = 1, H = 2 and $D=1+4t^2+4x^2$, which give the following envelope function (shown in Figure 2):

$${\psi }_l=\left[1-4\frac{1+2ix}{1+4t^2+4x^2}\right]e^{ix}$$

(11)

For the second-order solution, [19] identify the vectors r2 and s2 by solving Equations (6) and (7) by using the form of $\psi$ given in (10). This gives the second-order solution:

$${\psi }_2(x,t)=\left[{(-1)}^2+\frac{-2t^4-12t^2{x}^2+ix\left(8t^2{x}^2+2t^2-4x^4-2x^2+\frac{15}{4}\right)-3t^2-10x^4-9x^2+\frac{3}{8}}{\frac{16t^6}{3}+13t^2+\frac{16x^6}{3}+36x^4+33x^2+\frac{3}{4}}\right]e^{ix}$$

(12)

which is shown in Figure 3. The third and fourth order rational solutions are furthermore calculated and given in [19].

The same hierarchy-dependency is given in the approach by [17], for the transformed vci-NLSE, who define the general solutions for the NLSE in the hierarchical $\mathit{n}$-th order given by:

$${\psi }_n=\frac{1}{W}\sqrt{\frac{\beta }{\chi }}{[(-1)}^{-1}+\frac{G_n+i(Z-Z_0)H_n}{F_n}\times e^{i[(1-\frac{v^2}{2})(Z-Z_0)+vT+\varphi ]}]$$

(13)

where each factor is given for the first and second order rational solutions [17]. Similarly to the hierarchy solutions of Akhmediev [19], the increasing order gives higher and higher rogue waves, compared to their surrounding waves. The first and second order rational solutions given in [17] reflect respectively a 3× and 5× rogue wave height, compared to the surrounding waves. For plots of these, refer to [17].

The similarity between (13) and (10) is striking, and both retain the basic form of a complex polynomial multiplied by a complex exponential root function giving soliton solutions. The root functions of (10) and (13) are shown in their generic form in Figure 4, which depicts the distinction between the seed pulse for the regular NLSE and the vci-NLSE, as studied respectively by [17,19] for the rogue wave problem.

The root function for the vci-NLSE (Figure 4) shows its specific pattern of wave accumulation, which is similar to the formation of wave packets. This pattern is conserved with the physical behavior of rogue wave formation, where the rogue wave forms during a focusing phase [51]. Other approaches used to solve the NLSE have been given by [37], who used the inverse scattering method of transformation, which is a generalization of the Fourier analysis. Their solutions differ from the methods discussed above and are periodic and ascribed by a complex envelope function for the deep water train with added higher order terms from the perturbation procedure [37]. One of the solutions is shown in Figure 5, which shows the following variant of the Osborne models:

$$\psi =\frac{cos\left(\sqrt{2}x\right)sech\left(\sqrt{2}t\right)+i\sqrt{2}tanh\left(2t\right)}{\sqrt{2}-cos\left(\sqrt{2}x\right)sech\left(\sqrt{2}t\right)}{e}^{2it},$$

(14)

which is a periodic function in space, derived from the general form given in [37], shown in Figure 5.

The disadvantage of this system, compared to single-peak models derived from [19], lies in their periodicity and multiple peaks, while the rational solutions behind the single peak models of [18,19] are the first in general to serve as prototypes for rogue waves.

3. The Korteweg–de Vries Equation

Wave systems defined by higher order nonlinear PDEs, such as (2), can be solved also by the bilinearization technique [52]. This technique involves the step of transforming the differential equation into a more tractable form by replacing the unknown time- and position-dependent envelope function with a new form [52]. After this replacement has been performed, the bilinearization technique applies Hirota bilinear operators for a modified Bäcklund transformation technique [53], which assists in rewriting the original PDE into a simplified PDE composed of bilinear operators, from which exact soliton solutions can be identified. The most suitable example [52] for the application of the bilinearization technique is on the Korteweg–de Vries (KdV) equation:

$${\mathsf{\Psi }}_t+6\mathsf{\Psi }{\mathsf{\Psi }}_x+{\mathsf{\Psi }}_{xxx}=0$$

(15)

where the boundary conditions are that $\psi \to 0$ as $\left|x\right|\to \infty$. The real wavefunction is differentiated according to the spatial and temporal dimensions as denoted. In the bilinearization technique, a transformation of the wavefunction to another form is the first step, where an ideal steady-state form is proposed [52] to:

$$\psi (x,t)=(p^2/2)sec{h}^2(\eta /2),$$

(16)

where:

$$\eta =px-p^{3}t+{\eta }_0,$$

(17)

and ${\eta }_0$ and p are arbitrary constants. By the bilinearization technique [52], one can rewrite (16) in the form:

$$\psi (x,t)=2p^2(e^{\eta /2}+e^{-\eta /2}),$$

(18)

which is converted to its functional form:

$$\psi (x,t)=2\frac{{\partial }^{2}ln\left[f(x,t)\right]}{\partial x^2},$$

(19)

with f(x) = 1 + $e^{\eta }$.

The bilinearization technique [52] substitutes (19) into the original KdV Equation (15) and integrates it with respect to x:

$$f_{xt}f-f_x{f}_t+f_{xxxx}f-4f_{xxx}{f}_x+3{\left(f_{xx}\right)}^2=0$$

(20)

which is the original version of the bilinearized variant of the Korteweg–de Vries Equation (15) as derived by [53]. The solution to (20), f(x) = 1+$e^{\eta }$, is defined as a more fundamental quantity than $\psi$ in Equation (18) for the structure of the original nonlinear PDE in Equation (15). In the method of bilinearization, the Hirota bilinear operators are introduced. These are defined by the following definition [53]:

$$D_t^n{D}_x^{m}a·b={(\partial /\partial t-\partial /\partial t^{\prime })}^n{(\partial /\partial t-\partial /\partial t^{\prime })}^{m}a(x,t)b(x^{\prime },t^{\prime }){\mid }_{x=x^{\prime },t^{\prime }=t}$$

(21)

with $\mathit{m}$ and $\mathit{n}$ being arbitrary positive integers. At this stage, the converted form of the KdV Equation (20) is rewritten as a PDE composed of Hirota operators:

$$D_x(D_t+D_x^3)f·f=0,$$

(22)

which is a simplified form for the identification of exact solutions using the Bäcklund transformation for the original nonlinear PDE (15). The exact solution structure for the type of Hirota-operator based PDE form (22) of the KdV Equation (15) is given by:

$$\mathsf{\Psi }=1+ϵ(e^{{\eta }_1}+e^{{\eta }_2})-{ϵ}^2\frac{F({\mathsf{\Omega }}_1-{\mathsf{\Omega }}_2,p_1-p_2)}{F({\mathsf{\Omega }}_1+{\mathsf{\Omega }}_2,p_1+p_2)}{e}^{{\eta }_1+{\eta }_2}$$

(23)

which represents the two-soliton solution to the original KdV Equation (15). ${\eta }_1$ and ${\eta }_2$ are the functions with the independent variables x and t as given in (16) for each of the solitons, and ${\mathsf{\Omega }}_1=-p_1^3$ and ${\mathsf{\Omega }}_2=-p_2^3$, following the same definition for (16) for each soliton. $\eta$ represent the perturbations [52]. The KdV Equation (15) has also been solved by Matveev by identifying positon solutions [54], which exert the same behavior as solitons, such as conserved shape after collision and elastic collision behavior. The positon differs from the soliton in that it has an infinite energy and is therefore not a strong model for oceanic or optical rogue waves. Positons have however a tendency to represent smoother solutions than solitons to the KdV equation and can have very high peaks compared to the wave normal. The KdW equation has also been solved by a nonlinear Fourier method [55,56], which is represented by a superposition of nonlinear oscillatory modes of the wave spectrum. This model, developed by Osborne [55,56], has the capacity to include a large number of non-linear oscillatory patterns, also known as multi-quasi-cnoidal waves, which are used to form the rogue wave by superposition in constructive phases. These solutions to the original KdV Equation (15) include several solitons, depending on the number of degrees of freedom selected for the numerical simulation of the KdV equation. This yields a 3D wave complex composed of solitons and radiation components in the simulated wave train [55].

4. The Extended Dysthe Equation

In 1979, Dysthe [57] developed a modification of the perturbation-based NLSE by adding an additional term to the third-order perturbation variant originally developed by Higgins [58]. Dysthe’s method gave an NLSE variant, known as the extended Dysthe equation, which was shown to have a better agreement with the mean flow response to non-uniformities in deep-water waves. The extended Dysthe equation is given by:

$$\frac{i}{k}{\psi }_{xyy}+{\psi }_{yy}+2ik{\psi }_x+2{\psi }_z=o{ϵ}^4,$$

(24)

where the inhomogenous component is the fourth-order perturbation defined by Dysthe [57]. Dysthe transformed this equation to standard NLSE using dimensionless variables and added the following perturbation to the general solution:

$$\psi =c_0(1+\alpha )e^{i({\theta }^{\prime }-\frac{1}{2}i{c}_0^{2}t)}$$

(25)

where $\alpha$ and $\theta$ are small real perturbations of the amplitude and phase, respectively. After insertion of (25) in the dimensionless form of (24) and linearizing, Dysthe obtained a simplified system of two PDEs, where the respective plane-wave solutions are in the form:

$$\left(\begin{array}{c}{\alpha }^{\prime }\\ {\theta }^{\prime }\end{array}\right)\propto e^{i(\lambda x+\mu y-\mathsf{\Omega }t)}$$

(26)

and:

$$\overline{\varphi }\propto e^{[Kz+i(\lambda x+\mu y-\mathsf{\Omega }t\left)\right]}$$

(27)

where K = $\sqrt{({\lambda }^2+{\mu }^2)}$ and $\lambda$, $\mu$ and $\mathsf{\Omega }$ are selected parameters that satisfy a set of dispersion relations given by Dysthe [57].

The stability of the solutions derived by Dysthe shows that the Dysthe equation represents a more realistic model than the NLSE, given that it does not predict a maximum growth rate for all wavevectors, but only for some wavevectors. This displays that the fourth order perturbation term added to the NLSE gives a considerable improvement to the results relating to the stability of the finite amplitude wave. It is particularly the first derivative by the transformed variables in the x and z dimensions in Equation (24) that contributes to the excellent results of Dysthe. Dysthe and Trulsen [59,60] further developed this equation by including up to the fifth-order of the derivative of the wave amplitude describing the linear dispersive terms and simulated successfully [61] the New Year’s wave [4] using the extended Dysthe equation [57,61].

5. The MMT Model

The MMT equation is a one-dimensional nonlinear dispersion equation, which was originally proposed by Majda, McLaughlin and Tabak [62]. The MMT equation gives soliton-like solutions, which have been analyzed in detail by Zhakarov [63,64,65], and gives four-wave resonant interaction between waves, which, when coupled with large-scale forces and small-scale damping, yields a family of solutions that exhibit direct and inverse cascades [22]. The MMT equation is given by:

$$i{\psi }_t=|{\partial }_x{|}^{\alpha }+\lambda |{\partial }_x{|}^{\beta /4}\left(\right||{\partial }_x{|}^{\beta /4}{\psi |}^2|{\partial }_x{|}^{\beta /4}\psi )+iD\psi$$

(28)

where $\psi$ is a complex scalar and $|{\partial }_x{|}^{\alpha }$ is the pseudodifferential operator defined on the real axis through the Fourier transform:

$$|{\partial }_x{|}^{\alpha }\psi \left(k\right)={\left|k\right|}^{\alpha }\widehat{\psi }\left(k\right)$$

(29)

The last term in (28) is the dissipation term, which is tuned to fit ocean waves through the Laplacian operator, D$\psi$, defined in the Fourier space:

$$\widehat{D\psi }\left(k\right)=\left\{\begin{array}{cc}-(|k|-k^{\ast }{)}^2\widehat{\psi }\left(k\right)\hfill & \left|k\right|>k^{\ast }\hfill \\ 0\hfill & \left|k\right|⩽k^{\ast }\hfill \end{array}\right.$$

(30)

This dissipation term, used by [22], is similar to other dissipation models used by Komen and colleagues [66], who have developed concrete models for simulating large wave groups with focusing and defocusing effects. $\lambda$ is the nonlinearity coefficient and corresponds to the focusing phase when <0 and to the defocusing phase when >0. The MMT Equation (28) differs from the standard NLSE in that its family of solutions develops in a more exponential pattern, rather then the Gaussian bell-shaped pattern observed for the solutions for the NLSE [22]. The interesting aspect of this pattern of the spectrum of solutions of the MMT equation is in the mode of formation of the rogue wave, where their energy is transferred from and to the surrounding waves. The solutions are in other words induced by the intermittent formation from the localized rogue event arising out of the regular Gaussian background and collapsing into the surrounding waves. The energy of the rogue wave is transferred to the surroundings and experiences a complete zero-point state, merging completely in the background [22].

The MMT model shows also the formation of quasisolitons, which appear in triple-wave packets, as modeled by Zakharov and Pushkarev [63] and differ from regular solitons in that they radiate the energy backwards towards the preceding amplitudes. This behavior of the solutions may be particularly compatible with the simulation of rogue wave events occurring in regions with strong counter-wind currents, such as in the Agulhas Current [67] or in the regions of the Irish sea [2], which are heavily populated by rogue events, on the passage of the warm waters of the Gulf stream when encountering the frequent low-pressure systems over the Irish sea with counter-wave winds. The quasibreathers or quasisolitons [63] have a root function similar to the Dysthe-type solutions given in (25). Zakharov and Pushkarev [63] approach the solutions in the form:

$$\psi \left(t\right)=e^{i(\mathsf{\Omega }-kV)t}{\varphi }_k$$

(31)

where $\mathsf{\Omega }$ and V are constants ($\mathsf{\Omega }<0$ and $V>0$) and k is the wavenumber, which is an approximate solution to the soliton-like solution for the MMT model. In this approximation, [63] give $\varphi$k the following form:

$${\varphi }_k=\lambda \frac{\int T_{1234}{\varphi }_1^{\ast }{\varphi }_2{\varphi }_3\delta (k+k_1-k_2-k_3)d{k}_{1}d{k}_{2}d{k}_3}{-\mathsf{\Omega }+kV-{\left|k\right|}^{\alpha }},$$

(32)

which represents a form that gives quasi-soliton solutions [68] to the MMT equation. This form of the solutions to the MMT equation radiates energy backwards to the proceeding amplitudes and represents therefore an energy focusing, which is rather dissimilar from the focusing effects modeled by others for rogue patterns (vide supra). It is interesting to note that backward radiation plays also a central role in the dynamics of the quasi-solitons, and not only for their energy accumulation profile. Using the MMT model, Zakharov and Pushkarev [63] developed also a model for collapses of the rogue event, by using self-similar solutions, and modeled the formation of the wave wedge in the appearing and vanishing state, given by a Fourier space distribution of the wave function. Zakharov and Pushkarev [63] have also used the MMT model to develop a turbulence-based solution for the localized rogue event, using the initial condition in the form of an NLSE soliton:

$$\psi (x,0)=\frac{q}{2k_m^{9/4}}\frac{e^{i{k}_{m}x}}{cosh\left(qx\right)},$$

(33)

which shows a conserved action and momentum and an “inner turbulence” localized both in the real and Fourier spaces of the solutions to the modeled envelope function. This “intrinsic turbulence” is described by the authors as affecting the form of its wave spectra, which is irregular and has a stochastic behavior [63]. This model of the rogue wave shows quasi-periodic oscillations with slowly diminishing amplitudes over time, caused by the destruction of rogue wave by the surrounding interference, which the authors denoted as a “quasi-breather”.

6. The Hirota Equation

Multisolitons and breathers for rogue waves have been also successfully modeled [69] by applying the Darboux transformation on the Hirota equation [53]. In their approach, Tao and He [69] developed the Lax pair on the Hirota equation, by using the AKNS [70] procedure to get the Lax pair with the spectral parameters of the Hirota equation given below:

$$i{\psi }_t+{\alpha \left(2\right|\psi |}^2\psi +{\psi }_{xx})+i\beta ({\psi }_{xxx}+6{\left|\psi \right|}^2{\psi }_x=0$$

(34)

where the Lax pair is expressed as:

$${\varphi }_x=M\varphi ,{\varphi }_t=N\varphi$$

(35)

giving rise to the extended matrix representation of the operators in the Hirota equation as described by Tao and He [69]. Tao and He further applied the Darboux transformation [43] on the Lax-represented system by using the simple gauge transformation for spectral problems,

$${\varphi }^{\left[1\right]}=T\varphi$$

(36)

where T is the polynomial applied on the parameter $\lambda$ given in the Lax pair and $\varphi$ is the seed function. Tao and He [69] argue however that regular seed solution $\varphi =e^{ix}$ as described in the previous sections is too special and makes the rogue wave model not universal enough. Tao and He develop therefore a different seed function compared to Akhmediev and colleagues [19,44], for instance, and develop a more extended form of the seed function by starting from a zero seed solution and a periodic seed solution to construct the complete solutions for the breathers and solitons. At zero seed and with the parameter $\lambda$ from the Lax pair, they set the following Hermitian seed pair:

$$\varphi =e^{-i(\xi +i\eta )x-[4\beta i{(\xi +i\eta )}^3+2\alpha i{(\xi +i\eta )]}^{2}t}$$

(37)

and:

$${\varphi }^{\ast }=e^{i(\xi +i\eta )x+[4\beta i{(\xi +i\eta )}^3+2\alpha i{(\xi +i\eta )]}^{2}t}$$

(38)

back in the Darboux transformation to get the one-soliton solution:

$${\psi }_{soliton}^{\left[1\right]}=2\eta e^{2i(-i\xi x-4\beta {\xi }^{3}t-2\alpha {\xi }^{2}t+12\beta \xi {\eta }^{2}t+2\alpha {\eta }^{2}t)}\times sech(-2\eta x-24\beta {\eta }^{2}t+8\beta {\eta }^{3}t-8\alpha \eta \xi t).$$

(39)

Tao and He [69] further report the model for the two-soliton solution and finally give the form of the one-soliton breather solution:

$${\psi }_{breather}^{\left[1\right]}=e^{i\varphi }\left[c-\frac{2\eta [\eta cosh\left(2d_2\right)-i\sigma sinh\left(2d_2\right)-ccos\left(2d_1\right]}{ccosh\left(2d_2\right)-\eta cos\left(2d_1\right)}\right],$$

(40)

where $d_1$, $d_2$, $\sigma$ are given in [69]. Tao and He finally construct the Rogue wave solutions to the original Hirota equation in (34) by Taylor expansion on the breather solutions in (40). The Taylor expansion is carried out on the $\eta$ variable of the breather solution (40), which is given in [69], and forms the general form of the first order rogue wave of the Hirota equation:

$${\psi }_{roguewave}=k{e}^{i(-2\xi x+\beta t)}\left(1-\frac{2k_1+2k_2+i{k}_{3}t}{k_1-k_2}\right)$$

(41)

where the polynomials $k_1$, $k_2$, $k_3$ are given by Tao and He [69]. The rogue wave model resulting from this form is more general than the model given by Akhmediev and colleagues [71] on the Hirota equation. This difference is caused by the appearance of several parameters related to the eigenvalues of the Lax pairs and gives however a possibility to tune more finely the model to experiments on rogue waves. This advantage of the model by Tao and He increases the ability to modulate the precision of reproducing a rogue wave model by calculations. Tao and He’s method grants also the possibility of calculating higher order rogue wave solutions to the Hirota equation by determinant representation of the Darboux transform, which was carried out in a subsequent work [72].

7. The Ablowitz–Musslimani Models: Non-Local Rogue Waves

Another critical method for modeling rogue waves was developed by Ablowitz and Musslimani [73,74,75,76] and uses nonlocal integrable models of the NLSE (1) and KdV (15) equations, where the resulting wave is derived by reverse space-time symmetry. The model evolves by establishing integrability by an infinite number of constants of motion or an infinite number of conservation laws. By this, the method uses a compatible pair of linear equations (similar to the Lax pair in the Hirota equation in (35)) with the nonlinear integrable equation. The method by Ablowitz and Musslimani differs from the Hirota method in that the pair of linear equations represent the scattering problem and the evolution of the scattering data [73,76]. Furthermore, the method by Ablowitz and Musslimani is different from others in that it constructs an inverse scattering problem also known as a linear Riemann–Hilbert problem, which gives the solution to the nonlinear PDE with dependency on time.

The approach by Ablowitz and Musslimani [73] starts by linearizing the equation:

$$i{q}_t(x,t)=q_{xx}(x,t)\pm 2q(x,t)q^{\ast }(-x,t)q(x,t),$$

(42)

where one can immediately observe the existence of a Hermitian pair with reverse directional variables. This form, where reverse variables are used, defines the nonlocal property of the equation and has the advantage that the equation remains invariant in time and space, after the complex conjugate is taken. Hence, the nonlocal equation is parity- and time-symmetric ($\mathcal{PT}$-symmetric), which prevents the equation from yielding different results by a self-induced potential.

An exemplary Lax pair is given in [77] as:

$$v_x=\left(\begin{array}{cc}-ik& q(x,t)\\ r(x,t)& ik\end{array}\right)v,$$

(43)

$$v_t=\left(\begin{array}{cc}A& B\\ C& -A\end{array}\right)v,$$

(44)

where $\mathit{v}$ is the two-component vector and $\mathit{k}$ is a special parameter, while A and B are complex functions. Ablowitz and Musslimani use this step specific compatibility conditions [78] to transform the original PDE in (42) (i.e., ${\psi }_{xt}={\psi }_{tx}$) and gain the simplified PDE pair:

$$i{q}_t(x,t)=q_{xx}(x,t)-2r(x,t)q^2(x,t),$$

(45)

$$-i{r}_t(x,t)=r_{xx}(x,t)-2q(x,t)r^2(x,t).$$

(46)

which yields the original form in Equation (42). Ablowitz and Musslimani further define the nonlocality by using a specific symmetry reduction:

$$r(x,t)=\mp q^{\ast }(-x,t).$$

(47)

This step is particularly characteristic of Ablowitz–Musslimani models [73,74,75,76], from which the new class of nonlocal integrable evolution equations with the nonlocal NLSE hierarchy is directly derived.

The aforementioned property of conserved quantities and conservation laws is also characteristic of Ablowitz–Musslimani models [73,75]. Here, they define a set of eigenfunctions that obey specific boundary conditions [73]. The eigenfunctions are very similar to the seed functions used by other groups and, when inserted in the Lax pair, yield a Riccati model of the conservation quantities. This yields the global conservation laws, which are given in [73] by:

$${\mathcal{C}}_0={\int }_{-\infty }^{+\infty }q(x,t)q^{\ast }(-x,t)dx,$$

$${\mathcal{C}}_1={\int }_{-\infty }^{+\infty }[q_x(x,t)q^{\ast }(-x,t)+q(x,t)q_x^{\ast }(-x,t)]dx,$$

$${\mathcal{C}}_2={\int }_{-\infty }^{+\infty }[q_x(x,t)q_x^{\ast }(-x,t)+\sigma q^2(x,t)q^{\ast 2}(-x,t)]dx,$$

which are real integrable Hamiltonians. Ablowitz and Musslimani [73] derive furthermore local conservation laws defined by the equations:

$${\partial }_t\left[q(x,t)q^{\ast }(-x,t)\right]+i{\partial }_x[q(x,t)q_x^{\ast }(-x,t)+q^{\ast }(-x,t)q_x(x,t)]=0$$

(48)

$${\partial }_t\left[q(x,t)q_x^{\ast }(-x,t)\right]+i{\partial }_x[q_x^{\ast }(-x,t)q_x(x,t)+q(x,t)q_{xx}^{\ast }(-x,t)-\sigma q^2(x,t)q^{\ast 2}(-x,t)]=0$$

(49)

which are used to develop the framework for the direct scattering problem and the inverse scattering problem, where the scattering data are given by specific scattering matrices. The same symmetry is also in the problem of the potential and in the eigenfunction and leads naturally to the same symmetry relation in the scattering matrices, which is given by:

$$N(x,k)=\mathsf{\Lambda }{M}^{\ast }(-x,-k^{\ast })$$

(50)

and:

$$\overline{N}(x,k)={\mathsf{\Lambda }}^{-1}\overline{M}\ast (-x,-k^{\ast }),$$

(51)

where $\mathsf{\Lambda }$ is a 2 × 2 matrix with zeros in the diagonal and one, $\pm 1$ on the lower and upper diagonal, respectively. For the inverse scattering problem, Ablowitz and Musslimani [73] account for the symmetry condition by considering the set of basis terms as a left scattering problem and supplement these terms with the equivalent right-scattering problem, from which they formulate the Riemann–Hilbert problem and find the linear integral equations that govern the functions M and $\overline{M}$ in (50). These equations are given by:

$$M(x,k)=\left(\genfrac{}{}{0pt}{}{1}{0}\right)+\sum _{l=1}^{\overline{J}}\frac{{\overline{B}}_l{e}^{2i{\overline{k}}_{l}x}\overline{M}(x,{\overline{k}}_l)}{k-{\overline{k}}_l}-\frac{1}{2\pi i}{\int }_{-\infty }^{\infty }\frac{\overline{R}\left(\xi \right)e^{2i\xi x}\overline{M}(x,\zeta )}{\zeta -(k+i0)}d\zeta ,$$

(52)

$$\overline{M}(x,k)=\left(\genfrac{}{}{0pt}{}{0}{1}\right)+\sum _{l=1}^J\frac{B_l{e}^{2i{k}_{l}x}M(x,k_l)}{k-k_l}-\frac{1}{2\pi i}{\int }_{-\infty }^{\infty }\frac{R\left(\xi \right)e^{2i\xi x}M(x,\zeta )}{\zeta -(k+i0)}d\zeta ,$$

(53)

where $\mathit{R}\left(\mathit{k}\right)$ and $\overline{R}\left(x\right)$ are the reflection coefficients. The terms $B_l$ and ${\overline{B}}_l$ are the conservation law Hamiltonians applied symmetrically. From this stage, Ablowitz and Musslimani [73] derive a linear algebraic integral system of equations that solves the inverse problem for the eigenfunctions $\overline{M}(x,k)$ and $M(x,k)$.

The resulting soliton solutions of the Ablowitz–Musslimani model hence assume the form:

$$q\left(x\right)=-\frac{2(\eta +\overline{\eta })e^{i{\overline{\theta }}_1}{e}^{-4i{\overline{\eta }}_1^{2}t}{e}^{-2{\overline{\eta }}_{1}x}}{1+e^{i({\theta }_1+{\overline{\theta }}_1)}{e}^{4i({\eta }_1^2-{\overline{\eta }}_1^2)t}{e}^{-2({\eta }_1+{\overline{\eta }}_1)x}},$$

(54)

which represents a family of solutions defined by the four independent parameters, which have a dynamic relationship with the time variable and which gradually develop a singularity in a finite time period, $t_s$ at x = 0, where $t_s$ is given by:

$$t_s=\frac{(2n+1)\pi -{\theta }_1-{\overline{\theta }}_1}{4({\eta }_1^2-{\overline{\eta }}_1^2)}.$$

(55)

Equation (55) is a critical form of the time-variable, which distinguishes the method of Ablowitz–Musslimani [73,74,75,76,77,78] from other rogue wave models and adds a non-linear evolution of the rogue wave. Ablowitz and Musslimani have also most recently developed a new model that includes nonlocal rogue waves with nonzero background, which provide a more realistic view of the rogue wave, which focuses energy from neighboring waves [79].

Solutions to the Ablowitz–Musslimani model (42) were also developed by Yang and Yang [80], who used the Darboux transformation method on the PDE coupled with the Bäcklund transformation on the potential functions, identifying three types of rogue waves from the Ablowitz–Musslimani picture. Yang and Yang expanded the solutions to polynomials using Schur polynomials [80]. This analysis of the Ablowitz–Musslimani model showed greater variation in the rogue waves compared to the regular NLSE, where the variations were represented by the terms in the denominator of the soliton solutions.

The parity-time symmetry potential of the Ablowitz–Musslimani equations has also been studied by Yu [81] very recently, who obtained discrete rogue wave solutions with three free parameters (refer to Equation (54) for similarities). Yu studied in particular the effect that the dispersion of the parity-time symmetry has on the solutions, as well as the effect of the coefficients and the parameters. Yu [81] used the Darboux transformation method in a similar fashion to Yang and Yang [80] to derive different forms of solitons with different heights, which were defined by two of the three free parameters in the solution ($\eta$ and $\overline{\eta }$ in Equation (54)). Yu [81] also assessed the stability of rogue waves over a specific period of time and included a modulation instability coefficient that allows the modeling of several discrete solutions that represent various stages of rogue wave formation (appearing suddenly and disappearing suddenly), a property of rogue waves also reported by Akhmediev and colleagues [18,19,20]. Finally, Yu modeled rogue waves that appear rapidly and do not disappear. This latter model may be particularly relevant to describe rogue events during low-pressure systems at open sea, which have been reported in several cases to give stable rogue waves with a long lifetime (i.e., the rogue waves reported in the study by Munk [35]).

8. Conclusions

A survey of various mathematical models for representing rogue waves has herein been carried out to the maximum extent of including the most applied, as well as most recent models. The overall survey yields a perspective that delineates the common traits between the methods. This survey also shows how novel and emerging models allow for better modeling of rogue waves, by including several parameters associated with the evolution of the rogue waves, such as the duration, the height and other particular properties that a rogue event can display when occurring in oceans, optical systems or even in the atmosphere. The results show also that the new versatile forms of the older models (MMT and Dysthe) can be furthermore adapted and studied numerically in upcoming papers, particularly with relevance to counter-wave winds, which frequently occur on the Gulf Current off the Irish coast and on the western Norwegian coast during northeast and north winds, respectively. By the overall review of the various methods, a theorem is suggested for describing the origin of rogue waves in the ocean. The theorem suggests that a rogue wave in the ocean can be formed whenever there is a momentaneous surplus of energy perturbed on the momentum or in the kinetic term of a wave train, induced either by a sudden change in the atmosphere leading to strong winds appearing suddenly over large volumes of water, or induced by a collision of large volumes of water with highly different temperatures and densities, or finally, as often observed, a rogue event occurs by the constructive overlap of waves, in opposite directions, in traverse directions or running in the same direction, and its duration is determined, when occurring in the same direction, by the slight deviations in the momenta of the overlapping waves. Future work will be the submission of a project proposal for predicting rogue waves for off-shore structures, currently under development for the Norwegian Research Council.