On different aspects of the optical rogue waves nature

Rogue waves are giant nonlinear waves that suddenly appear and disappear in oceans and optics. We discuss the facts and fictions related to their strange nature, dynamic generation, ingrained instability, and potential applications. We present rogue wave solutions to the standard cubic nonlinear Schrödinger equation that models many propagation phenomena in nonlinear optics. We propose the method of mode pruning for suppressing the modulation instability of rogue waves. We demonstrate how to produce stable Talbot carpets—recurrent images of light and plasma waves—by rogue waves, for possible use in nanolithography. We point to instances when rogue waves appear as numerical artefacts, due to an inadequate numerical treatment of modulation instability and homoclinic chaos of rogue waves. Finally, we display how statistical analysis based on different numerical procedures can lead to misleading conclusions on the nature of rogue waves.


Introduction
Analytical and numerical solutions of the nonlinear Schrödinger equation (NLSE) of different orders have been widely analyzed for their importance in a number of mathematical and physical systems [1][2][3][4][5][6][7][8][9][10][11]. In this work, we focus on the simple cubic one-dimensional NLSE arising in many fields of nonlinear and fiber optics. Our main concern are the NLSE solutions that are unstable due to the influence of modulation instability (MI). Although MI is usually recognized as the Benjamin-Feir instability, it already appeared in 1947, in the work of Bogoliubov on the uniform Bose gas [12], which introduced and discussed the MI of the cubic NLSE. In general, the modulation instability can be regarded as a nonlinear optical process where the power of the fundamental pump wave is attenuated and redistributed to a finite number of spectral sidebands. These higher-order modes are very weak at the onset of nonlinear evolution but their power increases exponentially during propagation [3,4,13]. The notion of MI is widely spread over several fields of physics. The NLSEbased study of nonlinear modulated waves and modulation instability in real electrical lattices was presented in [14]. The properties of breathers and modulation instability in a discrete nonlinear lattice are investigated analytically in [15]. The issue of MI was also investigated in electronic wave packets in nonlinear chains modeled by the discrete NLSE [16].
Since MI is frequently indicated as the main cause of the rogue wave (RW) appearance in nonlinear optics [17], we are further motivated to investigate how this fundamental property of nonlinear systems affects various high-intensity NLSE solutions that can be regarded as rogue waves. This research is becoming more exciting since a new scheme for the excitation of rogue waves, via the electromagnetically induced transparency (EIT) [18][19][20][21], has been analyzed in [22].
In order to properly investigate the impact of MI on RWs, one needs to dynamically generate high-intensity peaks with narrow spatial profiles. We start from the well-known exact solution of the cubic NLSE, known as the higher-order Akhmediev breather (AB) [23][24][25][26][27]. A complication arises along the way, as these higherorder solutions also arise as the homoclinic orbits of unstable Stokes waves during the evolution of the system under NLSE [28][29][30][31]. Consequently, the long-time dynamics will become chaotic for ABs with two or more unstable modes. It is then questionable whether the chaotic behavior is intrinsic to the model equation or is induced by the numerical algorithm applied [28,29].
To address this question, we revisit our previous results [32] on the analytical and numerical NLSE solutions that are periodic both along the temporal and spatial axes, known as the Talbot self-images or carpets [25,33,34]. In optics, the Talbot effect is known as a near-field diffraction phenomenon, occurring when light beams undergo diffraction at some periodic structure and produce recurrent self-images at equidistant planes. The Talbot carpets refer to fractional and even fractal light patterns observed in-between the planes. Although Talbot effect is known for more than hundred years, it still finds applications in many fields, such as plasmonic nanolithography [35][36][37].
In reference [32], we presented a mode pruning procedure to mitigate the unavoidable impact of MI and obtain double-periodic solutions from Akhmediev breathers. Here, we briefly recap this effort: We calculate ABs of different orders using Darboux transformation (DT) technique and extract initial conditions in a wide box that is a multiple of the main breather period [38]. We next apply the mode pruning method combined with a chosen numerical NLSE integration scheme to preserve the self-imaging recurrences for as long as possible. In this paper, we extend this analysis by taking into account different values of the main breather parameter a (Sect. 2.1) and initially having more breathers in the box. But, the central theme of this paper is concerned with the conflicting opinions formed about the nature of RWs: Are they linear or nonlinear; random or deterministic; numerical or physical? A short and in our opinion correct answer to these nagging questions is as follows.
Rogue waves are essentially nonlinear, because their cause is the modulation or Benjamin-Feir instability, which is one of the basic nonlinear optical processes. They are deterministic, because modulation instability leads to homoclinic chaos, which by its nature is deterministic; random phenomena are probabilistic and may look chaotic but are not deterministic. RWs are physical, because they are observed in many experiments and media, with similar statistics [39][40][41].
Now, there exist reservations to these facts, depending on how one generates and analyzes RWs. This especially holds if the work is numerical and statistical in nature. Nevertheless, even in this case, we will display how widely used beam propagation method can produce spurious RWs when inappropriately employed. We will also demonstrate how statistical analyses of exactly the same dynamical systems, but produced using different numerical algorithms, may lead to different statistics of RWs. We will address these reservations later in the paper.

