The phenomenon of rogue waves (RWs)—unusually large waves that appear suddenly from moderate wave background–was intensively studied in the recent years. A number of mechanisms were suggested to explain their emergence, see e.g., the reviews [1–3]; with the most general idea stating that RWs could be related to breather-type solutions of the underlying nonlinear evolution equations [4–6]. Currently, ones of the most popular models for RWs are the Peregrine rational breather [7] and the higher-order rational breather [8] solutions of the one-dimensional nonlinear Schrödinger equation (1D-NLSE) of the focusing type, i ψ t + ψ x x + | ψ | 2 ψ = 0 (1)

These solutions belong to a family of localized in space and time algebraic breathers, which evolve on a finite background and lead to three-fold, five-fold, seven-fold, and so on, increase in amplitude at the time of their maximum elevation. Taking specific and carefully designed initial conditions, they were reproduced in well-controlled experiments performed in different physical systems [9–13].

The 1D-NLSE is integrable in terms of the inverse scattering transform (IST), as it allows transformation to the so-called scattering data, which is in one-to-one correspondence with the wavefield and, similarly to the Fourier harmonics in the linear wave theory, changes trivially during the motion. Thanks to its properties, the scattering data can be used to characterize the wavefield. For spatially localized case, the scattering data consists of the discrete (solitons) and the continuous (nonlinear dispersive waves) parts of eigenvalue spectrum, calculated for specific auxiliary linear system. For strongly nonlinear wavefields, such as the ones where emergence of rational breathers can be expected, the solitons provide the main contribution to the energy [14] and should therefore play the dominant role in the dynamics. In particular, as has been recently demonstrated in [15]; the modulationally unstable plane wave (the condensate) at its long-time statistically stationary state can be accurately modeled (in the statistical sense) with a certain soliton gas, designed to follow the solitonic structure of the condensate. The latter naturally raises a question of whether there is a difference between the RWs emerging in the two systems. Indeed, in a soliton gas all RWs are multi-soliton interactions by construction. Hence, if there is no significant difference, then we can draw a hypothesis that for the asymptotic stationary state of the MI (and, possibly, for other strongly nonlinear wavefields) the main mechanism of RW formation is interaction of solitons.

With the present paper, we contribute to the answer on this question by summarizing our observations of RWs for both systems. Specifically, we compute time evolution for 1,000 random realizations of the noise-induced MI of the condensate and also for 1,000 random realizations of 128-soliton solutions modeling the asymptotic state of the MI. For each realization, we analyze one largest RW emerging in the course of the evolution, thus focusing our analysis on the largest RWs. For both systems, we observe practically identical dynamical and statistical properties of the collected RWs. In particular, most of the RWs turn out to be very well approximated–simultaneously in space and in time–by the amplitude-scaled rational breather solution (RBS) of the second order. By measuring deviation between RWs and their fits with RBS as an integral of the difference in the ( x , t ) -space, we find that, in general, the larger the maximum amplitude of the RW, the better its convergence to the RBS of the second order (RBS2). The collected RWs for the two systems turn out to be identically distributed by their maximum amplitude and deviation from the RBS2. Additionally, we demonstrate that the observed quasi-rational profiles appear already for synchronized three-soliton interactions and discuss the next steps in the ongoing research of the RW origin.

Note that in the present paper we consider solutions of the 1D-NLSE for three different types of boundary conditions: the MI of the condensate for which we use the periodic boundary, the multi-soliton solutions with vanishing border conditions and the RBS having constant border conditions at infinity. Globally, these solutions are fundamentally different, and the different border conditions require application of separate IST techniques, see e.g., [5, 14, 16, 17]. For instance, formally our MI case corresponds to finite-band scattering data. However, the characteristic widths of the structures (RWs, solitons, RBS) are small compared to the sizes of the studied wavefields, so that the eigenvalue bands are very narrow and we neglect their difference from solitons. The similar idea was suggested in [18]; where, vice versa, the soliton gas was considered as a limit of finite-band solutions. Effectively, we assume that formation of a RW, as a local phenomenon, represents a similar process for all three cases of border conditions. As we demonstrate in the paper, this assumption is supported by the presented results, that raises an important problem that we leave for future studies–explanation of how the three models may exhibit locally similar nonlinear patterns.

