1. Introduction

Historically, reports of “monster” or “freak” waves [1–3] on the earth’s oceans have been seen largely as sea men’s tales [4,5]. However, the recent availability of reliable experimental observations [4, 6] has proved their existence and shown that these “rogues” are indeed rare events [7], governed by long tails in their probability distribution function (PDF) [8], and hence concurrent with very large wave amplitudes [9, 10]. As both deep water waves in the oceans and optical waves in fibres can be described by similar generalized non-linear Schrödinger equations (gNLSE) they both show rogue waves (RW) and long-tail statistics [8, 11, 12]. The case of RW generation in optical fibres during super-continuum generation has been observed experimentally [13–16]. Recently, experimental data of long tails in the PDF have been collected [17,18], as well as time correlations in various wave phenomena with RW occurrence studied [19]. RWs and long-tailed PDFs have also been found during high power femtosecond pulse filamentation in air [20], in non-linear optical cavities [21] and in the output intensity of optically injected semiconductors laser [22], mode-locked fiber lasers [23], Raman fiber lasers [24] and fiber Raman amplifiers [25]. However, it still remains largely unknown how RWs emerge [26–28] and theoretical explanations range from high-lighting the importance of the non-linearity [9, 29–31] to those based on short-lived linear superpositions of quasi-solitons during collisions [32,33].

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]. Energy-exchange in NLSEs has indeed been observed experimentally [36,37] and numerically [38–41] although no quantitative description of the process has been given. 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,42]. Energy exchange mechanisms have also been observed in the context of optical waves turbulence [43–45]. We derive a 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 [46], 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⁶ 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.

2. Rogue waves as cascades of interacting solitons

Analytical soliton solutions for the generalized non-linear Schrödinger equation

Without TOD β₃, the model (1) can be solved analytically and a u(z, t) describing the celebrated soliton solutions can be found [48]. 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,49]. When we introduce β₃, the system can no longer be solved analytically and no closed-form solutions are known in general. Numerical integration of Eq. (1) shows that stable pulses still exist and these propagate individually much like solitons in the β₃ = 0 case [50]. These quasi-solitons, and collisions between them, are in many cases as elastic as they are for integrable solitons. However, when two quasi-solitons 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 [51].

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 [47,52,53]. 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

3. 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)={\sum }_q{u}_q(z,t)$ using the full quasi-soliton waveform (2). The PDF for the complete set of ≳ 17 × 10⁶ pulses propagating over 1500m is shown in Fig. 1(a) and 1(b) for selected distances for cascade model and gNLSE, respectively, using, for the gNLSE, a highly-optimised, massively parallel and linearly-scaling numerical procedure [46]. After 100m, the PDF exhibits a roughly exponential distribution as seen in Fig. 1(a). With β₃ = 0 this exponential PDF remains stable from this point onwards (cp. inset). However, with β₃ ≠ 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 β₃ = 0 case has emerged and beyond 1000m, the characteristic L-shape of a fully-developed RW PDF has formed. The PDFs for both gNLSE and the cascade model in Figs. 1(a) and 1(b) then continue to evolve towards higher peak powers with some quasi-solitons becoming larger and larger. We remark that while fits descriptions of the PDF with Weibull, Pareto and stretched exponential are possible, no such fit is convincing across the full range 0–1000 W (see Figs. 1 (a) and 1(b) and Appendix C). 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.

 

Fig. 1 (a) PDFs of the intensity |u|² from the gNLSE (1) at β₃ = 2.64 × 10⁻⁴²s³m⁻¹ using a large time window of Δt = 3.865 × 10⁶ps. The PDFs have been computed at distances z = 100m, 200m, 500m, 1000m and 1500m. The left vertical axis denotes the values of the normalized PDF while the right vertical axis gives the event count per bin. The PDFs have been fitted using the Weibull function (full black lines). The brown points represnt the PDF calculate for a small time window of 200ps at 1500m. The inset shows results for β₃ = 0. (b) PDFs of the intensity |u|² from the cascade model for the same distances as in (a), using the same symbol and axes conventions. The PDFs have been fitted using the Weibull function (full black lines). The inset shows a comparison between the results from the gNLSE (colored lines) and the cascade model (black lines and symbol outlines) for z = 500m, 1000m and 1500m. Only every 50th symbol is shown.

Download Full Size | PDF

4. 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⁶ ps with β₂ = −2.6×10⁻²⁸s²m⁻¹ and γ = 0.01 W⁻¹m⁻¹. 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 2(c), we see that for β₃ = 0 the solitons interact elastically and propagate on average with the group velocity of the frame. However, for finite β₃ 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. appendix).

 

