Analysis of interaction dynamics and rogue wave localization in modulation instability using data-driven dominant balance

We analyze the dynamics of modulation instability in optical fiber (or any other nonlinear Schrödinger equation system) using the machine-learning technique of data-driven dominant balance. We aim to automate the identification of which particular physical processes drive propagation in different regimes, a task usually performed using intuition and comparison with asymptotic limits. We first apply the method to interpret known analytic results describing Akhmediev breather, Kuznetsov-Ma, and Peregrine soliton (rogue wave) structures, and show how we can automatically distinguish regions of dominant nonlinear propagation from regions where nonlinearity and dispersion combine to drive the observed spatio-temporal localization. Using numerical simulations, we then apply the technique to the more complex case of noise-driven spontaneous modulation instability, and show that we can readily isolate different regimes of dominant physical interactions, even within the dynamics of chaotic propagation.


INTRODUCTION
Modulation instability (MI) of the nonlinear Schrödinger equation (NLSE) describes the process whereby a weak perturbation experiences exponential growth at the expense of a strong input wave [1,2].MI (sometimes called the Benjamin-Feir or Bespalov-Talanov instability [3,4]) leads to complex spatio-temporal pattern formation, and is one of the fundamental nonlinear dynamical processes of nature.It has been observed in many different systems including hydrodynamics, plasmas, Bose-Einstein condensates, and fiber optics.Despite this large body of work over many years, its centrality to nonlinear science is such that it continues to be extensively studied from both experimental and theoretical perspectives.Recent work, for example, has explored its description in terms of integrable turbulence [5,6], its relationship with computational complexity [7], its thermodynamic link to the soliton-gas concept [8], and its intrinsic association with Fermi-Pasta-Ulam-Tsingou recurrence [9]; to cite only a small number of examples.
The dynamics of MI leads to the spontaneous emergence of localized structures that possess different spatial and/or temporal periodicities [10,11].These structures are intimately connected with known analytic solutions to the NLSE (including Peregrine soliton and Akhmediev and Kuznetsov-Ma breathers [12][13][14]), and understanding this correspondence has allowed experimentalists to excite a wide range of soliton and breather solutions in both optics and hydrodynamics [15][16][17][18].
Moreover, even under conditions where modulation instability is excited from noise, it has been shown that the random peaks developing from noise possess the expected characteristics of these analytic solutions [19][20][21].It is these nonlinear localization dynamics in particular that have attracted great interest as potentially underpinning the growth and decay of destructive rogue waves on the ocean [22,23].
These various studies have yielded significant insights into the properties of MI under diverse conditions, and in different systems.Somewhat surprisingly, however, although some aspects of MI localization can be interpreted precisely using mathematical methods such as the inverse scattering transform [24], the physics of nonlinear and dispersive interactions in MI is more often discussed in qualitative terms by a comparison with specific limiting cases or characteristic nonlinear and dispersive length scales [25].It would be highly desirable to have a means of interpreting the physics of MI that went beyond such a qualitative description, and yet avoided the formalism of the inverse scattering method.
In this paper, we show that the machine-learning technique of data-driven dominant balance can address this problem.Machine learning methods are currently of great interest in all areas of physics [26][27][28], and in the particular field of nonlinear optics, have been applied to the study of various NLSE propagation scenarios [29][30][31][32].The technique of dominant balance aims to automatically determine the contributing dominant physical processes at each step of propagation.As a subset of unsupervised learning techniques, it has been successfully applied to interpret the physics of a number of nonlinear propagation scenarios in hydrodynamics, as well as the more challenging case of broadband supercontinuum generation [33].
In this paper, we use a dominant balance approach to analyse modulation instability of the NLSE.We first apply the method to interpret known analytic solutions for Akhmediev breather, Kuznetsov-Ma, and Peregrine soliton structures, and for these spatio-temporal dynamics, we show how we can distinguish background regions of dominant nonlinear propagation from regions where nonlinearity and dispersion interact to drive localization.This is especially important in showing how dominant balance can provide complementary insights into the dynamics, because associating the nonlinear stage of evolution with the background may seem counter-intuitive as this is a region of low intensity.Following these studies of analytic SFB solutions we then use numerical simulations to study the more complex propagation case of noise-driven chaotic MI, and find again that we can automatically identify these different regimes of physical interaction.