For both the MI of the condensate and the soliton gas initial conditions, we solve Eq. 1 numerically in a large box x ∈ [ − L / 2 , L / 2 ] , L ≫ 1 , with periodic boundary. We use the pseudo-spectral Runge-Kutta fourth-order method in adaptive grid with the grid size Δ x set from the analysis of the Fourier spectrum of the solution; see [19] for detail. As an integrable equation, the 1D-NLSE conserves an infinite set of integrals of motion, see e.g., [14]. We have checked that the first ten integrals are conserved by our numerical scheme up to the relative errors from 10 − 10 (the first three invariants) to 10 − 6 (the 10th invariant) orders.

Without loss of generality, the initial conditions for the noise-induced MI of the condensate can be written as ψ | t = 0 = 1 + ε ( x ) (2) where ε ( x ) represents a small initial noise. We use statistically homogeneous in space noise with Gaussian Fourier spectrum, ε ( x ) = a 0 ( 8 π θ L ) 1 / 2 ∑ k e − k 2 / θ 2 + i ϕ k + i k x (3) where a 0 is the average noise amplitude in the x-space, k = 2 π m / L is the wavenumber, m ∈ ℤ is integer, θ is the characteristic noise width in the k-space and ϕ k are random phases for each k and each realization of the initial conditions; the average intensity of such noise equals to a 0 2 , 〈 | ε | 2 〉 = a 0 2 . For the numerical experiment, we take the box of length L = 256 π and small initial noise, a 0 = 10 − 5 , with wide spectrum, θ = 5 . Note that these parameters match those used in [19].

To generate soliton gas, modeling the asymptotic stationary state of the noise-induced MI, we use exact 128-soliton solutions of the 1D-NLSE. More precisely, we compute the corresponding wavefields using numerical implementation of the 1D-NLSE integration technique–the so-called dressing method [20]—with 100-digits precision arithmetics, as described in [21]. Each soliton has four parameters: amplitude a j , velocity v j , space position x 0 j and phase Θ j ; here j = 1 , … , M , M = 128 , and the one-soliton solution reads as ψ s ( x , t ) = a exp [ i υ 2 ( x − x 0 ) + i 2 ( a 2 − υ 2 2 ) t + i Θ ] cosh a ( x − x 0 ) − a υ t 2

Following [15]; we distribute soliton amplitudes according to the Bohr-Sommerfeld quantization rule, a j = 2 1 − ( j − 1 / 2 M ) 2 (4) and set soliton velocities to zero, υ j = 0 , using uniformly-distributed soliton phases Θ j in the interval ( 0,2 π ) and uniformly-distributed space position parameters x 0 j in a narrow interval at the center of the computational box. Zero velocities mean that these multi-soliton solutions are bound-state. For the 1D-NLSE in normalization Eq. 1, the Bohr-Sommerfeld rule describes amplitudes for the bound-state solitonic content of a rectangular box wavefield of unit amplitude ψ = 1 and width L o = 2 π M , calculated with the semi-classical Zakharov-Shabat direct scattering problem, see e.g., [14, 22, 23]. The generated 128-soliton solutions take values of unity order approximately within the interval x ∈ [ − L o / 2 , L o / 2 ] and remain small outside of it. For more detail on the soliton gas, we refer the reader to [15]; where it has been demonstrated that its spectral (Fourier) and statistical properties match those of the long-time statistically stationary state of the MI.

For the soliton gas, we gather the RWs by simulating time evolution of the 128-soliton solutions in the interval t ∈ [ 0,50 ] and then collecting one largest RW for each of the 1,000 realizations of initial conditions. For time evolution, we use the same pseudo-spectral Runge-Kutta numerical scheme as for the MI of the condensate, since application of the dressing method with evolving scattering data takes too much computational time and provides the same result. The pseudo-spectral scheme uses periodic boundary conditions, so that solution ψ ( x , t ) needs to be small near the edges of the computational box. We achieve this by taking the box of length L = 384 2 π , so that our 128-soliton solutions are of 10 − 16 order near its edges and take values of unity order, | ψ | ∼ 1 , only within its central 1 / 3 part ( ≡ L o / L ) .