Fig. 2 (a) Intensity |u(z, t)|² for β₃ = 2.64 × 10⁻⁴²s³m⁻¹ of the gNLSE Eq. (1) as function of the time t and distance z in a selected time frame of Δt = 15ps and distance range Δz = 1.5km. (b) |u|² with β₃ = 0 for a zoomed-in distance and time region, (c) |u|² with β₃ value as in (a) for a region of (a) with Δt and Δz chosen identical to (b). (d) Intensities |u|² as computed from the effective cascade model using the same shading/color scale as in (a). Note that we start the effective model at z₀ = 100m to mimic the effects of the modulation instability in (a).

Download Full Size | PDF

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 quasi-solitons 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 [47]. 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.

5. A threshold for RW emergence

We find that the emergence of RW behaviour relies on the presence of “enough” β₃ ≠ 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. 3 for two quasi-soliton collisions with β₃ = 2.64 × 10⁻⁴²s³m⁻¹.

 

Fig. 3 (a+b) Intensities |u(z, t)|² of the scattering between two quasi-solitons. The phase difference ϕ was chosen to correspond to (a) the minimum and (b) the maximum of ΔE₁/E₂(ϕ). (c) ΔE₁/E₂(ϕ) for various choices of initial speeds. The data points represents results of the gNLSE (1) while the lines denote the fit (5). (d) ∊ values have been obtained comparing the PDF from the gNLSE and the cascade model at short distance using a Kolmogorov-Smirnov-like (KS) test and minimizing the variance (cf. appendix B). A KS significance measure is also shown [60].

