Rogue solitons in optical fibers: a dynamical process in a complex energy landscape?

Nondeterministic giant waves, denoted as rogue, killer, monster, or freak waves, have been reported in many different branches of physics. Their physical interpretation is however still debated: despite massive numerical and experimental evidence, a solid explanation for their spontaneous formation has not been identified yet. Here we propose that rogue waves [more precisely, rogue solitons (RSs)] in optical fibers may actually result from a complex dynamical process very similar to well-known mechanisms such as glass transitions and protein folding. We describe how the interaction among optical solitons produces an energy landscape in a highly dimensional parameter space with multiple quasi-equilibrium points. These configurations have the same statistical distribution of the observed rogue events and are explored during the light dynamics due to soliton collisions, with inelastic mechanisms enhancing the process. Slightly different initial conditions lead to very different dynamics in this complex geometry; a RS turns out to stem from one particularly deep quasi-equilibrium point of the energy landscape in which the system may be transiently trapped during evolution. This explanation will prove to be fruitful to the vast community interested in freak waves.


INTRODUCTION
Observations of nondeterministic giant wave events have been reported in many different fields, from oceanic dynamics [1,2], to financial markets [3], and light propagation in optical fibers [4][5][6][7]. The origin of these phenomena, denoted as rogue, killer, monster, or freak waves, is still unknown. Since the first report of optical rogue solitons (RSs) in supercontinuum generation (SCG) occurring in optical fibers [4], a great effort has been put in characterizing these light pulses with extremely large amplitudes. Despite the massive numerical and experimental evidence [5,6,[8][9][10][11][12][13][14], a unified approach that describes their spontaneous formation has yet to be identified. Although it was suggested that RSs are exact solutions of the nonlinear Schrödinger (NLS) equation, such as the Akhmediev's breather [9,[15][16][17][18] or the Peregrine soliton [19], or rather solutions of a more general model [20][21][22], there are strong indications that they originate from a combination of higher-order effects such as Raman scattering, selfsteepening, and higher-order dispersions [8,23]. Actually, already in 1989, the appearance of transient pulses that strongly remember rogue wave solitons was identified in Raman systems [24] and, more recently, extreme events (in the form of RSs and rains of solitons) were found in dissipative systems [25][26][27][28][29] and even in integrated optical resonators on chips [30]. The exact role of those effects is still debated, because approaches based on integrable turbulence [17] or the Zakharov-Shabat problem [18] seem to provide insightful results, but a link between the regular NLS dynamics and the observed perturbed evolution is still unclear [31]. Exact solutions of integrable models are pivotal to explain the nonlinear stages of modulation instability (MI), which is observed under continuous wave excitation, but their properties (amplitude, finite background, periodic nature) are in contrast to the arbitrary strong pulses that are detected when exciting the optical fiber with ultrashort pulses. More importantly, this class of solutions does not give any information about the statistical distribution of high-intensity peaks, which, like many extreme natural and social phenomena, shows a heavy-tail behavior [4,11,32].
At variance with previous studies, we treat the nonlinear light propagation in an optical fiber as a dynamical system out of equilibrium. In a complex medium such as a molecular glass or softcolloidal matter, the local interactions among adjacent particles form a global potential, the energy landscape [33][34][35][36][37]. The energy landscape exhibits many local minima and proves a useful concept in different branches of physics to describe, for example, the supercooling phase transition from liquid to glass. The literature on kinetic equations for a soliton gas has a long history [38][39][40] and the soliton-bound states emerging from noise or random external forces were also studied extensively [41,42]. A similar statistical mechanics approach is found in the theory of mode locking [43][44][45] and turbulence in atmospheric propagation [46,47]. We prove instead that our system explores the energy landscape generated by the weak interactions between neighboring NLS solitons, while the dynamics is accelerated by the Raman effect, which leads to the continuous growth of an equivalent temperature. Thus, a RS turns out to have its origin in one of the quasi-equilibrium points of the energy landscape, which corresponds to a deep minimum in the interaction energy. This minimum also corresponds to collisions between high-energy pulses emerging during the SCG process. The role of collisions in RS formation was studied numerically elsewhere [10] and collision was shown to play a crucial role and to be an insightful mechanical analog in SCG from ultrashort pulses [48], where the effect of higher-order dispersion is similar to what we find for Raman scattering.
Finally, these quasi-equilibrium points in the landscape exhibit the same heavy-tail statistical distribution of the peaks detected in SCG experiments. This explanation will prove to be extremely useful to all communities dealing with the study of freak waves [49][50][51][52][53][54][55].

A. Energy Landscape
The energy landscape can be built by starting from the simplest model of soliton propagation in fiber optics: the dimensionless NLS equation in anomalous dispersion [56]: where uz; t is the normalized electric field envelope of the pulse, z is the longitudinal coordinate along the fiber, and t is time. We study a regime, commonly observed in experiments and numerical simulations, in which the highly nonlinear dynamics can be described as a system of N weakly interacting solitons. The conventional strategy consists of applying a variational approach and assuming for the solution of Eq. (1) the following Ansatz: where ν k , μ k , ξ k , and δ k represent, respectively, the amplitude, speed, temporal delay, and phase of each soliton, labeled by an index k 1; …; N . It is well established that one can model the nearest-neighbor interactions of NLS solitons as a system of ordinary differential equations for the 4N parameters in Eq. (3) [56][57][58][59][60][61][62][63]. We denote this model (detailed in Supplement 1) as the weak-interaction model (WIM). This model can be extended [64][65][66] to take into account a wide variety of nonintegrable perturbations; however, we show below that the basic WIM captures most of the fundamental features of the system dynamics, and a more sophisticated model will be developed and treated in a more technical work. The WIM can be reduced to a mechanical system under the hypothesis of nearly equal amplitudes ν k ≈ ν. This system is expressed in terms of complex canonical variables. It is outside the scope of the present work to apply the most powerful statistical mechanical machinery of, e.g., [67,68] to this complex Hamiltonian. Instead, we follow an approach similar to Ref. [37] where the appearance and stabilization of filaments in a nonlocal medium were explained as a complex phase transition of light, and we implement a numerical search for the fixed points (stable  A pictorial representation of the energy landscape of a complex system represented as a two-dimensional surface: we observe many local minima with various depths and locations, separated by saddle points. The dynamics is strongly influenced by the initial conditions: we highlight two different trajectories, one dropping into a deep minimum of the potential energy (white solid line) and one experiencing many shallow minima (yellow solid line). This explains why slightly different noise configurations around the input pulse generate completely different dynamics, leading only rarely to RSs. In panel (c), we compare this simplified image with the interaction of many solitons. It is characterized by evolution in a multidimensional complex topology with many minima and saddles. Rare events, due to collisions or higher-order effects (such as dispersion and Raman gain), trigger the system in regions o lower energy, which is in our case the Hamiltonian of the NLS. We represent the energy landscape generated by soliton interaction by different values of H at different numerically computed equilibria. The horizontal axis enumerates the solutions obtained from different initial random conditions. The solution marked by the filled circle is the exceptional solution, which we show in Fig. 2.
or unstable equilibria) of the WIM, i.e., those points for which _ ν k _ μ k _ ξ k _ δ k 0. The basic idea here is to map the highly nonlinear dynamics of rogue wave generation to a complex dynamical topography in a highly dimensional system, as it is performed in many other branches of physics, like supercooled liquids, chemical reactions, and protein folding [33][34][35]. In those systems, it is convenient to represent the interactions as a potential energy landscape, as sketched in Figs. 1(a) and 1(b), for a two-dimensional system. This landscape exhibits many minima of different depths, and saddle points, along which the system can escape and jump randomly from one configuration to another. In a multidimensional setting, the central parameter is the saddle order, which quantifies the ability of the system to escape from a specific minimum configuration.
The WIM corresponds to a 4N -dimensional phase space and we now discuss the properties of its equilibrium points in terms of saddle order and their relation to the minima of an equivalent global energy of the system. We discover that the statistical distributions of these points share many features with the rare RSs observed in the experiments.
The numerical procedure adopted here, and detailed in Supplement 1, allows us to find, for a fixed N , a large variety of fixed points, among which we identify exceptional solutions, which are characterized by a large amplitude soliton ν k , and that correspond to the RS-generating collisions. Moreover, we find many other solutions, which comprise in general smaller multiple unevenly spaced intensity peaks.
In Figs. 2(a) and 2(b), we show an example of an exceptional solution. We observe a typical pattern: a large-amplitude soliton Fig. 2(a1) is closely surrounded by two lower-amplitude solitons [see Fig. 2(a3)] and many others with smaller peak amplitudes. The velocity of the main soliton is much larger than those of its neighbors [ Fig. 2(a2)]; phases are ordered across the triplet. We thus observe a solution composed of a large peak surrounded by smaller oscillations. This is confirmed by the plot of the squared modulus and phase [ Fig. 2 The RS solution appears as a dominant peak superimposed by a smaller one with a different speed. This gives rise to oscillations and phase jumps across the profile. Any exceptional solution is actually composed of many solitons, but these are arranged in a nonperiodic, disordered fashion (a glass  squared amplitude and phase of the reconstructed field (the ellipses and arrows are meant to guide to the corresponding axis), which shows a three-peak structure and an apparent sign change across the pulse (the phase exhibits two subsequent jumps); in the inset, we compare the structure factor of the present solution (blue solid line) with a crystal, the lattice constant of which is Δξ ave (green dashed line). In contrast to (b), we show in (c) two different equilibrium solutions of the WIM, which are much less intense and present a multiple-peak structure: in (c1), the peak is one order of magnitude smaller than that in panel (b) and a trailing pulse appears, whereas (c2) exhibits two peaks. These solutions correspond to the shallower minima in Fig. 1(c). Finally, in (d), the statistical analysis of the ensemble of WIM solutions is presented: a histogram of amplitudes exceeding the Q 0.95 quantile of the distribution (normalized such that the total area of the bars sums up to unity). In the inset, the Weibull fit of the distribution is shown in double logarithmic scale (N 40).

Research Article
Vol. 2, No. 5 / May 2015 / Optica of solitons). To this aim, we define the structure factor Sq F P N n1 δξ − ξ n q, where F · denotes the Fourier transform. We consider the single exceptional realization, i.e., our special solution, and compare it to a perfect crystal of average lattice constant Δξ ave ≡ 1 We show in the inset in Fig. 2(b) the structure factor of our specific exceptional solution, and we notice that in contrast to regularly spaced equal peaks of an ordered arrangement, the plot exhibits irregularly spaced peaks, typical of glassy systems.
The WIM solutions of the kind shown in Fig. 2(b) seem to qualitatively resemble the Peregrine soliton [19] or a slice of an Akhmediev's breather [9,15]. This is, however, misleading, because in our case we do not have an infinite background, but a sum of pulses of different amplitudes, which decay exponentially. The ratio of amplitudes is not fixed and solutions with two peaks are also found. Instead, the Peregrine soliton and similar solutions of higher-order models [21] exhibit a fixed ratio of the peak versus the background. In order to further support the energy landscape hypothesis, we report in Figs. 2(c1) and 2(c2) two different equilibrium solutions: the peak intensity is one order of magnitude smaller than that in Fig. 2(b) and the field profile can have a single or multiple peaks.
We remark two fundamental facts: first, the solutions considered are unstable fixed points of the WIM, with saddle order (i.e., the number of eigenvalues with a positive real part; see Ref. [37]) N s ≈ 2N . Higher-order effects in optical fibers permit exploring a broad range of irregular solutions and the saddle order hints at the probability of falling into such a solution. Second, the exceptional solution corresponds to a minimum value (among the numerically found solutions) of the integrated Hamiltonian density of the NLS, with ℋz; t ≡ ℋ K ℋ NL , ℋ K ≡ ju t j 2 ∕2, and ℋ NL ≡ −juj 4 ∕2.
According to this argument, in Fig. 1(c) we show a representation of the energy landscape in the form of an enumeration of the values of H obtained numerically from each different initial guess. Each point in Fig. 1(c) is a quasi-equilibrium configuration. We notice that most solutions correspond to small or vanishing quantities, many solutions with H > −10 3 are found, and the exceptional solution of Fig. 2(b) clearly emerges. Shallower minima correspond to solutions such as those shown in Fig. 2(c), with smaller peak amplitudes and a multipeak structure.
A posteriori, we should notice that the hypotheses on which the WIM is based, i.e., nearly homogeneous amplitudes and small velocity discrepancies among solitons, are partly violated by our exceptional solution. However, as that class of solutions has a large saddle order, they are rapidly approached and departed from during evolution and can be thought as an ideal extreme state, not exactly achieved by the system.
The statistical analysis of the ensemble of solutions obtained by the numerical solver is shown in Fig. 2(d). We select values of ν exceeding the Q 0.95 quantile and show that it exhibits a heavy-tail behavior. To further confirm this, we fit the data to a Weibull distribution: this probability density function belongs to the class of heavy-tail distributions that occur in the description of many extreme phenomena in natural and social sciences. These include the Pareto, log-normal, and Lévy distributions, and, in contrast to Gaussian distribution, they exhibit at least one tail that is not bound exponentially. The Weibull function was used earlier to accurately fit the distribution of the RS amplitudes [8] and we adopt it as a benchmark. In the inset in Fig. 2(d), due to computational limits and thus due to the limited sampling of the minima, the frequency of rare events tends to exhibit discrepancies with respect to the fit for large values of ν; but this will progressively improve if we include more and more minima in the statistics. The mathematical details are reported in Supplement 1. The same can be done for −H > 0 (not shown). We again obtain heavy-tail behavior as expected. We remark that the statistical moments also are similar for the distributions of ν and H . Specifically, we use two parameters that summarize the deviation of the distribution from the Gaussian case: the skewness (third moment)γ and the kurtosis (fourth moment) κ. We computeγ ≈ 4 and κ ≈ 20 so that their product isγκ ≫ 10: this condition characterizes a heavy-tail statistics and corresponds exactly to what was found in previous RS numerical simulations [11]. This gives us further confidence that the WIM model is an appropriate description of the physics behind RS formation in fibers. In summary, we are able to draw a map of the topography of the energy landscape for the soliton interactions, which turns out to comprise many disordered solutions, among which a strongly peaked one emerges.

B. Simulations
Next we simulate light propagation in an optical fiber and compare the characteristic features of rogue wave formation with the exceptional solutions of the WIM. We consider the well-tested and universally used generalized NLS equation where RT 1 − f R δT f R h R T is the nonlinear response function, which includes Kerr and Raman components, Z and T are, respectively, the dimensional propagation distance and time in a frame moving at the group velocity of the input pulse, β k are the dispersion coefficients, γ is the nonlinear parameter, and τ shock represents the shock term, i.e., the coefficient of the first-order approximation of the nonlinearity dispersion [56]. Equation (5) is integrated by setting the following values, similarly to Ref. [8]: a pulse of wavelength λ 0 1064 nm of duration T FWHM 5 ps and peak power of P 100 W is injected into an optical fiber with nonlinear coefficient γ 15Wkm −1 and the chromatic dispersion is expanded up to the third order with coefficients β 2 −0.41 ps 2 ∕km and β 3 0.0687 ps 3 ∕km. We consider a fiber length of 40 m. Finally, we perform 10 3 iterations with different noise realizations (one photon per frequency bin with a random phase).
In order to clearly represent the initial stage of evolution, we first limit our observations to a fiber length of L 25 m and thus identify the two most important scenarios: among all the noise realizations, (i) the output exhibits the strongest peak-we define it as a type-I RS-and (ii) the output contains the most redshifted soliton-we denote it as a type-II RS. The latter corresponds more strictly to the original experimental scenario [4]. We report these two situations in Figs. 3(a) and 3(b), respectively. It is a

Research Article
well-established fact [32] that the strongest peak is not necessarily the most redshifted, due to the role played by collisions in locally enhancing intensity. As described in the following text, we find that the two scenarios correspond to two different kinetic processes in the energy landscape.
After the initial stage of MI and pulse train generation, we observe between 16 and 18 m that collisions start to occur. In the first scenario [see Fig. 3(a)], they do not result in a dramatic splitting of the spectrum into two parts (corresponding to the redshifted soliton and the main part of the pulse train), while many subsequent events lead to the observation of a peak at the fiber output, which is itself a collision. The second scenario [see Fig. 3(b)] entails a more dramatic early event, which transfers a large amount of energy to a single soliton and leads to the appearance of a strongly redshifted spectral part. We show that these sequences of inelastic collisions correspond to hopping from one metastable equilibrium point to another in the energy landscape.
Hereafter, in contrast to Fig. 3, we find that it is crucial to follow the evolution of type-I and type-II RS generation mechanisms beyond L 25 m and to extend our simulations up to L 0 40 m.
In Fig. 4, we show the details of three collision events for each of the scenarios shown in Fig. 3; the time scale is chosen according , which corresponds to a rescaling of distance, L d T 2 0 ∕jβ 2 j, amplitude U 0 γP −1∕2 , and HamiltonianH H ∕H 0 , with H 0 ≡ T 2 0 U 2 0 ∕L d .  This rescaling defined above permits comparing them easily with the solution of the WIM in Fig. 2.
With reference to the type-I event shown in Fig. 3(a), Figs. 4(a) and 4(b) show events occurring at Z ≈ 18 m (the first collision after the pulse train formation) and Z ≈ 25 m, respectively, which lead to the large output peak. Figure 4(c) shows a further collision at a longer propagation distance of Z ≈ 39 m with two main peaks, like the one in Fig. 2(c2).
With reference to Fig. 3(b), a type-II event, Figs. 4(d) and 4(e) deal with the generation of the most redshifted soliton, which stems from a collision at Z ≈ 18 m [Fig. 4(e)] where the power level crosses 2 × 10 3 W [or about 500 in normalized units], about twice as that in the case in Fig. 4(b). This is indeed the main event in the SCG (and the most studied after the original observations in Ref. [4]), as can be noted by comparing it to a collision occurring slightly before [at Z ≈ 17.5 m; Fig. 4(d)] and one at a much further distance, Z ≈ 30 m [ Fig. 4(f)].
We stress the similarity in the amplitude and temporal width of the solutions of the WIM [ Fig. 2(b)] to the temporal profile in Figs Fig. 2(c2)]. The phase profile is strongly distorted by Raman acceleration, i.e., a Zdependent phase slope proportional to the fourth power of the soliton amplitude. Apart from that and possibly an overall slope (due to the choice of the reference frame), the same characteristic jump sequence as in Figs. 2(b) and 2(c) may be recovered.
We thus found that subsequent collisions resemble more and more the solution of the WIM, up to a major collision event [Figs. 4(b) and 4(d)]. The collision at Z ≈ 18 [ Fig. 4(d)], a type-II RS generation, exhibits a much more energetic event than what we observe in Fig. 4(b) in a type-I RS generation. Then minor interactions occur and the RS freely propagates, giving rise to the signature redshifted part of the spectrum first observed experimentally in Ref. [4]. The statistics of the filtered spectra corresponds to the distribution in Fig. 4(d) and this strongly identifies these events with the minima of the energy landscape discussed above. In other words, in a type-II RS, a collision triggers the system into the exceptional point of the landscape and after that the system is nearly quenched in its dynamics. The dynamics of type-I events is apparently more gradual due to a sequence of minor collisions.
The last observation of our work is that after splitting the Hamiltonian into two parts, according to Eq. (4) and as detailed in Supplement 1, a kinetic H K z (due to group velocity dispersion) and a potential H NL (originating from nonlinear interaction) part, see above, we notice that the collision events correspond to the minima of the potential part of the Hamiltonian, H NL z. This is clearly visible in Fig. 5. In order to compare the simulation results with the energy landscape in Fig. 1(c), we normalize the y axis according to the above definitions. While the kinetic part in Fig. 5(a) plays the role of the effective temperature that grows mainly because of Raman acceleration, the potential part H NL z exhibits a sequence of minima. In this regard, the two scenarios considered in Fig. 3 are fundamentally different: in the first [ Fig. 3(a)], the deepest minimum occurs at Z ≈ 25 m, after exploring shallower minima; the second [ Fig. 3(b)] exhibits a much deeper minimum at Z ≈ 18 m, and minor events follow. This is consistent with the collision sequence described above. Moreover, the depth of the first local minimum in the type-II scenario is of the same order of H as obtained for the exceptional solutions of the WIM. The occurrences of these H NL z minima are easily predictable from the definition of Hamiltonian density in Eq. (4); nevertheless, this element reinforces our energy landscape picture and, to our knowledge, was never discussed before. Notice that at the initial stage, at around Z 10 m, the evolution of the two components is nearly opposite: this corresponds approximately to the nonlinear evolution of the MI process, e.g., Akhmediev's breather, which is suddenly broken in favor of the kinetic process, which leads to rogue waves. Finally, we observe that contrary to the expected quadratic growth (see Supplement 1), H K shows only a linear increase: the inelastic collisions among solitons behave as an equivalent friction, which puts an upper limit to the speed of the "particles" in the system. To summarize, it is possible to observe a collision event near the fiber output: this corresponds to a regime characterized by jumps from one of the many shallow minima of the energy  We report the two scenarios described in the text and shown in Fig. 3: the blue solid line corresponds to a type-I RS generation, whereas the red dashed line corresponds to the type-II mechanism. In order to thoroughly characterize the evolution of the two scenarios, we show the evolution up to L 0 40 m and normalize these quantities in order to compare H NL with the energy landscape in Fig. 1(b) (see the main text). We notice that H NL exhibits many local minima, which correspond to step increases in H K and the latter tends to settle on a linear growth. Importantly, each minimum corresponds to a soliton collision. The arrows highlight the main events at Z ≈ 18 m and Z ≈ 25 m, for the two scenarios. The minima correspond to an increase in H K , which then settles on a steady linear growth. The deeper minimum in the type-II scenario leads to a sudden increase to larger values of H K , whereas type-I is characterized by smaller steps.
Research Article landscape to another. On the contrary, in the presence of a largeamplitude redshifted soliton, an early collision lets the system fall inside a very deep minimum in the energy landscape and this is the commonly denoted RS. Upon propagation, Raman effect leads to the growth of the effective temperature of the system, i.e., an increase of kinetic energy, and inhibits the stabilization at one equilibrium point. We have also explored the effect of higher-order dispersion in the absence of the Raman effect. The convective instabilities [23] are undoubtedly an important ingredient for RSs to be excited. Anyway, the collisions in this case correspond to the much shallower minima of H NL (not shown), so that the appearance of extremely strong peaks is substantially hampered.

CONCLUSIONS
We have established a possible link between RSs appearing in SCG via successive soliton collisions and in the minima of the energy landscape that can be described by the WIM, i.e., the most simplified NLS soliton interaction model. We explored a scenario of slow energy buildup, type I, which is explained by hopping among many different shallow equilibria of our energy landscape, and a scenario of quick RS emergence, type-II, in which the system sinks early into a deep minimum. In essence, RS formation is a rare event since slightly different initial configurations for the input noise can lead to very different trajectories in the phase space. When such trajectories intersect the deep minima of the energy landscape, collisions occur that form strongly peaked RSs. The exceptional solution of the WIM represents an extremely rare collision event, whereas other minor equilibria correspond to less dramatic interactions. This also supports the notion that the existence of these equilibria in an integrable model is crucial to the occurrence of extreme events, while other mechanisms such as external forces [42], partial incoherence [17], inhomogeneity [54], or, in our case, Raman scattering are forces that drive the random walk among equilibria. This statistical mechanical framework sheds new light on the numerical and experimental observation that very similar noise configurations lead to completely different pulse dynamics. We conclude that the topography of the energy landscape is thus crucial to explain the formation of RSs and to classify their behavior in terms of collision profiles and statistical properties. Understanding and controlling the energy landscape will lead to a better control of the rogue wave formation. To the best of our knowledge, this result represents the first actual explanation for RS formation in optics, which goes beyond the mere description of their appearance. The generality and universality of the concept will prove to be extremely useful to all communities dealing with the study of freak waves.
German Max Planck Society for the Advancement of Science (MPG). See Supplement 1 for supporting content.