NLSE SOLUTIONS
We consider MI occurring in the focusing NLSE which is written in normalised form as follows: Here ψ(ξ, τ ) is a field envelope evolving in distance ξ and co-moving time τ .Dimensionless variables ξ and τ are related to the usual notation of nonlinear optics by ξ = z/L NL and where L NL = (γP 0 ) −1 .Here z and t are dimensional distance and time, P 0 is power (usually that of the input continuous wave), and β 2 and γ are the usual dimensional fiber group velocity dispersion and nonlinearity parameters respectively [25].The field envelope ψ(ξ, τ ) is normalized with respect to P 1/2 0 .
The NLSE possesses a number of known analytic solutions [11,34].Those associated with MI are the solitons on finite background (SFB), that can be written in compact form as follows: The physical behaviour of the solution is determined by the single governing parameter a through , ω m = b = 0 and the solution is the limiting rational Peregrine soliton, double-localized in ξ and τ [14].For a < 1/2, ω m and b are real, and we obtain the τ -periodic Akhmediev breather, with ω m and b taking on physical significance of a modulation frequency and exponential growth/decay rate respectively.
When a > 1/2, ω m and b become imaginary, and we obtain the ξ-periodic Kuznetsov-Ma solution.
These various SFB structures are well known, and have been observed in a range of experiments since 2010 [35][36][37].

IMPLEMENTING THE DOMINANT BALANCE TECHNIQUE
In this section, we give a general overview of how the dominant balance technique and algorithm are applied to nonlinear propagation in the NLSE.Further details and references are given in the Methods section.The dominant balance technique aims to automate the process of identifying the key interacting physical processes associated with different spatio-temporal regions of evolution.
The technique involves several steps.The first is to determine the evolution of the field ψ(ξ, τ ), and this is straightforward here as we have access to the analytic result in Eq. ( 2).However, as we see below for noise-driven MI, the evolution can also be obtained using numerical integration of the NLSE.Indeed, in the most general case, this could also involve analysis of experimental data when access to full field information is available [38,39].
The second step analyses the evolution ψ(ξ, τ ) in its associated "equation space," where each coordinate axis corresponds to a physical process defined by one of the terms in the governing NLSE (see Methods).Specifically, for each point (ξ, τ ), the NLSE terms {iψ ξ , ψ τ τ , ψ|ψ| 2 } are separately computed, and we search for a "dominant balance" regime where the NLSE is approximately satisfied by only a subset of terms (the other terms contributing only negligibly.)As shown in Ref. [33], machine learning tools can automate this search, using cluster detection (Gaussian mixture modelling) and sparse regularization to identify regions where different combinations of terms drive the dynamics.These are standard tools of unsupervised learning and optimization, and allow robust detection of clusters even when they overlap (see Methods) [28,40].When different clusters are found to possess the same sparcity pattern (significantly reduced variance in the same directions of equation space), these are grouped together to form a particular candidate "balance model."In the case of the NLSE with three possible interacting terms, this process has a simple geometric interpretation: two interacting terms will be associated with a cluster falling on a line in the three-dimensional equation space, three interacting terms will be associated with a cluster in a plane.
When the data is fully grouped into balance models, the final step is to re-map the clusters back onto the (ξ, τ ) space for comparison with the standard evolution dynamics.Visually, we do this by segmenting the original domain using a color key describing each balance model.In our analysis, we used the code package described in Ref. [33], and available at the online repository [41].We also note that since we are dealing with complex fields, we stacked real and imaginary components as input to allow grouping of regions of significant variance irrespective whether identified in the real or imaginary components [41].