Rogue wave solutions to the NLSE
To recap, we discuss the nature of optical rogue waves in the cubic NLSE, in view of conflicting opinions expressed in the literature. In particular, as already mentioned, we address three pairs of opposing suppositions on their nature: Linear versus nonlinear [3,5]; random versus deterministic [42,43]; and numerical versus physical [28,29]. A short answer to these suppositions is that rogue waves in optics are nonlinear, deterministic, and physical. To stress again, they are nonlinear because the major cause of rogue waves is the modulation or Benjamin-Feir instability. They are deterministic because modulation instability leads to deterministic chaos. They are physical because they appear in many experiments and media. Our opinion is supported by extensive numerical simulations of the nonlinear Schrödinger equation in different regimes that touch upon the aspects of all three conflicting suppositions.
Disturbingly, in numerical simulations optical rogue waves may appear fictitiously, as numerical artefacts. One such published instance will be displayed below. Different numerical algorithms for exactly the same input may provide different evolution pictures and different statistics when modulation instability sets in and when numerical grid parameters are not set precisely. An example of different statistics of the probability density distributions, obtained by two different numerical algorithms applied to the same input, is also provided below. There, the standard beam propagation method (the split-step fast Fourier transform method), predominantly used in the literature, predicts the appearance of thousands of rogue waves, whereas the more precise symplectic algorithm of higher order predicts significantly fewer. Hence, owing to the vague definition of rogue waves, the way they are produced, and an exponential amplification of numerical errors, there are situations in which optical rogue waves may appear as linear, random, and numerical.