For the MI of the condensate, we collect the RWs similarly, but in the time interval t ∈ [ 174,200 ] . From the one hand, the end of this interval is far enough, so that the system is sufficiently close to its asymptotic stationary state, see [19] where the same initial conditions were used. From the other hand, a chance to detect a large RW is higher in larger simulation boxes and if we wait longer. To make RW events for the two systems comparable, we impose a restriction L ( M I ) ⋅ Δ T ( M I ) = L ( S G ) ⋅ Δ T ( S G ) on the lengths L ( M I , S G ) of the regions where RWs may appear and on the time intervals Δ T ( M I , S G ) during which we wait for the largest RW. For the soliton gas, the collected RWs appear approximately in the space interval x ∈ [ − 210,210 ] with practically uniform distribution of their position, so that L ( S G ) = 420 . We believe that this property is connected with behavior of the ensemble- and time-averaged intensity I ( x ) = 〈 | ψ ( x , t ) | 2 〉 , which remains flat I = 1 inside this interval and starts to deviate from unity at its edges. For the MI, the RWs may appear anywhere within the computational box L ( M I ) = 256 π ; together with the observation time for the soliton gas case Δ T ( S G ) = 50 , this yields Δ T ( M I ) = 26 and the time interval t ∈ [ 174,200 ] for the MI.

The rational breather solution of the first order (RBS1)—the Peregrine breather [7]—reads as ψ p ( 1 ) ( x , t ) = e i t [ 1 − 4 ( 1 + 2 i t ) 1 + 2 x 2 + 4 t 2 ] (5)

The RBS of the second order (RBS2) ψ p ( 2 ) is too complex and we refer the reader to [8] where it was first found. Both solutions describe localized in space and in time perturbations that evolve on a finite background–the condensate–and lead to three-fold and five-fold increase in the overall amplitude at the time of their maximum elevation. For approximation of a RW with a RBS, we use the scaling, translation and gauge symmetries of the 1D-NLSE: indeed, if u ( x , t ) is a solution of Eq. 1, then A 0 e i Θ ⋅ u ( χ , τ ) , where χ = A 0 ( x − x 0 ) , τ = A 0 2 ( t − t 0 ) , and A 0 , Θ ∈ ℝ , is also a solution. Technically, we detect the maximum amplitude A of a RW together with its position x 0 and time t 0 of occurrence, and also the phase at the maximum amplitude Θ = arg ψ ( x 0 , t 0 ) , and then use the scaling coefficient A 0 = − A / 3 for the RBS1 and A 0 = A / 5 for the RBS2.

Note that, in general, RBS may have nonzero velocity v ≠ 0 . To account its influence, one can make a transformation u ( x , t ) → e i v x / 2 − i v 2 t / 4 ⋅ u ( x − v t , t ) , which also prompts a simple way to find the velocity. Indeed, at the time of the maximum elevation t 0 , a RBS with zero velocity, v = 0 , has constant phase arg ψ ( x , t 0 ) = c o n s t in the region between the two zeros closest to the maximum amplitude. In contrast, a RBS with nonzero velocity, v ≠ 0 , has constant phase slope, arg ψ ( x , t 0 ) − i v x / 2 = c o n s t , in the same region. Hence, by computing the phase slope one can approximate RWs with RBS of nonzero velocity. For all the RWs studied in this paper, we have checked that taking into account velocity improves our approximations only very slightly, and for this reason we have decided to use RBS with zero velocity only.

Also note that in addition to the RBS1 and the RBS2, we have examined approximation with the RBS of the 3rd order [8]; as well. However, we have found that it works worse than either the RBS1, or the RBS2 for all 2000 examined RWs, and thereby excluded it from the analysis.

As we have mentioned in [21]; some two- and three-soliton collisions at the time of their maximum elevation have space profiles remarkably similar to those of the RBS1 and the RBS2. Moreover, we have presented an example of a phase-synchronized three-soliton collision, for which both the space profile and the temporal evolution of the maximum amplitude are very well approximated by the RBS2. The solitons in [21] had nonzero velocities; here we modify the two- and three-soliton examples for the case of zero velocities and also examine the local deviations d p ( 1,2 ) ( x , t ) Eq. 6 together with the integral deviations D p ( 1,2 ) Eq. 7.

Figure 3 shows an example of three-soliton interaction with solitons having amplitudes a 1 = 1 , a 2 = 1.5 and a 3 = 2 , zero velocities v j = 0 , zero space position parameters x 0 j = 0 and, at the initial time t = 0 , zero phases Θ j = 0 . The space profile | ψ ( x , t 0 ) | at the time of the maximum elevation t 0 = 0 is remarkably similar to that of the RBS2, and the local deviation d p ( 2 ) ( x , t ) remains well within 5 % for most of the area presented in the figure as well. The integral deviation Eq. 7 equals to D p ( 2 ) ≈ 0.016 , that is even smaller than for the RW presented in Figure 1.