Download Full Size | PDF

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 Fig. 3(a) is (nearly) elastic, they simply exchange their positions, while retaining their individual power. This situation retains much of the dynamics from the β₃ = 0 case. In contrast, in Fig. 3(b), we see that after collision three pulses emerge: an early, much weaker quasi-soliton (∼ 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. [54, 55] 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. We emphasize that although the dispersive wave is fundamental for the energy transfer, its peak power is 5 orders of magnitude smaller than the solitons power, therefore negligible for the implementation of the effective model.

In Fig. 3(c), we show that agreement between Eq. (5) and the numerical simulations is indeed remarkably good. Indeed, Eq. (5) was used to generate the cascade model results for Figs. 1 and 2. We have further verified the accuracy of Eq. (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 appendix for an analytical justification assuming two-particle scattering). It depends on β₃ as shown in Fig. 3(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 β₃, and hence in the local composition of the optical fibre itself [28,56], can lead to dramatic changes for the emergence of RWs. Indeed, we find that once a RW has established itself in the large-β₃ region, it continues to retain much of its amplitude when entering a β₃ = 0 region.

6. 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 [19] can be observed in Fig. 4. In Fig. 4(a) we can clearly see an asymmetry in the normalized power |u(Δt)|²/〈|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. Physically, this effect can be understood as follows: a RW moves more slowly than the non-RW solitons. Upon interaction with a soliton, the RW gains in energy, but due to Eq. (3) and Eq. (4), slows down even more (cp. Fig. 2); the soliton, having overtaken the RW, has lost some of its energy, therefore leading to a reduction in intensity before the RW.

 

Fig. 4 Normalized averaged powers |u(Δt)|²/〈|u(Δt)|²〉 for times Δt in the vicinity of a RW event at Δt = 0. Panel (a), (b) corresponds to 200 and 500m, respectively. Solid lines in both panels indicate averaged results for two gNLSE runs (with parameters as in Fig. 2), while dashed lines show the corresponding results for the cascade model. In both panels, we identify RWs as corresponding to powers equal to or larger than 150W. The colours are chosen to indicate distances compatible with a full set of results z = 150, . . ., 1500m given in the appendix. Note that |u(0)|²/〈|u(0)|²〉 > 10 in both panels.

Download Full Size | PDF

This finding is further supported by Fig. 4(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 note that the apparently more regular oscillations in the cascade model at 200m is an artefact of our starting condition with solitons chosen perfectly equidistant at z = 100m to reproduce the same density as in the gNLSE.

7. Methods

The numerical simulations of Eq. (1) were performed using the split-step Fourier method [47] in the co-moving frame of reference. A massively parallel implementation based on the discard-overlap/save method [60] was implemented to allow for simulations with 2³¹ intervals of Δt = 1.8fs and hence long time windows up to 3.865 × 10⁶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 scales linearly up to 98k cores.

We start the simulations with a continuous wave of P ₀ = 10W power at λ₀ = 1064nm. For the fibre, we assume the parameters β₂ = −2.6 × 10⁻²⁸s²m⁻¹, γ = 0.01 W⁻¹m⁻¹, and varying β₃ up to 1.7 × (2.64 × 10⁻⁴²s³m⁻¹), see Fig. 3. Due to the modulation instability, we observe, after seeding with a small 10⁻³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 remains conserved. The PDF of |u|² 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 quasi-solitons/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⁻³ps and we simulate the propagation in 1500 replicas of time windows of 4000ps duration. This gives an effective duration of 6 × 10⁶ps. For all such 30 × 10⁶ 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_qT_q and the algorithm proceeds to find the next collision. The PDF of |u|² is computed assuming that the shape of each quasi-soliton is given by (2).

8. Conclusions

Our results emphasize the crucial role played by quasi-soliton interactions in the energy exchange underlying the formation of RWs via the proposed cascade mechanism. While interactions are known to play an important role in RW generation [37–41,43–45,57,58], the elucidation of the full cascade mechanism including its resonance-like quasi-soliton pair scattering and details such as the “calm before the storm”, might 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.

RWs emerge when β₃ is large enough as shown in Fig. 3. Their appearance is very rapid in a short range 0.8 ≲ β₃/2.64 × 10⁻⁴²s³m⁻¹ ≲ 1. This can be understood as follows: the dispersion relation, in the moving frame, is β(ω) = β₂(ω − ω₀)²/2 + β₃(ω − ω₀)³/6. The anomalous dispersion region of β(ω) < 0, with soliton-like excitations, ends at ω_c − ω₀ ≥ −3β₂/β₃ beyond which β(ω) ≥ 0 and dispersive waves emerge. From (3), we can estimate the spectral width as $2({\omega }_c-{\omega }_0)=2\pi \sqrt{\gamma P/\left|{\beta }_2\right|}$. This leads to the condition

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, in agreement with recent work [28], 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 β₃ = 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.

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. Last, we emphasise that the underlying cause for the RW emergence due to TOD or Raman terms is of course the absence of integrability. Our results suggest that, e.g., a suitable dispersion engineering for Eq. (1) even without TOD (or Raman term) could destabilize the solitons and lead eventually to the formation of RWs.

A. Further details of the cascade model

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 β₃ = 0. We select a distance of z = 1500m such that PDF(|u|²) has stabilized. We find that the resulting PDF(P_q) can be described as

(7)$$\rho (P_q)=\frac{b}{P_0}{\left(\frac{P_q}{P_0}\right)}^{b-1}\text{exp}\left[-{\left(\frac{P_q}{P_0}\right)}^b\right],$$

where P₀ = 31.4±0.6W and b = 1.69±0.04. The PDF(P_q) and its fit are shown in Fig. 5(a) (the fit has been performed taking the log of the data and the log of the fitting function). The value P₀ can alternatively be estimated using energy conservation. We start with a continuous wave (CW) power of P_CW = 10W. From the autocorrelation, 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 contained in ΔT, using Eq. (3) and Eq. (7), is $E_{\text{final}}={\int }_0^{\infty }2P_q{T}_q(P_q)\rho (P_q)d{P}_q\simeq 1.796\sqrt{\left|{\beta }_2\right|P_0/\gamma }$. From energy conservation, E_final = E_init, we find P₀ ≃ γ (P_CWΔT/1.796)² /|β₂| = 34 ± 3W.

 

Fig. 5 (a) Normalized peak power distribution PDF(P) for β₃ = 0 at 1.5 km. The data points denotes (blue squares) denote the data while the solid (red) line shows the fit with Eq. (7). The dashed black line is at the fitted value P₀ = 31.4 W (b)Relative variance r (filled squares) and largest difference D (open circles) calculated for different values of ∊_eff at β₃ = 2.64 × 10⁻⁴²s³m⁻¹. Parabolic fits to the data are shown as lines. The vertical dotted line denotes the estimated ∊_eff = (1.23 ± 0.05)fs/m at which D is minimal, the grey region indicates the error of that estimate. The vertical dashed-dotted line denotes the estimate ∊_eff = (1.32 ± 0.05)fs/m from r.

Download Full Size | PDF

Derivation of the effective quasi-soliton description

In section 2, we argued that the shape of a β₃ ≠ 0 quasi-soliton can be approximated by Eq. (2). The argument follows Ref. [47]. The solution of Eq. (1) can be approximated as a soliton-like pulse

(8)$$u(z,t)=\sqrt{P}\text{sech}\left[\frac{t-q(z)}{T}\right]\times \text{exp}\left\{-i\mathrm{\Omega }\left[t-q(z)\right]-iC\frac{{\left[t-q(z)\right]}^2}{2T^2}\right\},$$

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 [47,62,63]. This gives

(9)$$\begin{array}{ll}\frac{dT}{dz}=({\beta }_2+{\beta }_3\mathrm{\Omega })\frac{C}{T},\hfill & \frac{dC}{dz}=\left(\frac{4}{{\pi }^2}+C^2\right)\frac{({\beta }_2+{\beta }_3\mathrm{\Omega })}{T^2}+\frac{4\gamma P}{{\pi }^2},\hfill \\ \frac{dq}{dz}={\beta }_2\mathrm{\Omega }+\frac{{\beta }_3}{2}{\mathrm{\Omega }}^2+\frac{{\beta }_3}{6T^2}\left(1+\frac{{\pi }^2}{4}{C}^2\right),\hfill & \frac{d\mathrm{\Omega }}{dz}=0.\hfill \end{array}$$

The equations can be solved for C = 0, resulting in Eq. (2) with ϕ = −Ω(t − v⁻¹z), while T and v⁻¹ have been given in Eq. (3) and Eq. (4).

B. Determining the energy gain

Deriving the energy gain formula

In Fig. 3 we observe an energy transfer due to inelastic scattering. This energy gain of quasi-soliton 1 from quasi-soliton 2, can be written as ΔE₁ = ∫ 𝒢(P₁, Ω₁, P₂, Ω₂; z, t, ϕ)dzdt [61], where 𝒢 is an energy density and ϕ is the phase difference between the two quasi-solitons. It is convenient to change variables s = t − z/v₁, w = t − z/v₂. Such that

(10)$$\mathrm{\Delta }{E}_1=\frac{1}{\left|v_1^{-1}-v_2^{-1}\right|}\int 𝒢\left(P_1,{\mathrm{\Omega }}_1,P_2,{\mathrm{\Omega }}_2;s,w,\varphi \right)dsdw.$$

Eq. (10) 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. Fourier-expanding Eq. (10), we can write

(11)$$\frac{\mathrm{\Delta }{E}_1}{E_2}=\frac{1}{\left|v_1^{-1}-v_2^{-1}\right|}\sum _{n=0}^{\infty }{∊}_{P_1,{\mathrm{\Omega }}_1,P_2,{\mathrm{\Omega }}_2}(n)\text{cos}\left[n(\varphi -{\varphi }_0)\right],$$

where ∊(n) are Fourier coefficients (we suppress the P_q and Ω_q indices for a moment) and ϕ₀ 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₁ > 0), hence ∊(0) ≥ |∊(1)| and because for a certain ϕ the energy gain is zero, we have ∊(0) = −∊(1). Thus we can write

(12)$$\frac{\mathrm{\Delta }{E}_1}{E_2}\simeq \frac{{∊}_{P_1,{\mathrm{\Omega }}_1,P_2,{\mathrm{\Omega }}_2}}{\left|v_1^{-1}-v_2^{-1}\right|}{\text{sin}}^2\left(\frac{\varphi -{\varphi }_o}{2}\right),$$

where we have used $1-\text{cos}(\varphi )=2{\text{sin}}^2\left(\frac{\varphi }{2}\right)$ and defined ∊ = 2∊(0). The value of ∊ is yet undetermined while the dependence on the group velocities and the phase difference is clear. As shown in Fig. 3(c), Eq. (12) provides an excellent description of the energy gain in pair-wise quasi-soliton collisions. In order to determine the energy transfer we estimate the individual pulse energies in the gNLSE via ${\int }_{\mathrm{\Delta }t}{\left|u(z,t)\right|}^{2}dt$ with Δt = 1ps symmetric w.r.t. the maximum peak power. In the cascade model, we use E_q = 2P_qT_q. In principle, a collision can shift the soliton frequency Ω → Ω + ΔΩ and hence the dispersion via β₂ + β₃ΔΩ. The relative correction is ∼ (β₃/β₂)ΔΩ ∝ 10⁻²ΔΩ THz⁻¹. In the gNLSE simulations, we observe at most ΔΩ ≤ 3THz, and the relative variation will be ≤ 3%. We thus ignore the effect and only consider the change in energy due to the peak power variation.

Effective coupling constant ∊_eff

We want to model a system of many colliding quasi-solitons as shown in Fig. 2(a). The parameter ∊_{P1,Ω1,P2,Ω2} 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 use a constant trial value for ∊_eff and apply Eq. (5) to all quasi-soliton collisions in the cascade model computing data similar to Fig. 2(d) and Fig. 1(b). We 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. 3(d) for the different β₃ 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 7. We emphasize the good agreement of the PDFs for z ≠ 500m.

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|² > 150W. We therefore define the relative variance

(13)$$r({∊}_{\text{eff}})=\frac{{\sum }_i{\left[\text{log}\hspace{0.17em}{\text{PDF}}_{\text{gNLSE}}\left({\left|u_i\right|}^2\right)-\text{log}\hspace{0.17em}{\text{PDF}}_{\text{CM}}\left({\left|u_i\right|}^2,{∊}_{\text{eff}}\right)\right]}^2}{{\sum }_i{\left[\text{log}\hspace{0.17em}{\text{PDF}}_{\text{gNLSE}}\left({\left|u_i\right|}^2\right)\right]}^2}$$

and minimize it with respect to ∊_eff as shown in Fig. 5. The ∊_eff at minimum is our estimate with accuracy ${∊}_{\text{eff}}\sqrt{r({∊}_{\text{eff}})}$. The results for ∊_eff calculated using this variance minimization are shown in Fig. 3(d) (red line).

As a second test, we perform a Kolmogorov-Smirnov (KS)-like two-sample test [60] between log PDF_gNLSE(|u_i|²) and log PDF_CM(|u_i|², ∊_eff). We need to renormalize the data count as Ñ_i = log(1+ N_i) with each j denoting a |u_i|² bin and overall $\tilde{N}={\sum }_{j=1}^{N_{\text{bins}}}\text{log}(1+N_i)$. Hence the effective number of data is given by ${\tilde{N}}_e=\left({\tilde{N}}_{\text{gNLSE}}{\tilde{N}}_{\text{CM}}\right)/\left({\tilde{N}}_{\text{gNLSE}}+{\tilde{N}}_{\text{CM}}\right)$. Following the KS prescription, we then define the difference $D=\underset{i\in \left[1,N_{\text{bins}}\right]}{\text{max}}\left|{\text{CDF}}_{\text{gNLSE}}(i)-{\text{CDF}}_{\text{CM}}(i)\right|$ and minimize it with respect to ∊_eff as shown in Fig. 5 (b). As usual, with $\lambda =\left(\sqrt{{\tilde{N}}_e}+0.12+0.11/\sqrt{{\tilde{N}}_e}\right)D$, a KS-like accuracy can be given as $Q_{\text{KS}}(\lambda )=2{\sum }_{j=1}^{\infty }{\left(-1\right)}^{j-1}{e}^{-2j^2{\lambda }^2}$, although it should no longer be interpreted probabilistically. The results for ∊_eff calculated using this KS-like test are also shown in Fig. 3(d) with Q_KS given for each β₃ value. It is important to note that Q_KS remains roughly constant for all β₃ values, indicating a comparable level of similarity between PDF_CM and PDFgNLSE across the full β₃ range. In Fig. 5 (b), we show the ∊_eff dependence of the test.

C. Fitting the PDFs

The PDFs of Fig. 1 have been fitted using the Weibull function

(14)$$W\left({\left|u\right|}^2\right)=b{a}^{-b}{\left({\left|u\right|}^2\right)}^{b-1}\text{exp}\left[-{\left(\frac{{\left|u\right|}^2}{a}\right)}^b\right]$$

The fits are shown in Fig. 1. Every fit has been performed taking the log of the PDF(|u|²) and W(|u|²), the resulting coefficients are in Table 1. 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|²) with a stretched exponential

(15)$$F\left({\left|u\right|}^2\right)=F_0\text{exp}\left[-{\left(\frac{{\left|u\right|}^2}{a}\right)}^b\right],$$

and, following Ref. [40], with a Pareto function

(16)$$Q\left({\left|u\right|}^2\right)=\frac{a{b}^a}{{\left|u\right|}^{2(a+1)}}.$$

The resulting coefficients are given in Table 1. For all fits, the results are not convincing as documented by the rather large χ² values found. This hints towards a continuing development of the shape of the PDF as the propagation continues.

Table 1. Values of the coefficients a, b and F₀ needed to fit the various PDF(|u|²) distributions from Fig. 1 at three representative distances z = 500 m, 1000 m and 1500 m, using (W) the Weibull function of Eq. (14), (F) the stretched exponential of Eq. (15) and (P) the Pareto function as in Eq. (16). The reduced χ² value is also displayed.

View Table

Acknowledgments

We thank George Rowlands for stimulating discussions in the early stages of the project. We are grateful to the EPSRC for provision of computing resources through the MidPlus Regional HPC Centre (EP/K000128/1), and the national facilities HECToR (e236, ge236) and ARCHER (e292). We thank the Hartree Centre for use of its facilities via BG/Q access projects HCBG055, HCBG092, HCBG109.

References and links

1. L. Draper, “‘Freak’ ocean waves,” Oceanus 10, 13–15 (1964).

2. L. Draper, “Severe wave conditions at sea,” J. Inst. Navig. 24, 273–277 (1971). [CrossRef]  

3. J. Mallory, “Abnormal waves in the south-east coast of South Africa,” Int. Hydrog. Rev. 51, 89–129 (1974).

4. S. Perkins, “Dashing Rogues,” Science News 170, 328 (2006). [CrossRef]  

5. M. Erkintalo, “Rogue waves: Predicting the unpredictable?” Nat. Photonics 9, 560–562 (2015). [CrossRef]  

6. M. Hopkin, “Sea snapshots will map frequency of freak waves,” Nature 430, 492 (2004). [CrossRef]   [PubMed]  

7. C. Kharif and E. Pelinovsky, “Physical mechanisms of the rogue wave phenomenon,” Eur. J. Mech. B/Fluids 22, 603–634 (2003). [CrossRef]  

8. M. Onorato, A. R. Osborne, M. Serio, L. Cavaleri, C. Brandini, and C. T. Stansberg, “Observation of strongly non-Gaussian statistics for random sea surface gravity waves in wave flume experiments,” Phys. Rev. E 70, 067302 (2004). [CrossRef]  

9. T. A. A. Adcock, P. H. Taylor, and S. Draper, “Nonlinear dynamics of wave-groups in random seas: unexpected walls of water in the open ocean,” Proc. Roy. Soc. A: Math. Phys. Eng. Sci. 471, 20150660 (2015). [CrossRef]  

10. M. Onorato, S. Residori, U. Bortolozzo, A. Montina, and F. Arecchi, “Rogue waves and their generating mechanisms in different physical contexts,” Phys. Rep. 528, 47–89 (2013). [CrossRef]  

11. J. M. Dudley, F. Dias, M. Erkintalo, and G. Genty, “Instabilities, breathers and rogue waves in optics,” Nature 8, 755–764 (2014).

12. A. Chabchoub, N. Hoffmann, M. Onorato, G. Genty, J. M. Dudley, and N. Akhmediev, “Hydrodynamic supercontinuum,.” Phys. Rev. Lett. 111, 054104 (2013). [CrossRef]   [PubMed]  

13. D. R. Solli, C. Ropers, P. Koonath, and B. Jalali, “Optical rogue waves,” Nature 450, 1054–1057 (2007). [CrossRef]   [PubMed]  

14. D. R. Solli, C. Ropers, and B. Jalali, “Active control of rogue waves for stimulated supercontinuum generation,” Phys. Rev. Lett. 101, 233902 (2008). [CrossRef]   [PubMed]  

15. M. Erkintalo, G. Genty, and J. M. Dudley, “Rogue-wave-like characteristics in femtosecond supercontinuum generation,” Opt. Lett. 34, 2468 (2009). [CrossRef]   [PubMed]  

16. B. Kibler, C. Finot, and J. M. Dudley, “Soliton and rogue wave statistics in supercontinuum generation in photonic crystal fibre with two zero dispersion wavelengths,” Eur. Phys. J. Spec. Top. 173, 289–295 (2009). [CrossRef]  

17. S. Randoux, P. Walczak, M. Onorato, and P. Suret, “Intermittency in integrable turbulence,” Phys. Rev. Lett. 113, 113902 (2014). [CrossRef]   [PubMed]  

18. P. Walczak, S. Randoux, and P. Suret, “Optical rogue waves in integrable turbulence,” Phys. Rev. Lett. 114, 143903 (2015). [CrossRef]   [PubMed]  

19. S. Birkholz, C. Brée, A. Demircan, and G. Steinmeyer, “Predictability of Rogue Events,” Phys. Rev. Lett. 114, 213901 (2015). [CrossRef]   [PubMed]  

20. J. Kasparian, P. Béjot, J.-P. Wolf, and J. M. Dudley, “Optical rogue wave statistics in laser filamentation,” Opt. Express 17, 12070 (2009). [CrossRef]   [PubMed]  

21. A. Montina, U. Bortolozzo, S. Residori, and F. T. Arecchi, “Non-Gaussian statistics and extreme waves in a nonlinear optical cavity,” Phys. Rev. Lett. 103, 173901 (2009). [CrossRef]   [PubMed]  

22. C. Bonatto, M. Feyereisen, S. Barland, M. Giudici, C. Masoller, J. R. R. Leite, and J. R. Tredicce, “Deterministic optical rogue waves,” Phys. Rev. Lett. 107, 053901 (2011). [CrossRef]   [PubMed]  

23. C. Lecaplain, P. Grelu, J. M. Soto-Crespo, and N. Akhmediev, “Dissipative rogue waves generated by chaotic pulse bunching in a mode-locked laser,” Phys. Rev. Lett. 108, 233901 (2012). [CrossRef]   [PubMed]  

24. S. Randoux and P. Suret, “Experimental evidence of extreme value statistics in Raman fiber lasers,” Opt. Lett. 37, 500–502 (2012). [CrossRef]   [PubMed]  

25. K. Hammani, C. Finot, J. M. Dudley, and G. Millot, “Optical rogue-wave-like extreme value fluctuations in fiber Raman amplifiers,” Opt. Express 16, 16467 (2008). [CrossRef]   [PubMed]  

26. N. Akhmediev and E. Pelinovsky, “Editorial – Introductory remarks on “Discussion &Rogue waves–towards a unifying concept: discussion and debates,” Eur. Phys. J. Spec. Top. 185, 1–4 (2010). [CrossRef]  

27. V. Ruban, Y. Kodama, M. Ruderman, J. Dudley, R. Grimshaw, P. V. E. McClintock, M. Onorato, C. Kharif, E. Pelinovsky, T. Soomere, G. Lindgren, N. Akhmediev, A. Slunyaev, D. Solli, C. Ropers, B. Jalali, F. Dias, and A. Osborne, “Rogue waves – towards a unifying concept?: Discussions and debates,” Eur. Phys. J. Spec. Top. 185, 5–15 (2010). [CrossRef]  

28. N. Akhmediev, B. Kibler, F. Baronio, M. Belić, W.-P. Zhong, Y. Zhang, W. Chang, J. M. Soto-Crespo, P. Vouzas, P. Grelu, C. Lecaplain, K. Hammani, S. Rica, A. Picozzi, M. Tlidi, K. Panajotov, A. Mussot, A. Bendahmane, P. Szriftgiser, G. Genty, J. Dudley, A. Kudlinski, A. Demircan, U. Morgner, S. Amiraranashvili, C. Bree, G. Steinmeyer, C. Masoller, N. G. R. Broderick, A. F. J. Runge, M. Erkintalo, S. Residori, U. Bortolozzo, F. T. Arecchi, S. Wabnitz, C. G. Tiofack, S. Coulibaly, and M. Taki, “Roadmap on optical rogue waves and extreme events,” J. Opt. 18, 063001 (2016). [CrossRef]  

29. B. Kibler, J. Fatome, C. Finot, G. Millot, F. Dias, G. Genty, N. Akhmediev, and J. M. Dudley, “The Peregrine soliton in nonlinear fibre optics,” Nat. Phys. 6, 790–795 (2010). [CrossRef]  

30. A. Chabchoub, N. P. Hoffmann, and N. Akhmediev, “Rogue wave observation in a water wave tank,” Phys. Rev. Lett. 106, 204502 (2011). [CrossRef]   [PubMed]  

31. B. Kibler, J. Fatome, C. Finot, G. Millot, G. Genty, B. Wetzel, N. Akhmediev, F. Dias, and J. M. Dudley, “Observation of Kuznetsov-Ma soliton dynamics in optical fibre,” Sci. Rep. 2, 463 (2012). [CrossRef]   [PubMed]  

32. M. Onorato, A. R. Osborne, M. Serio, and L. Cavaleri, “Modulational instability and non-Gaussian statistics in experimental random water-wave trains,” Phys. Flu. 17, 078101 (2005). [CrossRef]  

33. R. Höhmann, U. Kuhl, H.-J. Stöckmann, L. Kaplan, and E. J. Heller, “Freak waves in the linear regime: A microwave study,” Phys. Rev. Lett. 104, 093901 (2010). [CrossRef]   [PubMed]  

34. N. Akhmediev, J. M. Soto-Crespo, and A. Ankiewicz, “Could rogue waves be used as efficient weapons against enemy ships?” Eur. Phys. J. Spec. Top. 185, 259–266 (2010). [CrossRef]  

35. G. Weerasekara, A. Tokunaga, H. Terauchi, M. Eberhard, and A. Maruta, “Soliton’s eigenvalue based analysis on the generation mechanism of rogue wave phenomenon in optical fibers exhibiting weak third order dispersion,” Opt. Express 23, 143 (2015). [CrossRef]   [PubMed]  

36. F. Luan, D. V. Skryabin, A. V. Yulin, and J. C. Knight, “Energy exchange between colliding solitons in photonic crystal fibers,” Opt. Express 14, 9844 (2006). [CrossRef]   [PubMed]  

37. A. Mussot, A. Kudlinski, M. Kolobov, E. Louvergneaux, M. Douay, and M. Taki, “Observation of extreme temporal events in CW-pumped supercontinuum,” Opt. Express 17, 17010–17015 (2009). [CrossRef]   [PubMed]  

38. A. Armaroli, C. Conti, and F. Biancalana, “Rogue solitons in optical fibers: a dynamical process in a complex energy landscape?” Optica 2, 497 (2015). [CrossRef]  

39. C. Brée, G. Steinmeyer, I. Babushkin, U. Morgner, and A. Demircan, “Controlling formation and suppression of fiber-optical rogue waves,” Opt. Lett. 41, 3515 (2016). [CrossRef]   [PubMed]  

40. A. Demircan, S. Amiranashvili, C. Brée, C. Mahnke, F. Mitschke, and G. Steinmeyer, “Rogue wave formation by accelerated solitons at an optical event horizon,” Appl. Phys. B 115, 343–354 (2014). [CrossRef]  

41. G. Genty, C. de Sterke, O. Bang, F. Dias, N. Akhmediev, and J. Dudley, “Collisions and turbulence in optical rogue wave formation,” Phys. Lett. A 374, 989–996 (2010). [CrossRef]  

42. V. V. Voronovich, V. I. Shrira, and G. Thomas, “Can bottom friction suppress ‘freak wave’ formation?” J. Flu. Mech. 604, 263–296 (2008). [CrossRef]  

43. V. Zakharov, A. Pushkarev, V. Shvets, and V. Yan’kov, “Soliton turbulence,” Pis’ma v Zhurnal Eksperimental’noi i Teoreticheskoi Fiziki 48, 79–82 (1988).

44. U. Bortolozzo, J. Laurie, S. Nazarenko, and S. Residori, “Optical wave turbulence and the condensation of light,” J. Opt. Soc. Am. B 26, 2280 (2009). [CrossRef]  

45. A. Picozzi, J. Garnier, T. Hansson, P. Suret, S. Randoux, G. Millot, and D. Christodoulides, “Optical wave turbulence,” Phys. Rep. 542, 1–132 (2014). [CrossRef]  

46. M. Eberhard, “Massively parallel optical communication system simulator,” https://github.com/Marc-Eberhard/MPOCSS.

47. G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, 2013).