The influence and suppression of modulation instability
Before proceeding to these more subtle points, let us start by discussing the more standard behavior of RWs in the standard NLSE. We study double-periodic solutions of the simplest cubic nonlinear Schrödinger equation The wave function ψ ≡ ψ (x, t) corresponds to the slowly varying envelope that could be optical, plasmonic or other in nature. In a fiber optics notation, the transverse variable t is the retarded time in the frame that moves with a pulse group velocity. The evaluation variable x is considered as the distance along the fiber. Let us consider the first-order AB solution, which is single-periodic along t-axis [24,25]: (2. 2) The period L, angular frequency ω, and growth rate λ of an AB (first-or higher-order) are solely determined by parameter a, with 0 < a < 0.5 [13]: (2.5) When a = 0.5, AB turns into the Peregrine RW, and for a > 0.5, it becomes the Kuznetsov-Ma (KM) soliton. The procedure for dynamic generation of the first or higher-order Akhmediev breathers is explained in detail in Section 2 of [32]. Here, we briefly state that the initial condition for numerics is derived from an exact AB solution at a certain value x 0 of the evolution variable, using Darboux transformation [44]. To obtain double-periodic breathers resembling the Talbot carpets solutions, we need to adjust the size of the transverse interval to an integer multiple of the fundamental breather's period L, apply periodic boundary conditions and use the mode pruning method. Namely, when the box size is equal to the breather's fundamental period, the Fourier harmonics will form the set S 1 of N stable fundamental breather modes (N denotes the number of grid points along the transversal t-axis). If we extend the box to an integer multiple of the main period, the mode spacing will get smaller and a new set S 2 will be obtained. The way to generate nonlinear Talbot carpets is to suppress modulation instability caused by the exponential growth of unstable Fourier modes belonging to S 2 and not to S 1 .
In Fig. 1a, we present the numerical evolution of the first-order Akhmediev breather (a = 0.33) when the box size is equal to five breather's periods. The intensity peaks at x = 0 are repeated along x-axis at the Talbot periods. This array is shifted for half-a-period along the t-axis during evolution, forming the secondary Talbot image at half the Talbot period. The corresponding Fourier spectrum is shown in Fig. 1b. This result is obtained when unstable subharmonics are eliminated after each iteration, leaving only those that form the fundamental AB, with indices 0, ±5, ±10, ±15, etc. This is in agreement with our previous work [32] showing the same behavior for a = 0.36. Our motivation to study this problem for NLSE also originates from the similarity of our results shown in Fig. 1 to those presented in [35] (their Fig. 2) and [36] (their Figs. 2, 3, and 4) which are related to the Talbot effect in lithography. Next, we want to verify these results for even larger period, having seven breathers in the box. We apply the simple pruning algorithm to Fourier modes, setting all unstable mode amplitudes to zero except the modes indexed 0, ±7, ±14, ±21, and so on. The result is an extended Talbot carpet with alternate shifting of intensity maxima along the x-and t-axes, as shown in Fig. 1c. The Fourier spectrum of the mode amplitudes is shown in Fig. 1d. If the pruning algorithm is not applied, the chaotic behavior ruins the carpet after three Talbot cycles, as shown in Fig. 1e. The unstable modes having the index not divisible by 7 grow exponentially due to MI, and prevent the homoclinic orbit from returning to itself after more than three cycles. The other view on the Talbot carpet disintegration is the irregular buildup of the Fourier modes spectrum, shown in Fig. 1f. As a result of the interaction or beating of different AB modes in the chaotic region, four second-order ABs are formed around x = 41, x = 61, and x = 67. We may consider these maxima as channels for the RW production through MI.
Another procedure to curb MI is the Gaussian pruning algorithm [32], in which the unstable modes are not eliminated completely but suppressed by a Gaussian factor. Thus, the unstable modes are multiplied by a number that depends on their strength: if the modes grow more, then the level of suppression is higher. In turn, the modes can grow only up to a certain value, determined by the Gaussian distribution, and cannot exceed the predefined value in the range of 10 −11 to 10 −9 (not shown).
To finalize, in Fig. 2, we present the nonlinear Talbot carpet consisting of the second-order breathers with higher intensities and narrower profiles. The box size is seven times the fundamental breather period, having a = 0.41. This is an additional analysis of our previous example in [32], but here with more breathers in the box and decreased spatial resolution. Figure 2a represents an unstable situation without pruning, while Fig. 2b is the stabilized carpet. The simple pruning algorithm was used, which left only each seventh Fourier mode intact.
In the end, one may wonder what such pruning procedures are worth for and how they can be utilized in the real world? After all, suppressing unstable Fourier modes numerically is a relatively simple task, but what its relevance might be in actual experiments is another matter. However, it is not difficult to envision situations in which it might be useful: for example, in experiments where signal transmission is interrupted periodically (and electronically) for data analysis. This typically happens when you have an experiment in an optical loop, where the beam circulates around. It also happens in long fiber experiments, when periodically you have to interrupt the propagation and ensure that it is still stable. The question for experimenters is then, how viable or advisable is to prune the data at the interruption, and inject back into the experiment, to force the beam and control the instability?