RESULTS
We first apply this technique to identify locally-dominant interactions during the evolution of the three classes of SFB described above.Figure 1 shows results for the Peregrine soliton.
Specifically, Fig. 1(a-i The color-coded clusters are then mapped back onto a segmented dominant balance plot shown in Figure 1(a-ii), and the particular intensity profile at ξ = 0 is also plotted in Fig. 1(c-i) using the same color key.At ξ = 0, it is also instructive to plot the different contributions of each term of the equation space as shown in Fig. 1(c-ii), clearly revealing how different combinations of terms contribute to satisfy the NLSE (i.e.add to zero) in different regions.Note that at ξ = 0 all the three terms {iψ ξ , ψ τ τ , ψ|ψ| 2 } corresponding to the SFB solution are purely real.
These results reveal the key physical features of NLSE dynamics.For example, considering the Peregrine soliton and comparing Figs.1(a-i) and 1(a-ii), the orange region reveals how the strong spatio-temporal localization around (ξ = 0, τ = 0) arises from the interaction between all terms in the NLSE, as both nonlinearity and dispersion combine to drive spatio-temporal compression.
In contrast, the surrounding background region (blue) is dominated only by nonlinear evolution, and whilst this might be considered counter-intuitive since the background is where the intensity is lowest, this result actually highlights how interpreting NLSE physics requires comparison of the relative contributions of dispersion and nonlinearity.Specifically, a plane wave with no τ -structure can not "experience" dispersion, and thus it is only nonlinear self-focussing that can initially influence the evolution of the background.It is only after temporal structure develops from this nonlinear stage of evolution that dispersion and nonlinearity interact.In fact, this approach to visualizing the evolution very clearly illustrates the well-known "nonlinear" stage of the instability [34,42].The ability of the dominant balance analysis to identify this nonlinear stage explicitly (even though perhaps counter-intuitive from a naive perspective) is an example of how it can yield important insights into nonlinear evolution.The results in Figs 2 and 3 for the Akhmediev and Kuznetsov-Ma breathers respectively have similar interpretation.Here we see again see how regions of background associated only with dominant nonlinearity (blue) have been clearly identified, but we also clearly see how the contributions of all terms (orange) leads to the expected spatio-temporal localization characteristics.We also note how for the particular case of the Akhmediev breather, the ξ = 0 profile plot in Fig. 2(c) shows how all terms contribute to the dynamics in the lower amplitude regions between the localized peaks.These analytical SFB solutions, of course, do not exhaust the full variety of localised structures appearing in MI such as higher-order solutions [43], breather or soliton collisions [44], ghost interactions [45] etc.However, these key examples provide a clear indication of how the dominant balance approach can complement existing techniques such as inverse scattering transform [24,34,46] in interpreting NLSE dynamics.We now apply the dominant balance approach to interpret the more complex dynamics of noisedriven MI.For this case, the NLSE is solved numerically for a plane wave input with an imposed low level broadband noise background.We used a common optical noise model corresponding to a one photon per mode background [47], but in fact similar chaotic dynamics in MI can be seen with essentially any class of random amplitude and/or phase fluctuation on the input [25].The spatiotemporal intensity dynamics of |ψ(ξ, τ )| 2 for this case are shown in Fig. 4(a) and for completeness, we also show in Fig. 4(b) the associated spectral evolution [19].
We clearly see how the input plane wave evolves into a series of localized peaks, displaying both random temporal (transverse) and spatial (longitudinal) periodicity.Maximum gain for the spontaneous instability is at sideband frequency Ω = 1 and is associated with the initial emergence of Akhmediev breathers with temporal periodicity of ∆τ = 2π.After this initial stage, subsequent evolution is plotted up to ξ = 120.We also see how the incoherent temporal evolution is reflected in the frequency domain with chaotic spectral expansion and contraction, as the random emergence of particularly high intensity temporal peaks of ultrashort duration is associated with broader spectra.
Analyzing the evolution in terms of dominant balance yields the results shown in Fig. 4(c).
The color scale is the same as the previous figures.Comparing these results with the analytic SFB structures above allows us to distinguish the emerging localised structures.Indeed, in this case of highly random MI dynamics data-driven dominant balance successfully finds the Akhmediev breathers with period ∆τ ≈ 2π (for example, at ξ ≈ 11 and ξ ≈ 38), and we also see how propagation is associated with various ξ-periodic structures, breather collisions and Peregrine soliton-like rogue wave structures (e.g. the isolated feature in Fig. 4(c) at ξ ≈ 93).Being based on unsupervised clustering of contributing terms to the evolution equation rather then simple intensity thresholding, the technique successfully identifies developing localised structures even in low-intensity regions.This suggests the further application of the method in automated identification of emerging rogue wave structures [48].

DISCUSSION AND CONCLUSIONS
In conclusion, these results have clearly shown how the dominant balance approach provides a powerful tool for studying the interactions between dispersion and nonlinearity in the context of breather and modulation instability dynamics.In particular, even though these processes have been the subject of much previous study, visualising the dynamics with dominant balance segmentation clearly provides valuable insights into the relative contributions of different physical processes at different points in the evolution map.
We stress here, however, that data driven methods are not designed to replace existing techniques of analysing nonlinear dynamics, but should be seen as complementary tools to assist the use of physical considerations.For example, of particular interest is the way in which the dominant balance technique correctly associates the evolution of the plane wave background with a nonlinear stage of propagation.illustrates how simplistic interpretations such as associating nonlinear evolution with intensity thresholding could be misleading, and it is always necessary to consider the relative contributions of nonlinearity and dispersion in discussing the dynamics of the NLSE.
Moreover, whilst with experience, inspection of spectral and temporal evolution maps of the NLSE can allow some processes (such as collisions) to be readily associated with combined nonlinear and dispersive interactions, such interpretations can sometimes be misleading.This is particularly the case with generalized forms of the NLSE where multiple processes combine, as previous studied in Ref. [33] for the case of optical fibre supercontinuum generation.The strength of the dominant balance approach is that it provides additional information in an unsupervised manner (i.e.not based on intuition or experience).When applied in parallel with other analysis techniques, this provides important complementary information to yield the best possible physical interpretation of complex evolution.
Finally, we note that the NLSE describes propagation in many systems other than optical fiber, and there has been a strong recent focus on studying novel NLSE dynamics in deep-water hydrodynamics [22].In this context, we anticipate an important area of future application will be the case of MI induced by localized perturbations [49], and the associated emergence of rogue wave statistics [50,51].There is clearly much potential for data-driven discovery methods to be applied in NLSE-related systems.[25].

Equation space representation
The methodology of identifying a dominant balance model for a physical system at a particular stage of propagation aims to find a subset of terms of a more broadly applicable propagation model that locally dominates the dynamics.Following the approach and notation of Ref. [33], we consider a general evolution equation on a domain (ξ, τ ) written as follows: where K is the number of terms, and the terms f i can be constructed in various ways from the spatio-temporal field ψ(ξ, τ ).As discussed in Ref. [33] (and its accompanying Supplementary Information), the advantage of this implicit form of the propagation equation is that it stresses the balance that must be present to satisfy the equality: the sum of all the terms must be zero."Dominant balance" describes the situation when only a subset p of the K terms dominate the equality such that the contributions from the other K − p terms are small or negligible.Geometrically, the equation space is described by a vector:

Finding Dominant Balance models through clustering
The search for dominant combinations of terms within a higher-dimensional equation space is an ideal problem for unsupervised clustering algorithms [28,40].In particular, we use the algorithm and code package described in Ref. [33] and Ref. [41] respectively which are based on a probabilistic Gaussian Mixture Model (GMM) framework.GMM seeks to locate clustered subpopulations within an overall population of data, under the assumption that the data consists of a mixture of Gaussian distributions with specified weights, means and covariance matrices (usually denoted π k , µ k , Σ k respectively, where k is the cluster index).The covariance matrix here generalizes the usual variance of a one-dimensional Gaussian distribution to higher dimensions.In contrast to simpler techniques such as k-means associated with hard partitions between clusters, GMM describes membership of a clusters in a probabilistic sense, allowing the algorithm to fit and return clusters that overlap.The GMM algorithm is based on the expectation-maximisation technique, a standard approach that is fully described in e.g.Ref. [40].The particular GMM algorithm used here is GaussianMixture from the scikit-learn Python package [52], as implemented in Ref. [41].
A key motivation to use the GMM is that the covariance matrices can be interpreted physically to identify combinations of terms that dominate the dynamics.In particular, clusters associated with directions (dimensions) with significant variance correspond to physical terms that contribute actively to the dynamics (see the discussion of the results in Figs 1-3 above).However, there are some important additional factors that need to be considered to apply this approach successfully.
In particular, since the data points in the equation space may not actually approximate a mixture of Gaussian distributions, the algorithm will usually return a number of clusters greater that the number of physical balance regimes.As described in detail in Ref. [33], this problem can be overcome using Sparse Principal Component Analysis (Sparse PCA) which uses l 1 -regularisation to determine a sparse approximation to the leading principal component of each cluster [53,54].
In this case, when a particular cluster is associated with a dominant balance regime, it should be well described by the particular direction of its maximum variance.Note that l 1 -regularisation in this context is a standard approach in machine learning using the l 1 norm as the penalty in the PCA regression-optimization problem [53].
There are two key parameters that need to be selected to ensure that the returned models correspond as accurately as possible to physical regimes of dominant balance.The first is the particular number of clusters used in the Gaussian Mixture Model.Although in principle we can already anticipate the maximum number of potential clusters based on the number of terms in the propagation model, it is usually advantageous to initially choose a greater number, as the l 1 -regularisation step will later group together clusters found to possess the same sparcity patterns (i.e.reduced variances in the same directions of equation space) [33].The second parameter is associated with the sparse regularisation of the PCA that describes the tradeoff between accuracy and sparsity in the returned models.A procedure for this selection process is described in detail in the Supplementary information of Ref. [33], and is based on considering a returned Pareto-type curve that plots the residual error of the inactive terms (accuracy) against the regularization parameter (sparsity).It is generally straightforward to see from this plot the most suitable parameter to generate the returned balance model.The very last step of the algorithm involves re-mapping the sparse clusters back onto the original spatio-temporal domain, and it is at this point we can directly compare the initial field distribution with the identified cluster map (as in Figs 1-4).
It is useful to give further numerical details for our results.For the three classes of soliton on finite background considered in Figs 1-3, the evolution maps ψ(ξ, τ ) were computed over (N ×M ) = (501 × 1024) in ξ and τ respectively.For the noise-driven map considered in Fig. 4, evolution was computed over (N ×M ) = (5001×1024) in ξ and τ respectively.The GMM search was based on an initial selection of up to 5 clusters and the sparse regularisation parameter α (used in the Python function SparsePCA [41]) was in the range 50-100.We also note the computation time associated with the GMM clustering and SPCA analysis, which was typically 6 and 21 minutes respectively for solitons on finite background and noise-driven MI, running on a standard Windows PC with 3.00 GHz 6 MB cache double-core CPU.
Fig.1(b-i) plots the identified clusters in the three-dimensional space of the real parts of coordinates {iψ ξ , ψ τ τ , ψ|ψ| 2 }, whereas Fig.1(b-ii) and Fig.1(b-iii) show two projections as indicated.The color key corresponds to two different dominant balance models that are found: one where only the nonlinear and propagation terms contribute (blue) and another where all NLSE terms contribute (orange).No cluster is found that involves only the dispersive and propagation terms.Note that for convenience we plot dependencies only for the real field components, but similar results are found for the imaginary components.The results in Fig.1(b)show that all the points assigned to the blue cluster (nonlinear and propagation terms) are strongly localised in the equation space forming a dense distribution that manifests nearly zero variation with respect to the ψ τ τ axis (see particularly Fig.1(b-iii)).In contrast, the orange cluster (all terms) is distributed throughout the equation space with no reduced variance with respect to any of three axes.This illustrates the geometrical interpretation lying behind the dominant balance approach.

FIG. 4 .
FIG. 4. (a) Spatio-temporal evolution of the normalized intensity for spontaneous modulation instability.(b) The associated spectral evolution.(c) The results of dominant balance revealing the different interaction regions according to the colormap shown (same as in previous figures).
τ m ), ...]] T where each of the dimensions (directions) corresponds to a specific term in the evolution equation (here indices n ∈ [1, N ] and m ∈ [1, M ] represent the discretization of ψ(ξ, τ ), where N and M are the number of points in the ξ and τ directions respectively).A dominant balance regime then has a direct geometrical interpretation -dynamical points attributed to a certain dominant balance regime will be restricted to p directions of the full K-dimensional space.In other words, when plotting the different terms in the equation space, the points associated with the dominant p terms will have significantly reduced variance with respect to other K − p directions.In the case of the NLSE, the dimensionality K = 3 and each dynamical pointψ(ξ n , τ m ) is associated with a vector [iψ ξ (ξ n , τ m ), ψ τ τ (ξ n , τ m ), |ψ(ξ n , τ m )| 2 ψ(ξ n , τ m )] T .In geometrical terms, dominant balance between the propagation and nonlinear Kerr terms (iψ ξ , |ψ| 2 ψ) will be represented by an ensemble of points restricted on a line with near-zero variance with respect to the dispersion term ψ τ τ (e.g. the blue clusters in Figs.1-3).In contrast, the ensemble of points distributed throughout the iψ ξ + ψ τ τ + |ψ| 2 ψ = 0 plane will represent the full dynamics that involves the interplay of all three dynamical terms (e.g. the orange clusters in Figs1-3).