48. V. E. Zakharov and A. B. Shabat, “Exact theory of two-dimensional self-focusing and one dimensional self-modulation of waves in nonlinear media,” Sov. Phys. JETP 34, 62–69 (1972).

49. V. E. Zakharov, “Turbulence in Integrable Systems,” Stud. Appl. Math. 122, 219–234 (2009). [CrossRef]  

50. M. Taki, A. Mussot, A. Kudlinski, E. Louvergneaux, M. Kolobov, and M. Douay, “Third-order dispersion for generating optical rogue solitons,” Phys. Lett. A 374, 691–695 (2010). [CrossRef]  

51. N. Akhmediev and M. Karlsson, “Cherenkov radiation emitted by solitons in optical fibers,” Phys. Rev. A 51, 2602–2607 (1995). [CrossRef]   [PubMed]  

52. K. Tai, A. Hasegawa, and A. Tomita, “Observation of modulational instability in optical fibers,” Phys. Rev. Lett. 56, 135–138 (1986). [CrossRef]   [PubMed]  

53. A. Hasegawa, “Generation of a train of soliton pulses by induced modulational instability in optical fibers,” Opt. Lett. 9, 288 (1984). [CrossRef]   [PubMed]  

54. V. Karpman and V. Solov’ev, “A perturbational approach to the two-soliton systems,” Phys. D: Nonlin. Phenom. 3, 487–502 (1981). [CrossRef]  