The influence of numerics on the dynamics and statistics of rogue waves
In this subsection, we address the questions pertaining to the numerics and statistics of RWs. But, before going to the details of numerics and statistics, let us briefly recap how we got to this point. To recall, the generating mechanism of optical rogue waves is the Benjamin-Feir or modulation instability. It is the basic nonlinear optical process in which a weak perturbation of the background pump wave produces an exponential growth of spectral sidebands that constructively interfere to build RWs. We have produced RWs in numerical simulations of the cubic nonlinear Schrödinger equation with inputs on the flat background determined by the DT. Optical RWs represent homoclinic orbits of unstable sideband modes that, due to MI, generate homoclinic chaos. The question is then whether the chaos seen belongs to the model itself or is induced by the numerical procedure applied. After all, even though the model is the same, different numerical algorithms for its solution represent different dynamical systems and may approach chaos differently. As it will be discussed below, different numerical algorithms lead to similar values of the wave function once the evolution step in numerics dx for each algorithm is carefully chosen, but only up to a certain value x cr . If the evolution is continued beyond x cr , different dynamics will be generated using different numerical schemes, due to combined effects of finite numerical precision and modulation instability. It is common occurrence that in the region of instability and chaos different numerical methods offer different solutions to the same partial differential equation, under the same initial and boundary conditions. Sensitivity to small differences in initial conditions is one of the hallmarks of deterministic chaos. We present such an occurrence by simulating the same NLSE, with exactly the same parameters, inputs and boundary conditions, but using different numerical algorithms. We compare the standard second-order beam propagation method (BPM) with various higherorder symplectic algorithms. In our work, we prefer symplectic integrators for their utility in treating chaotic Hamiltonian systems that require long-time integration of noisy inputs. We demonstrate that distressingly, different algorithms may provide different evolutions and different statistics of the RW peaks, owing to homoclinic chaos and imprecise choice of the evolution step dx for each numerical procedure. Even the same algorithms may provide different evolutions when the numerical step size is changed.
An appropriate example is the appearance of the third-order RW solution as a result of the collision of three ABs [45] for an initial condition given by ψ (t, x = 0) = e iπ t/10 [1 + 0.002 cos (π t/10)]. It is known that the collision of three localized pulses at the same point is unlikely to happen and difficult to observe, so we were intrigued by this finding in an otherwise excellent paper. Hence, we reproduced the same results by taking dx = 10 −4 (Fig. 3a). However, we discovered that just by taking half of the original numerical step size in the same numerical method, the third-order RW disappears (Fig. 3b). To check the result, we recreated the solution using a fourth-order symplectic algorithm, which confirmed the disappearance of the RW (dx = 10 −4 , Fig. 3c). It appears that the third-order RW in this example is just a numerical artefact.
When one solves the same equation by two numerical methods and finds different answers, an obvious question arises: Which of the two (if any) is correct? The question is especially relevant when one deals with the RWs in homoclinic chaos, because in chaos sooner or later all algorithms give different answers. Hence, the appropriate question to ask is not which solution is correct (because in chaos there exist no correct-meaning precise-solutions) but which algorithm is "more" correct. Then, the answer is obvious-the one that keeps to the presumed correct solution for longer (or conversely, the one shadowed by the correct solution for longer). Let us clarify this by presenting an example.
In Fig. 4, we analytically generate the fifth-order Akhmediev breather, represented by a very sharp and narrow peak at x = 13, using Darboux transformation. The main breather's frequency ω 1 is determined by parameter a = 0.4850173 (Eq. (2.4)). The remaining four higher-order modes (i = 2, 3, 4, 5) have harmonics with higher frequencies ω i = iω 1 . We take the initial condition from DT at x = 0 (well before the peak) and numerically evolve the wave function, in order to see whether RW will appear at x = 13 and what happens after the peak. In the figure, we compare four different algorithms, by solving the cubic NLSE for the same initial condition around a highly unstable homoclinic orbit-the fifth-order RW. Such a RW, towering 81-times over the background intensity, is never seen in experiments and rarely even in numerical simulations.
One can see that all numerical algorithms correctly resolve the appearance of the RW, but afterward they all sooner or later give different answers. This happens because of the developing MI, which entails different evolutions in different algorithms. However, one can notice that the fourth-order and the sixth-order symplectic algorithms provide visually consistent solutions over the whole integration interval of 100 propagation units (10 5 integration steps!), while others do not. One may infer that these two numerical pseudo-orbits are shadowed by the true orbit over a longer distance. For this reason, in simulations that require long-distance integrations of noisy inputs that are necessary for the formation of statistics of RWs, we employ the fourthorder symplectic algorithm.
Here, it is worth mentioning that almost identical wave functions are observed up to some evolution point when running three different algorithms starting from the same initial condition, but with different integration steps: BPM for dx = 0.00001, sixth-order symplectic for dx = 0.0002, and eighth-order symplectic for dx = 0.0005. Quite expectedly, the lower-order algorithms require smaller steps for similar accuracy, whereas higher-order allow for larger steps. Nevertheless, as the evolution x-coordinate advances, the differences emerge due to both finite numerical precision and the intrinsic MI (not shown).
The fact that individual RWs may appear in different algorithms at different places and thus may represent fictitious structures or numerical artefacts is not the worst feature of dealing with RWs in homoclinic chaos. The more troublesome is the fact that the statistics of RWs obtained by different numerical methods may also be different, because by fiat, these statistics come from performing long-distance evolutions of a solution (or a class of solutions) to the NLSE, which are different for different numerical methods. It is desirable that they all provide the same statistics for the same underlying model, but this is not written in stone. As already mentioned, different numerical schemes represent different dynamical models. Worst of all, the long tail statistics of high intensity waves, i.e., the RWs, seems to be most affected by different numerical methods. Again, this is best explained by an example.
An example of two different evolution scenarios, starting from white noise with amplitude A = 5% around the background intensity B = 1, is shown in Fig. 5. We use two different methods: The second-order BPM (top) and the fourth-order symplectic (bottom). We intentionally use the same evolution step of dx = 10 −3 and the same transverse box (−69.42, +69.42) (divided into N = 2048 intervals of the same width) for both figures, in order to spot the differences in the intensity graphs for identical initial conditions and numerical grids; the only difference is different algorithms. The evolution of intensity by the two algorithms appears consistent to about 100 longitudinal units (each unit is covered by 1000 numerical iterations), but after that the distributions become different. The wave function evolution is observed up to x = 6500, but for clarity the intensity is shown only from x = 0 to x = 250 in both figures. If certain regions of intensity maps are To settle a fine but important point that will become apparent later, we examine numerical evolution of the wave function ψ(x, t), using the fourth-order symplectic algorithm, when the initial condition is chosen so that ψ is normalized at x = 0. The module of the wave function is The last integral is approximated by numerical formula where L indicates the integration interval (box size), N is the number of points along the transverse t-axis, ψ i ≡ ψ (x = 0, t i ), and Δt i = Δt = L/N . The initial wave function ψ B=1, A=5% in Fig. 5 is not normalized and we can chose two different ways to set its module to one. In the first case, one can divide all N initial complex numbers with so that ψ norm1 = ψ B=1, A=5% / √ C. In the second case, one can hold these numbers unchanged and modify the grid size L so that the module is set to one according to Eq. (2.7): ψ norm2 = ψ B=1, A=5% and L = L/C. In the first case of lower intensity initial condition, the numerical evolution is very slow and the first structures that resemble Akhmediev breathers appear at x ≈ 800 (Fig.  6a). The most similar structure appears at t ≈ −47 and x ≈ 2450, but evidently it is hugely extended along the evolution x-axis, with an approximate width of Δx ≈ 300. This is far more stretched compared with the similar features obtained for non-normalized initial conditions, shown in Fig. 5b. In the second case, when the numerical t−box is significantly decreased (from L ≈ 140 to L ≈ 1) the evolution gave us no char- acteristic NLSE structures (Fig. 6b). One can easily spot the uniform intensity over the entire xt-grid with no vertical traces (Kuznetsov-Ma solitons) or localized maximums (Akhmediev breathers of the first or second order). Numerically, this is ill-advised.
What one can learn from these figures is that the ψ evolution is determined by the box size and the intensity of random wave function at the beginning of numerical iterations. Namely, the intrinsic modulation instability is the phenomenon that forms the characteristic NLSE structures during evolution. The amplitudes of higher Fourier modes that are very weak at x = 0 grow exponentially and produce KM solitons and ABs. If there are more unstable modes at beginning, then MI and NLSE features will appear sooner. The number of unstable modes N um is related to the background noise level B and box size L via the following expression: (2.8) In case of low B and large L (see Fig. 6a, B = 0.072 and L ≈ 140) we get N um = 3.2. This means we have only a few unstable modes and MI is weak. That is the reason for obtaining very slow evolution and stretched NLSE features. In the second case, we have B = 1 and the small box L ≈ 1 and compute N um = 0. 16. This means that all Fourier modes are stable and the evolution is uniform without solitons and breathers formation, in agreement with Fig. 6b. We conclude that the normalization requirement imposes unacceptable conditions for noise background level and box size that adversely affect the dynamics of white noise under the NLSE model. Finally, we run the second-order BPM algorithm on the same initial condition and got the same results (not shown). The overall conclusion is that the integration with normalized ψ should be avoided, as it leads to different model behavior that is difficult to efficiently integrate. We now go back to the case of non-normalized initial conditions (B = 1 and different values of A). We analyze the statistics of intensity maxima for evolutions obtained using the second-order BPM and the fourth-order symplectic algorithms. During numerical integration, our software records intensity maxima I i, j = |ψ i, j | 2 at (x i , t j ) point of the numerical grid, if this intensity is the highest value in the 8-connected neighborhood region. All maxima are being recorded in a file and then plotted in the histogram form, where each bin height corresponds to the number of maxima within the bin's intensity interval, divided by the total number of peaks. The additional requirement after each computational iteration is to eliminate the high-frequency Fourier modes forming the fast zig-zag glitches in the wave function ψ that are considered as false maxima. We should note that although our computational procedure appears to be rather stringent, it affects the chosen algorithms equally, so that the resulting differences are not caused by this feature.
First, we discuss the probability density of peaks' intensities from Fig. 5, where the white noise amplitude is A = 5%. The results for BPM are presented in Fig. 7a, and for the fourth-order symplectic algorithm in Fig. 7b. It turns out that the results in general are different. Both distributions follow the familiar exponential decay of MI-driven systems, which may not be the case for the density of high-intensity peaks, forming RWs. At this point, we pose the question of what is the criterion for announcing an intensity peak as a rogue wave? Here we simultaneously adopt two criteria and make comparative analysis. The first is relatively simple: the peak with the intensity greater than that of the Peregrine soliton (I ≥ I PS = 9) is considered as a rogue wave. The second is statistical in nature and is widely adopted in the literature [46][47][48][49][50]: the intensity threshold I RW above which one obtains a RW is a mean value of the largest third of intensity peaks multiplied by two (the same definition applies to the rogue wave amplitude in oceanography). Therefore, in Figs. 7 and 8, we indicate two vertical lines, dividing the intensity scale on RW and non-RW regions.
In Figs (with c = 0.34151, I 0 = 1.25399), respectively. Thus, while the distributions are well-described by exponentially decaying curves, in detail they differ most in the heavy-tailed region of higher intensities, where RWs reside.
The next relevant question is what is the influence of the white noise amplitude A (at x = 0) on the statistics of NLSE evolution? We therefore run six simulations using two algorithms (BPM and fourth-order symplectic) for three values of A. The corresponding histograms of intensity maxima are shown in Fig. 8a (BPM alg. and A = 1%), 8b (4S alg. and A = 1%), 8c (BPM alg. and A = 2%), 8d (4S alg. and A = 2%), 8e (BPM alg. and A = 3.5%), and 8f (4S alg. and A = 3.5%). Here, as in Fig. 5, we removed fast glitches from the wave function that present false maxima, by eliminating high-frequency Fourier modes. The statistics are similar to those obtained with 5% amplitudes, but the differences in the results between the two algorithms for the same value of A still persist.
All data for two algorithms and different amplitudes are summarized in Table 1. The first observation is that for the background level of B = 1, the statistical RW threshold I RW is lower than I PS = 9. Namely, the majority of intensity peaks are located at lower intensities for all eight histograms. Thus, the largest third of the peaks begins well below the Peregrine soliton intensity. Another conclusion is that the larger the amplitude A, the higher the statistical threshold I RW . As for the overall intensity maximum I max during numerical evolution, we observe that it depends on the algorithm and grid parameters.
We have also increased the value of the background level (B = 1.5) and realized that both I RW and I max increase significantly, so that I RW > 9. However, once we magnify the breather structures in such an evolution, we realize that very high intensities at the breather's center do not correspond to the n th -order AB, as expected (n ≥ 3), but retain the intensity distribution pattern characteristic of the second-order AB. Such dynamics can deceive the observer since it produces Akhmediev breathers of intrinsically larger intensity that could be generated through Darboux transformation using the seed function with modulus greater than one (which usually is not the case in literature). Therefore, one should stick to the B = 1 case.
The third conclusion is that the higher the noise amplitude A is, the larger is the total number of maxima N TP . In Fig. 8, we also observe millions of peaks and thousands of RWs (N RW and N PS , according to the adopted RW definition). The fractions of peaks that can be considered as rogue waves (F RW and F PS ) is considerably lower than 1% and relatively weekly depend on the algorithm and noise amplitude A.
Our final conclusion is that the most prominent difference between statistical results obtained by two algorithms for the same A value is observed at higher intensities. The second-order breathers that approach the limiting intensity value of 25 at its center arise in the collision of highly unstable modes during evolution. By switching from BPM to 4S algorithm one changes the way and the number of mathematical operations needed for the calculation of ψ values. This difference affects more the collision processes of highly unstable modes that then lead to high intensity maxima. Overall, these conclusions strengthen the claim that different numerical algorithms produce different statistics, especially in the long tail region of high intensities, where RWs exist.
We got similar results and confirmed our conclusions by using two times smaller evolution step (dx = 0.0005) or the transversal box size (L/2 ≈ 35, instead of 70) (results not shown).

