Rogue wave generation by inelastic quasi-soliton collisions in optical fibres

We demonstrate a simple cascade mechanism that drives the formation and emergence of rogue waves in the generalized non-linear Schr\"{o}dinger equation with third-order dispersion. This conceptually novel generation mechanism is based on inelastic collisions of quasi-solitons and is well described by a resonant-like scattering behaviour for the energy transfer in pair-wise quasi-soliton collisions. Our results demonstrate a threshold for rogue wave emergence and the existence of a period of reduced amplitudes - a"calm before the storm"- preceding the arrival of a rogue wave event. Comparing with ultra-long time window simulations of $3.865\times 10^{6}$ps we observe the statistics of rogue waves in optical fibres with an unprecedented level of detail and accuracy, unambiguously establishing the long-ranged character of the rogue wave power-distribution function over seven orders of magnitude.

Here we will show that there is also a process to generate rogue waves through an energy-exchange mechanism when collisions become inelastic in the presence of a third-order dispersion (TOD) term [34,35]. Energyexchange in NLSEs has indeed been observed experimentally [36,37] albeit not yet for TOD. Recent studies confirmed experimentally and numerically that the presence of TOD in optical fibers turns the system convectively unstable and generates extraordinary optical intensities [7,38]. We derive a novel cascade model that simulates the RW generation process directly without the need for a full numerical integration of the gNLSE. The model is validated using a massively parallel simulation [39], allowing us to achieve an unprecedented level of detail through the concurrent use of tens of thousands of CPU cores. Based on statistics from more than 17 × 10 6 interacting quasi-solitons, find that the results of the full gNLSE integration and the cascade model exhibit the same quantitative, long-tail PDF. This agreement highlights the importance of (i) a resonance-like two-soliton scattering coupled with (ii) quasi-soliton energy exchange in giving rise to RWs. We furthermore find that the cascade model and the full gNLSE integration exhibit a "calm-before-the-storm" effect of reduced amplitudes prior to the arrival of the RW, hence hinting towards the possibility of RW prediction. Last, we demonstrate that there is a sharp threshold for TOD to be strong enough to lead to the emergence of RWs.