55. A. V. Buryak and N. N. Akhmediev, “Internal friction between solitons in near-integrable systems,” Phys. Rev. E 50, 3126–3133 (1994). [CrossRef]  

56. H. Degueldre, J. M. Jakob, T. Geisel, and R. Fleischmann, “Random focusing of tsunami waves,” Nat. Phys. 12, 259–262 (2016).

57. A. Slunyaev and E. Pelinovsky, “Role of Multiple Soliton Interactions in the Generation of Rogue Waves: The Modified Korteweg-de Vries Framework,” Phys. Rev. Lett. 117, 214501 (2016). [CrossRef]   [PubMed]  

58. Y.-H. Sun, “Soliton synchronization in the focusing nonlinear Schrödinger equation,” Phys. Rev. E 93, 052222 (2016). [CrossRef]  

59. M. Eberhard, A. Savojardo, A. Maruta, and R. A. Römer, “Rogue wave generation by inelastic quasi-soliton collisions in optical fibres,” http://dx.doi.org/10.17036/5067df69-bf43-4bea-a5ff-f0792001d573 (2017).

60. W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes in C, 2nd ed. (Cambridge University, 1992).

61. Y. Nagashima, Elementary Particle Physics, vol. 1, (Wiley-Vch, (2010). [CrossRef]  

62. Z. Chen, A. J. Taylor, and A. Efimov, “Soliton dynamics in non-uniform fiber tapers: analytical description through an improved moment method,” J. Opt. Soc. Am. B 27, 1022 (2010). [CrossRef]  

63. J. Santhanam and G. P. Agrawal, “Raman-induced spectral shifts in optical fibers: general theory based on the moment method,” Opt. Commun. 222, 413–420 (2003). [CrossRef]  

64. A. Savojardo, “Rare events in optical fibers”, PhD Thesis (2017).