Conclusion
In conclusion, we have discussed the facts and fictions related to the strange nature, dynamic generation, ingrained instability, and potential applications of RWs. We have proposed the method of mode pruning for suppressing the modulation instability of rogue waves. Using this procedure, we have computed and demonstrated the stable numerical Talbot carpets-the recurrent images of light and plasma waves-by rogue waves, for possible use in nanolithography. The pruning procedure was found indispensable in the production of stable recurring periodic images over wide but still finite windows.
We have also discussed the nature of optical rogue waves, in view of conflicting opinions expressed in the literature. In particular, we have addressed the three pairs of opposing suppositions on their nature: Linear versus nonlinear; random versus deterministic; and numerical versus physical. In summary, a correct answer to the three suppositions is that the rogue waves in optics are essentially nonlinear, deterministic, and physical. They are nonlinear because the major cause of rogue waves is the modulation or Benjamin-Feir instability, which by its nature is the basic nonlinear optical process. Rogue waves are deterministic because modulation instability leads to deterministic chaos; random phenomena are probabilistic and may look chaotic but are not deterministic. Rogue waves are physical because they appear in many experiments and media, with similar statistics.
Nevertheless, in numerical simulations optical rogue waves may appear fictitiously, as numerical artefacts. Different numerical algorithms for exactly the same inputs may provide different evolution pictures anddistressingly-different statistics of RWs, caused by imprecisely chosen integration steps and intrinsic modulation instability. An overall conclusion is that the numerical and statistical treatments of NLSE by different algorithms represent different dynamical systems that may introduce different long-time behaviors when they are unstable and under the influence of MI. Therefore, owing to a vague definition of rogue waves, different ways in which they are generated, and an exponential amplification of numerical errors in chaos, there are situations in which optical rogue waves may appear as linear, random, and numerical. Data Availability Statement All data generated or analyzed during this study are included in the published article.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.