II. ROGUE WAVES AS CASCADES OF INTER-ACTING SOLITONS
Analytical soliton solutions for the generalized nonlinear Schrödinger equation are only known for β 3 = 0. Here u(z, t) describes a slowly varying pulse envelope, γ the non-linear coefficient, β 2 (< 0) the normal group velocity dispersion and β 3 the TOD [40]. Due to the short distances of only a few kilometers in a super-continuum experiment the effects of absorption, shock term and delayed Raman response can be neglected [11]. Our aim hence is not the most accurate microscopic description possible, but the most simplistic numerical model that is still capable of generating RWs.
Without TOD β 3 , the model (1) can be solved analytically and a u(z, t) describing the celebrated soliton solutions can be found [41]. Depending on the phase difference between two such solitons, attracting or repulsing forces exist between them. This then leads to them either moving through each other unchanged or swapping their positions. The solitons emerge unchanged with the same energy as they had before the collision. This type of collision is elastic and it is the only known type of collision for true analytic (integrable) solitons [26,42]. When we introduce β 3 , the system can no longer be solved analytically and no closed-form solutions are known in general. Numerical integration of (1) shows that stable pulses still exist and these propagate individually much like solitons in the β 3 = 0 case [43]. These quasi-solitons, and collisions between them, are in many cases as elastic as they are for integrable solitons. However, when two quasisolitons have a matching phase, energy can be transferred between them leading to one gaining and the other losing energy [34]. In addition, the emerging pulses have to shed energy through dispersive waves until they have relaxed back to a stable quasi-soliton state [44].
Based on this observation we can now replace the full numerical integration of the gNLSE with a phenomenological cascade model that tracks the collisions between quasi-solitons. Our starting point are quasi-solitons with exponential power tails as found in a super-continuum system after the modulation instability has broken up the initial CW pump laser input into a train of pulses [40,45,46]. We generate a list of such quasi-solitons as our initial condition by randomly choosing the power levels P q , phases φ q and frequency shifts Ω q in accordance with the statistics found in a real system at that point of the integration. The shape of a quasi-soliton is approximated by with effective period and inverse velocity for each quasi-soliton labelled by q. Note that due to (3), the velocity v q of a quasi-soliton depends on its power in addition to the frequency shift for β 3 = 0. We calculate which quasi-solitons will collide first, based on their known initial times t q and velocities v q . The phase difference φ q − φ p between two quasi-solitons is drawn randomly from a uniform distribution in the interval [0, 2π[. Then we calculate the energy transferred from the smaller quasi-soliton to the larger via where ∆E 1 is the energy gain for the higher energy quasisoliton and E 2 is the energy of the second quasi-soliton involved in the collision. A detailed justification of (5) and a discussion of the cross-section coefficient eff will be given below. At this point we estimate the next collision to occur from the set of updated v q and continue as described above until the simulation has reached the desired distance z. Obviously, this procedure is much simpler than a numerical integration of the gNLSE (1). The main strength of the model is a new qualitative and quantitative understanding of RW emergence and dynamics.

III. RARE EVENTS WITH ULTRA-LONG TAILS IN THE PDF
To compare the results with the full numerical integration we generate a power-distribution function (PDF) of the power levels at fixed distances ∆z from u(z, t) = q u q (z, t) using the full quasi-soliton waveform (2). The PDF for the complete set of 17 × 10 6 pulses propagating over 1500m is shown in Fig. 1 (a) and (b) for selected distances for cascade model and gNLSE, respectively, using, for the gNLSE, a highly-optimised, massively parallel and linearly-scaling numerical procedure (see Methods section). After 100m, the PDF exhibits a roughly exponential distribution as seen in Fig. 1(a). With β 3 = 0 this exponential PDF remains stable from this point onwards (cp. inset). However, with β 3 = 0 the inelastic collision of the quasi-solitons leads to an ever increasing number of high-energy RWs. After 500m, a clear deviation from the exponential distribution of the β 3 = 0 case has emerged and beyond 1000m, the characteristic L-shape of a fullydeveloped RW PDF has formed. The PDFs for both gNLSE and the cascade model in Figs  continue to evolve towards higher peak powers with some quasi-solitons becoming larger and larger. In the inset of Fig. 1(b), we compare the long tail behaviour of both PDFs directly. We see that the agreement for PDFs is excellent taking into account that we have reduced the full integration of the gNLSE to only discrete collision events between quasi-solitons. Thus, we find that the essence of the emergence of RWs in this system is very well captured by a process due to inelastic collisions.

IV. MECHANISMS OF THE CASCADE
In Fig. 2, we show a representative example for the propagation of u(z, t), in the full gNLSE and the cascade model, in a short 15ps time range out of the full 3.865×10 6 ps with β 2 = −2.6×10 −28 s 2 m −1 and γ = 0.01 W −1 m −1 . A small initial noise leads to differences in the pulse powers and velocities and hence to eventual collisions of neighbouring pulses. In the enlarged trajectory plots Figs. 2(b) and (c), we see that for β 3 = 0 the solitons interact elastically and propagate on average with the group velocity of the frame. However, for finite β 3 in Fig. 2(c), most collisions are inelastic and one quasi-soliton, with higher energy, moves through the frame from left to right due to its higher energy and group velocity mismatch compared to the frame. It collides in rapid succession with the other quasi-solitons travelling, on average, at the frame velocity. From Fig. 2(b), we see that our cascade model, in which we have replaced the intricate dynamics of the collission by an effective process, similarly shows quasi-solitons starting to collide inelastically, some emerging with higher energies and exhibiting a reduced group velocity. Note that in the cascade model we use an initial pulse power distribution that mimics the PDF of the gNLSE at 100m, and thus only generate data from 100m onwards (cp. Supplement).
The main ingredient of the cascade model, the energy exchange, is obvious in the gNLSE results: in almost all cases, energy is transferred from the quasi-soliton with less energy to the one with more energy leading to the cascade of incremental gains for the more powerful quasi-soliton. This pattern is visible throughout Fig. 2(a) where initial differences in energy of quasisolitons become exacerbated over time and larger and larger quasi-solitons emerge. These accumulate the energy of the smaller ones to the point that the smaller ones eventually vanish into the background. In addition, the group velocity of a quasi-soliton with TOD is dependent on the power of the quasi-soliton [40]. Thus, the emerging powerful quasi-solitons feature a growing group velocity difference to their peers and this increases their collision rate leading to even stronger growth. This can clearly be seen from Fig. 2(a) where larger-energy quasi-solitons start to move sideways as their velocity no longer matches the group velocity of the frame after they have acquired energy from other quasi-solitons due to inelastic collisions. Indeed, the relatively few remaining, soliton-like pulses at 1500m can have peak powers exceeding 1000W. They are truly self-sustaining rogues that have increased their power values by successive interactions and energy exchange with less powerful pulses.

V. CALM BEFORE THE STORM
Looking more closely at the temporal vicinity of waves with particularly large power values, we find that these tend to be preceded by a time period of reduced power values. This "calm before the storm" phenomenon can be observed in Fig. 3. In Fig. 3(a) we can clearly see an asymmetry in the normalized power |u(∆t) relative to the RW event at ∆t = 0 (∆t < 0 denotes events before the RW). The average includes all RWs, defined here as large power events above a threshold of 150W (thresholds 200 and 300W show similar behaviour) and also two independent simulations of the gNLSE, both with parameters as in Fig. 2. The period of calm in power before the RW occurs lasts about ∆t = 1.5ps at z = 200m. It broadens for larger distances, but an asymmetry is retained even at 500m. This finding is further supported by Fig. 3(b), where we note that the "calm before the storm", already observed for the gNLSE, is even clearer and more pronounced for the cascade model. We observe strong oscillations away from ∆t = 0. These describe the simple quasi-soliton pulses which we used to model u(t) in the cascade model. For the gNLSE, these oscillations are much less regular, although still visible. The time interval of the period of calm appears shorter in the cascade model while the amplitude reduction is stronger. We stress that the z = 100m starting position for the cascade model remains a convenient choice.

VI. A THRESHOLD FOR RW EMERGENCE
We find that the emergence of RW behaviour relies on the presence of "enough" β 3 = 0 TOD -or similar additional terms in Eq. (1). Even then, the conditions for the cascade to start are subtle as we show in Fig. 4 for two quasi-soliton collisions with β 3 = 2.64 × 10 −42 s 3 m −1 . The initial conditions for both collisions have been chosen to be identical apart from the relative phase between the two quasi-solitons: an earlier quasi-soliton with initially large power (200W) is met by a later, initially weak pulse (50W). Clearly, the collision of these two quasi-solitons as shown in (a) is (nearly) elastic, they simply exchange their velocities, while retaining their individual power.
This situation retains much of the dynamics from the β 3 = 0 case. In contrast, in (b), we see that after collision three pulses emerge: an early, much weaker quasisoliton (∼ 24W), a later very high power quasi-soliton (∼ 245W), and a final, very weak and dispersive wave (∼ 0.002W) -the collision is highly inelastic. We note that this process is similar to what was described in Refs. 47, 48 for other NLSE variants. Systematically studying many such collisions, we find that the outcome can be modelled quite accurately using Eq. (5) which is similar as in a two-body resonance process. In Fig. 4(c), we show that agreement between Eq. (5) and the numerical simulations is indeed remarkably good. Indeed, (5) was used to generate the cascade model results for Figs. 1 and 2. We have further verified the accuracy of (5) by simulating a large number of individual collision processes with varying relative phases and varying initial quasi-soliton parameters. In Eq. (5), eff is an empirical cross-section coefficient (see the supplement for an analytical justification assuming two-particle scattering). It depends on β 3 as shown in Fig. 4(d). The transition from a regime without RWs to a regime with well pronounced RWs appears rather abrupt. This indicates that already small, perhaps only local changes in β 3 , and hence in the local composition of the optical fibre itself [28,49], can lead to dramatic changes for the emergence of RWs. Indeed, we find that once a RW has established itself in the large-β 3 region, it continues to retain much of its amplitude when entering a β 3 = 0 region.

VII. DISCUSSION AND CONCLUSIONS
Our results emphasize the crucial role played by quasisoliton interactions in the energy exchange underlying the formation of RWs via the proposed cascade mechanism. While interactions are known to play an impor- tant role in RW generation [37,50,51], the elucidation of the full cascade mechanism including its resonance-like quasi-soliton pair scattering and details such as the "calm before the storm", will be essential ingredients of any attempt at RW predictions. In addition, these features are quite different from linear focusing of wave superpositions [32,33] and allow the experimental and observational distinction of both mechanisms.

This leads to the condition
which is in very good agreement with the numerical result of Fig. 4(c). In Eq. (6) the β 3 threshold that leads to fibres supporting RWs depends on the peak power P .
In deriving the numerical estimate in (6) we have used a typical P ∼ 50W as appropriate after about ∼ 100m (cp. Fig. 1 and also the movies in [52]). Once such initial, and still relatively weak RWs have emerged, the condition (6) will remain fulfilled upon further increases in P due to quasi-soliton collisions, indicating the stability of large-peak-power RWs. Our estimation of the effective energy-transfer cross-section parameter eff can of course be improved. However, we believe that (5) captures the essential aspects of the quasi-soliton collisions already very well. Thus far, we have ignored fibre attenuation. Clearly, this would cause energy dissipation and eventually lead to a reduction in the growth of RWs and hence give rise to finite RW lifetimes. But as long as the fibre contains colliding quasi-solitons of enough power, RWs will still be generated.
Up to now, we have used the term RW only loosely to denote high-energy quasi-solitons as shown in Figs. 2 and 1. Indeed, a strict definition of a RW is still an open question and qualitative definitions such as a pulse whose amplitude (or energy or power) is much higher than surrounding pulses are common [28]. Our results now suggest a quantifiable operational definition at least for normal waves in optical fibres described by the gNLSE: a large amplitude wave is not a RW if it occurs as frequently as expected for the PDF at β 3 = 0 (cp. Fig. 1). We emphasise that both high spatial and temporal resolution are required to obtain reliable statistics for RWs in optical fibres for reliable predictions of the PDF. A small time window in the simulation can severely distort the tails of the PDF, and hence their correct interpretation (see supplement).
We find that RWs are preceded by short periods of reduced wave amplitues. This "calm before the storm" has been observed previously [19] in ocean and in optical multifilament RWs, but not yet in studies of optical fibres. We remark that we first noticed the effect in our cascade model, before investigating it in the gNLSE as well. This highlights the usefulness of the cascade model for qualitatively new insights into RW dynamics. More results are needed to ascertain if the periods of calm can be used as reliable predictors for RW occurrence, i.e., reducing false positives.

VIII. METHODS
The numerical simulations of (1) were performed using the split-step Fourier method [40] in the co-moving frame of reference. A massively parallel implementation based on the discard-overlap/save method [53] was implemented to allow for simulations with 2 31 intervals of ∆t = 1.8fs and hence long time windows up to 3.865 × 10 6 ps with several kilometres in propagation distance. We assume periodic boundary conditions in time and, as usual, a coordinate frame moving with the group velocity. The code was shown to scale linearly up to 98k cores (see Supplement).
We start the simulations with a continuous wave of P 0 = 10W power at λ 0 = 1064nm. For the fibre, we assume the parameters β 2 = −2.6 × 10 −28 s 2 m −1 , γ = 0.01 W −1 m −1 , and varying β 3 up to 1.7×(2.64×10 −42 s 3 m −1 ), see Fig. 4. Due to the modulation instability, we observe, after seeding with a small 10 −3 W Gaussian noise, a break-up into individual pulses within the first 100m of the simulation with a density of ∼ 5.88 pulses/ps. Throughout the simulation, we check that the energy re-mains conserved. The PDF of |u| 2 is computed as the simulation progresses. For the two-quasi-soliton interaction study, we use the massively-parallel code as well as a simpler serial implementation. The collision runs are started using pulses of the quasi-soliton shapes (2) with added phase difference exp(iφ) in the advanced pulse. The effective cascade model assumes an initial condition of quasi-solitons of the same density of 5.88 quasisolitons/ps. Their starting times are evenly distributed with separation ∆t ∼ 0.17ps while their initial powers are chosen to mimic the distribution observed in the gNLSE at z ∼ 100m. The time resolution is 2 × 10 −3 ps and we simulate the propagation in 1500 replicas of time windows of 4000ps duration. This gives an effective duration of 6 × 10 6 ps. For all such 30 × 10 6 quasi-soliton pulses, we compute their speeds, find the distance at which the next two-quasi-soliton collision will occur and compute the quasi-soliton energy exchange via (5). The power of the emerging two pulses is calculated from E q = 2P q T q and the algorithm proceeds to find the next collision. The PDF of |u| 2 is computed assuming that the shape of each quasi-soliton is given by (2).

Rogue wave generation due to inelastic quasi-soliton collisions in optical fibres
Marc Eberhard, Antonino Savojardo, Akihiro Maruta and Rudolf A Römer S1. SPECTRA OF RWS FOR β3 = 0 In Fig. S1, we show the spectra |u(λ)| 2 for the gNLSE at long distance z = 1500m for β 3 = 0 and β 3 = 2.64 × 10 −42 s −4 m −1 . The much broader spectrum for β 3 = 0 shows how the TOD has led to wave excitation across a broad range of wave lengths.

A. PDFs for small time window simulations
The ultra-long time window simulations of the gNLSE (1) unambiguously establish the long-range nature of RW PDFs. As shown in Fig. 1, the tails for large powers |u| 2 are clearly visible. In contrast, if the time window of a numerical simulation were chosen too small, it would lead to only a few (and eventually only one) powerful quasi-soliton(s) to extract all the energy from the system. Hence the resulting PDF would have a bump at some high peak power given by the total available energy and reflecting the small numerical system size. This effect is demonstrated in Fig. S2 where we show PDFs derived for a short, 200ps time window simulation, a factor ∼ 20000 shorter than our high-precision parallel simulation. We observe a PDF with an artificial "knee" for high peak powers. Only by increasing the length of the time window, such that even the most powerful quasi-solitons can not pass through the time window more than once throughout a simulation, does one recover the true PDF   (Fig. 1 a).

B. Fitting distributions and their tails
The full PDFs have also been fitted using the Weibull function the PDFs with the respective fits are shown in Fig. S3 (a) for the gNLSE and Fig. S3 (b) for the cascade model. Every fit has been performed taking the log of the PDF(|u| 2 ) and W (|u| 2 ), the resulting coefficients are in Tab. S1. Our results suggest that while a Weibull fit [19] is indeed possible for the tails, a systematic and consistent variation of the fitting parameters with distance travelled is not obvious. While, e.g., the PDF for 200, 500 and 1500m as shown in Fig. 1(a) appears sub-exponential in the tails, we find that the tail of the PDF for 1000m is super-exponential for z 600m. We also fitted the tails of PDF(|u| 2 ) with Both the PDFs and the respective fits are shown in Fig.  S4 while the resulting coefficients are given in Tab. S1. As for the Weibull fit, the results are not convincing and hint towards a continuing development of the shape of the PDF as the propagation continues.

S3. COLLISIONS: DETERMINING THE EN-ERGY GAIN
A. Deriving the energy gain formula In Fig. 4 we observe an energy transfer due to inelastic scattering. This energy gain of quasi-soliton 1 from quasi-soliton 2, can be written as [54] ∆E 1 =ˆG(P 1 , Ω 1 , P 2 , Ω 2 ; z, t, φ)dzdt, where G is an energy density and φ is the phase difference between the two quasi-solitons. It is convenient to change variables (S4) Eq. (S4) shows that a large difference between the (inverses of the) v q results in a reduced energy gain for the larger quasi-soliton in agreement with the results in the main paper. In section VI, we showed the importance of φ for the energy transfer. Fourier-expanding (S4), we can write where (n) are Fourier coefficients (we suppress the P q and Ω q indices for a moment) and φ 0 is the phase difference for which the gain has a maximum. We approximate the above formula with just the first two coefficients. These two coefficients are related; indeed the larger quasi-soliton always gains energy (∆E 1 > 0), hence (0) ≥ | (1)| and because for a certain φ the energy gain is zero, we have (0) = − (1). Thus we can write where we have used 1 − cos (φ) = 2 sin 2 φ 2 and defined = 2 (0). The value of is yet undetermined while the dependence on the group velocities and the phase difference it is clear. As shown in Fig. 4(c), (S6) provides an excellent description of the energy gain in pair-wise quasi-soliton collisions.

B. Effective coupling constant eff
We want to model a system of many colliding quasisolitons as shown in Fig. 2(a). The parameter P1,Ω1,P2,Ω2  Table S1: Fitted coefficients using the Weibull function (W) and F (|u| 2 ) (F).
depends on the individual power and frequency shift of each quasi-soliton pair. In order to devise a tractable model, we have to find an effective eff that describes the average properties of the u-amplitudes well. We therefore choose a distance z = 500, where from Fig. 2(a) we see that well-developed quasi-soliton pulses exist, while the situation is not yet RW dominated as shown in Fig.  1(a). We then use a constant trial value for eff and apply (5) to all quasi-soliton collisions in the cascade model computing data similar to Fig. 2(d) and Fig. 1(b). We then repeat the calculation with another trial eff . For different eff , we compare the PDF created from the cascade model with the PDF obtained from the gNLSE and choose eff such that the agreement is best (see below). We note that this process was followed for the eff values shown in Fig. 4(d) for the different β 3 values. Furthermore, we have checked that similar results can be obtained by using z = 300 as the starting point of the analysis. Once eff is determined, we use it to compute the results for the cascade model, starting at z = 100m and "propagating" all the way to 1500m as described in section VIII. We emphasize the good agreement of the PDFs for z = 500m.

C. Estimating eff
We are interested in determining eff such that the PDFs of the gNLSE and the cascade model (CM) agree at z = 500m. Since we are interested in RWs we want that agreement to be good in the tail region of |u| 2 > 150W. We therefore define the relative variance (S7) and minimize it with respect to eff as shown in Fig.  S5. The eff at minimum is our estimate with accuracy As a second test, we perform a Kolmogorov-Smirnov (KS)-like two-sample test [53] between log PDF gNLSE (|u i | 2 ) and log PDF CM (|u i | 2 , eff ).
We need to renormalize the data count as N i = log(1 + N i ) with each j denoting a |u i | 2 bin and overall N = N bins j=1 log(1 + N i ). Hence the effective number of data is given by Following the KS prescription, we then define the differ- ence and minimize it with respect to eff as shown in Fig. S5.
As usual, with λ = ( N e +0.12+0.11/ N e )D, a KS-like accuracy can be given as although it should no longer be interpreted probabilistically. The results for eff calculated using this KS-like test are also shown in Fig. 4(d) with Q KS given for each β 3 value. It is important to note that Q KS remains roughly constant for all β 3 values, indicating a comparable level of similarity between PDF CM and PDF gNLSE across the full β 3 range. In Fig. S5, we show the eff dependence of the test while Fig. S6 displays the CDFs .

S4. INDEPENDENT QUASI-SOLITONS
In Ref. [19], the authors studied RW events in three data sets, one from ocean waves and two based on optical devices.
They compute, e.g., the time-series autocorrelation function C z (τ ) = dt |u(z, t )| 2 − |u(z, t)| 2 t |u(z, τ − t )| 2 − |u(z, t)| 2 t (see section VI). In Fig. S7, we show C z (τ ) for our gNLSE data as well as for the cascade model at various distances. At 100m quasi-solitons are clearly correlated because of the initial conditions in both models, but from 200m onwards C z (τ ) is nearly zero after a small fraction of picoseconds supporting the notion of largely independently travelling quasi-solitons that interact only when in close spatial and temporal proximity. In agreement with Ref. [19], we hence find a quick decay of C z (t) after z > 100m which supports the notion of well-separated individual quasi-solitons (cp. Fig. S7). Furthermore, the agreement between gNLSE results and the cascade model is very good. Also, the correlation C z (τ ) is close to what has been reported in Ref. [19].

S5. FURTHER DETAILS OF THE CASCADE MODEL
A. Choice of initial PDF As initial condition for the cascade model we compute the PDF of the soliton peak power, P q , in the NLSE case β 3 = 0. We select a distance of z = 1500m such that PDF(|u| 2 ) has stabilized. We find that the resulting PDF(P q ) can be described as where P 0 = 31.4 ± 0.6W and b = 1.69 ± 0.04. The PDF(P q ) and its fit are shown in Fig. S8 (the fit has been performed taking the log of the data and the log of the fitting function). The value P 0 can alternatively be estimated using energy conservation. We start with a continuous wave (CW) power of P CW = 10W. From the autocorrelation C z (τ ), c.p. Fig. S7, we measure the average time between two peaks as ∆T = 0.170 ± 0.008 ps. Hence the initial energy contained in the time window ∆T is E init = P CW ∆T = 1.70 pJ. At distances when quasi-solitons have been created the average energy con- Figure S8: Normalized peak power distribution PDF(P ) for β3 = 0 at 1.5 km. The data points denotes (blue squares) denote the data while the solid (red) line shows the fit with Eq. (S11). The dashed black line is at the fitted value P0 = 31.4 W tained in ∆T , using (3) and (S11), is (S12) From energy conservation, E final = E init , we find In section II, we argued that the shape of a β 3 = 0 quasi-soliton can be approximated by (2). The argument follows Ref. [40]. The solution of (1) can be approximated as a soliton-like pulse where P , T and C represent the amplitude, duration and chirp. The other two parameters are the temporal shift q of the pulse envelope and the frequency shift Ω of the pulse spectrum. The distance-dependence of the parameters can be obtained using the momentum method [40,55,56]. This gives dΩ dz = 0. (S15d) These equations can be solved for C = 0, resulting in (S16) where T and v −1 have been given in (3) and (4).

S6. CALM BEFORE THE STORM
We studied the temporal vicinity of RWs and found that these are preceded by a period of reduced power values. We will refer to this phenomenon as "calm before the storm" [19]. First we consider the intensity |u| 2 at a certain distance z, then we find quasi-solitons with a peak power P larger than a certain minimum power P min . Next, we consider the time window that goes from 2ps before to 2ps after every of these events. Finally we take the average of all the events found. We choose P min = 150W, while the distances z go from 150m to 500m for the gNLSE and up to 1500m for the cascade model. The results for the gNLSE and the cascade model are shown in Fig. S9. In both cases we can see a dip in the normalized power |u(∆t)| 2 / |u(∆t)| 2 before the RW event at ∆t = 0. In the cascade model, the period of calm is more pronounced, but for longer distances (z ≥ 1000m, light gray lines in the plot) the intensity starts to flatten as in the gNLSE case. For distances z ≤ 200m, we observe strong oscillations away from ∆t = 0. In the cascade model these are more regular than in the gNLSE. The cause is in the initial conditions of the cascade model, where at the beginning quasi-solitons are equally spaced, resulting in the observed correlations in the intensity. For z > 200m the oscillations disappear because the quasi-solitons are becoming more uncorrelated. The number of RWs depends on the minimum power chosen, therefore the best statistics is for smaller RWs with, e.g., P min =150W, while there are more fluctuations for P min = 300W. At large distances the statistic improves because more and more RWs are produced.

S7. PARALLEL SCALING
The results presented here were run on various massively parallel high-performance computing (HPC) machines and architectures. These included HPC clusters at Warwick's Centre for Scientific Computing, the UK national facilities HECToR and ARCHER as well as the BlueGene/Q machine of the Hartree Centre. As explained in section VIII, the code was MPI-parallized and  can take optimal advantage of these HPC architectures. In Fig. S10, we show the speedup and efficiency of a test run on the 98k cores at the Hartree Centre's BG/Q. The results show linear scaling and nearly 100% efficiency.