To analyze how often phase-synchronized interactions of two and three solitons of various amplitudes may lead to such quasi-rational profiles, we have created 20 two-soliton and 20 three-soliton interactions with solitons of random amplitudes, zero velocities v j = 0 , zero space positions parameters x 0 j = 0 and phases Θ j = 0 . For the two-soliton interactions, the minimum deviations from the RBS1 and the RBS2 turned out to be D p ( 1 ) ≈ 0.077 and D p ( 2 ) ≈ 0.061 , and the average ones— 〈 D p ( 1 ) 〉 ≈ 0.14 and 〈 D p ( 2 ) 〉 ≈ 0.075 , respectively. For the three-soliton case, the minimum deviations were D p ( 1 ) ≈ 0.18 and D p ( 2 ) ≈ 0.003 , and the average ones— 〈 D p ( 1 ) 〉 ≈ 0.23 and 〈 D p ( 2 ) 〉 ≈ 0.022 ; the maximum deviation from the RBS2 equaled to D p ( 2 ) ≈ 0.03 , that is still very good for comparison with the RBS2. Hence, we conclude that quasi-rational profiles very similar to that of the RBS2 appear already for three-soliton interactions, provided that the solitons are properly synchronized (that is, have coinciding positions and phases).

We think that the presented elementary three-soliton model might provide an explanation of RW formation inside multi-soliton solutions. The most direct way for future studies might be a demonstration of a RW for synchronized many-soliton solution. Here, however, we face a new question, that is, whether formation of a RW is a collective phenomenon that requires synchronization of all the solitons, or a “local” event that can be achieved by synchronizing of a few with arbitrary parameters of the others. Note that even the latter case represents a challenging problem. Indeed, the few solitons that generate a RW acquire space and phase shifts due to the presence of the remaining solitons, that should influence their optimal synchronization condition. For remote solitons, the shifts can be computed analytically using the well-known asymptotic formulas, see e.g., [14]; which however do not work for our case of a dense soliton gas where all solitons effectively interact with each other. This leaves us two options: 1) a local numerical synchronization of a small group with “trial and error” method and 2) calculation of generalized space-phase shifts for closely located solitons. Both ways seem challenging at the moment.

Also note that our study is limited with respect to statistical analysis of RWs, as we have focused on the largest RWs, while the “common” RWs may have different dynamical and statistical properties. Nevertheless, we believe that, since the largest RWs for the two systems show identical properties, the “common” RWs have identical properties too. Identification of all the RWs according to the standard criterion A ≥ 2.8 is a nontrivial problem by itself, and we plan to return to it in the near future.

In this brief report we have presented our observations of RWs within the 1D-NLSE model for 1) the modulationally unstable plane wave at its long-time statistically stationary state and 2) the bound-state multi-soliton solutions representing the solitonic model of this state. Focusing our analysis on the largest RWs, we have found their practically identical dynamical and statistical properties for both systems. In particular, most of the RWs turn out to be very well approximated–simultaneously in space and in time–by the amplitude-scaled rational breather solution of the second order (RBS2), see Figures 1, 2 and the two sets of the collected RWs are identically distributed by their maximum amplitude and deviation from the RBS2, see Figures 2C,D. Additionally, we have demonstrated the appearance of quasi-rational profiles very similar to that of the RBS2 already for synchronized three-soliton interactions, see e.g., Figure 3.

The main messages of the present paper can be summarized as follows. First, a quasi-rational profile very similar to a RBS does not necessarily mean emergence of the corresponding rational breather, as it can be a manifestation of a multi-soliton interaction. Second, the identical dynamical and statistical properties of RWs collected for the two examined systems strongly suggest that the main mechanism of RW formation should be the same, i.e., that RWs emerging in the asymptotic stationary state of the MI (and, possibly, in other strongly nonlinear wavefields) are formed as interaction of solitons. However, more study is necessary to clarify exactly how interaction of solitons within a large wavefield may lead to formation of a RW, and we plan to continue this research in the future.

As future directions of our studies, we consider two problems. The first is about examination of whether formation of a RW in a soliton gas is a collective phenomenon that requires synchronization of all the solitons, or a “local” event that can be achieved by synchronizing of a few. The second problem relates to statistical analysis of all the RWs according to the standard criterion A > 2.8 , in contrast to the largest RWs analyzed in this paper.