Thresholds between modulational stability, rogue waves and soliton regimes in saturable nonlinear media

We study a third-order nonlinear Schrödinger equation for saturable nonlinear media, which can represent either optical pulse propagation in 1D nonlinear waveguide arrays or electron dynamics in one-dimensional lattices. Here, we describe how stable uniform solutions can evolve into regimes in which they are localized and the influence of saturation on this transition. We observed regular and chaotic-like breathing dynamics as intermediate regimes, which have their domain of existence significantly altered by the saturation parameter. Critical nonlinear strengths separating the existing regimes are shown in phase diagrams, evidencing the role played by saturation. Numerical data and analytical approach show the nonlinear strength above which uniform solutions become breather solutions increasing with the saturation parameter. In the regular breathing regime, we reveal the breathing frequency for media with saturable nonlinearity displaying faster growth as it moves away from the critical point of uniform solutions. Furthermore, critical nonlinear strength separating the regimes of regular and chaotic-like breather solutions exhibits a decreasing behavior as the saturation parameter increases. The latter ones present clear signatures of emerging rogue waves, such as peaks showing long-tailed statistics. Thresholds of this regime are increased by the saturable nonlinearity. Thus, the regime of localized solutions, which we have shown to be well-described by bright soliton-like structures, becomes less accessible with increasing saturation.

tions increasing with the saturation parameter. In the regular breathing regime, we reveal the breathing frequency for media with saturable nonlinearity displaying faster growth as it moves away from the critical point of uniform solutions. Furthermore, critical nonlinear strength separating the regimes of regular and chaotic-like breather solutions exhibits a decreasing behavior as the saturation parameter increases. The latter ones present clear signatures of emerging rogue waves, such as peaks showing long-tailed statistics. Thresholds of this regime are increased by the saturable nonlinearity. Thus, the regime of localized solutions, which we have shown to be well-described by
Modulational instability refers to the stability, with respect to any perturbation, of a wave with constant amplitude that propagating through a nonlinear dispersive medium [27,28]. Such phenomenon has been widely studied in different frameworks and is usually seen as the starting point for regimes of breathing modes and localized solutions [29][30][31][32][33][34][35][36]. Significant changes in the crossover between modulational instability and breather solutions have been reported by exploring the influence of next-nearestneighbor coupling in nonlinear discrete lattices [29]. Breather solutions with dark-and bright-type profiles are observed depending on the strength of such next-nearest-neighbor interaction. From the perspective of Kerr-like media, modulational instability has been described as the mechanism by which bright and dark soliton-like wave functions can be generated in Bose-Einstein condensates in one-dimensional optical lattices [30]. Soliton solutions coming from modulational instability have experimentally been reported in real nonlinear electrical lattices, with results exhibiting full agreement with theoretical analysis based on the nonlinear Schrödinger equation [31]. Bright and dark soliton solutions have been reported in a medium with cubic-quintic nonlinearity, in which an increase in modulational instability gain and a decrease in the amplitude of solitons are shown as a result of increasing nonlinearity [32]. Studies on 1D and 2D nonlinear lattices show the critical nonlinear coupling, above which the modulational instability occurs, remaining finite for square lattices and inversely proportional to the length for 1D lattices. A direct transition between stable uniform wave functions and asymptotically localized solutions is reported for square lattices, while linear chains exhibit an intermediate regime in which the wave function dynamics develop complex breathing patterns [35]. The interplay between modulational instability and localized solutions has been evaluated by considering the presence of long-ranged hopping amplitudes in one-dimensional lattices, a prospect that influences the effective dimensionality of the system [36]. Results show breathing and chaotic-like solutions, emerging as intermediary patterns between stable uniform and localized solutions, are fully suppressed for hopping amplitudes decaying slower than 1/r 2 . Persisting localized solutions and universal features of oscillatory structures have been experimentally observed in nonlinear Kerr-like media undergoing modulational instability, demonstrating the potential of optical fiber experiments in proving such phenomena [37,38].
Previous results show that stable uniform solutions evolve for localized solutions by undergoing different intermediate dynamic behaviors. However, such observations were restricted to systems with cubic (Kerrlike) nonlinearity, giving rise to the following question: How do systems in which the nonlinearity cannot be described as Kerr-like behave? Given the range of systems in which the effective evolution is governed by nonlinear equations, we observe other possible scenarios of nonlinearity, such as cubic-quintic law [39,40], dual-power law [41,42], log [43,44] and saturating laws [45][46][47]. This latter, which is well-known in the modeling of semiconductor-doped glass [45][46][47], has been shown to have nontrivial impacts on wave dynamics. For example, it may lead to the formation of two stable soliton pulses with the same time duration [47] or promote opposite impacts on the rectifying action for nonreciprocal transmission through an asymmetric dimer [48]. Thus, we study how saturable nonlinearity influences the interplay between modulational stability and localized regimes. The set of nonlinear equations governing the system is described in Sect. 2, as well as the initial conditions and solving process. In Sect. 3, we show regular and chaotic-like breathing dynamics emerging as intermediate regimes between uniform and localized solutions. In addition to characterizing the existing dynamical regimes, in which emerging rogue waves and bright soliton-like structures are unveiled, we offer phase diagrams that reveal critical nonlinear strengths separating them, evidencing the role played by saturation. Summary and final considerations are shown in Sect. 4.

Model
The problem consists of analyzing a set of onedimensional discrete nonlinear Schrödinger equations, written as (1) Such equations have been applied for optical pulse propagation in 1D equidistant nonlinear waveguide arrays [25] and electron dynamics in one-dimensional lattices [49]. Hereafter, the propagation coordinate t is time, the parameter α denotes the on-site energy at site n and β stands for the coupling between adjacent sites. χ represents the strength of a third-order nonlinear response, while ζ is the parameter that controls the saturation of this nonlinear response. We are interested in low saturable media ζ 1.5. We also consider the absence of disorder, avoiding related features to Anderson localization [50]. Thus, α and β are constant, with α taken as the reference energy (α = 0) without any loss of generality. β = 1, which implies that χ is in units of β.
We concentrate our study by considering the initial wave packet composed of a uniform wave of amplitude ψ 0 = 1/ √ N superposed to a very weak noise = 10 −3 ψ 0 , such that the initial amplitudes at each site n of the lattice with N sites are randomly distributed in the interval [ψ 0 − , ψ 0 + ]. After generating all site amplitudes, a proper normalization is further employed to warrant the unitary norm of the resulting wave packet (total power conservation). The above set of differential equations was solved by using a standard eighth-order Runge-Kutta method, with a time step that is small enough to maintain the wave function norm conservation (|1 − n |ψ n | 2 | ≤ 10 −14 ) along the entire time interval considered. Just as in Ref. [35,36], periodic boundary conditions are assumed.

Results and Discussions
We start by following the time-evolution of the wave packet probability distribution |ψ n (t)| 2 on a chain with N = 100 sites. In Fig. 1, we show the dynamics in systems ruled by representative values of the nonlinear parameter χ and saturation ζ . In Fig. 1a-c, we illustrate different regimes in ascending order of nonlinearity in absence of saturation coefficient (ζ = 0). For very small nonlinearities (Fig. 1a), the wave packet remains uniformly extended over the entire lattice, signaling the stability of the uniform solution. A wave packet dynamics with regular breathings are obtained by increasing the nonlinear parameter (see Fig. 1b). This result shows the critical nonlinearity, above which the uniformly stable solution does not survive, as being 0.195 < χ BS < 0.200, in full agreement with previous reports [35,36]. A further increase of χ can drive the wave packet into a localized regime, as shown in Fig. 1c. Such behavior is also in full accordance with previous studies [35,36], which report this interplay. Here, by considering saturable nonlinear media, we explore the influence of saturation (ζ ) and observe significant changes in the dynamic behavior. For a nonlinear response χ just below the critical point of modulational instability of non-saturable media (Fig. 1a), we observe the survival of a uniform regime (see Fig. 1d, g). However, for a χ just above the critical point of modulational instability of non-saturable media (χ = 0.200), the wave packet can exhibit a weakening or even a suppression of regular breathing dynamics as a function of ζ (see Fig. 1e, h). The saturation effect can lead to a change in dynamical regimes even for strong nonlinearities, exemplified here by the localized regime ( Fig. 1c) that gave way to the regime of irregular (chaotic-like) breathings (Fig. 1f, i).
Thus, we observed the saturation of third-order nonlinearity exhibiting a significant contribution, with clear changes in the thresholds of the existing dynamical regimes: uniform, breathing, chaotic-like, and localized. To better understand and characterize such changes, we compute the time-evolution of the participation function, defined as ( This measure gives an estimate of the number of lattice sites over which the wave packet is spread at time t, with P(t) ∝ N 0 describing a localized wave packet and P(t) ∝ N a wave packet extended over the lattice. In Fig. 2a, we show the time-dependent normalized participation number [P(t) = P(t)/N ] for the same configurations used in Fig. 1. We observe the participation function recovering relevant aspects reported before. The participation function is equal to the lattice size along the entire time-evolution for the regime of stable uniform solution, while the breathing regime is signaled by regular breaths of the participation function which always visits the fully extended wave packet condition. In the chaotic-like regime, the participation function exhibits a series of irregular breathings, with different minimum and indefinite frequency. The reference for the breathing dynamics is no longer the fully  extended wave packet. The localized solution is characterized by P(t) close to zero emerging after an initial transient.
At first, we focused on determining the crossover from a stable uniform regime to a regular breathing regime. By evaluating the participation function for systems with different χ and ζ , we report in Fig. 2b a phase diagram with the critical nonlinear strength (χ BS ) separating such regimes. We observe χ BS scaling linearly with the saturable parameter. Such behavior is in agreement with the analytical approach, in which we follow the standard procedure for investigating the stability of a continuous wave solution [27,34]. We start by observing that the continuous version of Eq. 1 supports a continuous-wave solution Just like in the numerical analysis, we impose a slight perturbation in its amplitude: Then, the perturbation (x, t) is examined by utilizing a standard linear stability analysis. By considering the linearization as a function of the initial state [ (x, t) ψ 0 ], we obtain Here, * denotes the complex conjugate of , and the abbreviation = (x, t) has been used. Since a white noise random perturbation exhibits a wide spectral range, which includes all harmonic contributions, we study solutions (x, t) = Ae i(kx− t) + Be −i(kx− t) , where A and B are the amplitudes of the weak modulation, k is the modulation wave number, and is the frequency. Thus, substituting the last one into Eq. 4 and considering its nontrivial solution, we obtain the dispersion relation Such expression determines how time oscillations are associated with spatial oscillations of wave number k.
We observe the stability of the stationary regime for disturbances with large wave vectors since remains real. In contrast, ∈ I m reveals an exponential amplification of all wave vectors with which results in the breakup of continuous waves and therefore describes the uniform solution becoming unstable. Given the boundary conditions, the minimum wave number for a harmonic perturbation is of the order of 2π/N . Since the wave-function normalization for a uniform solution provides ψ 0 = 1/ √ N , we obtain the characteristic nonlinear strength (χ BS ) above which the solution is unstable as with the term proportional to 1/N 3 being neglected.
Although we use a continuous approach to derive an analytical evaluation, such expression corroborates numerical results shown in Fig. 2b, which displays the critical nonlinearity of uniform solutions regime increasing linearly as a function of the saturation parameter. The Eq. 7 also informs a scaling behavior with the lattice size, which is in full agreement with numerical data (see Fig. 2c). This critical nonlinearity decreases for larger system sizes, providing modulational instability of the uniform solution in the limit N → ∞ for any finite nonlinear strength. Given the dependence of the lattice size, the influence of saturation ζ on the modulation instability threshold becomes more pronounced for smaller N . Just above the modulational instability threshold, the wavefunction develops regular breathing oscillations between the uniform and the localized state in a finite segment. Figure 3a shows such behavior by using the normalized participation function. Breathing patterns exhibit larger amplitudes and shorter periods (T) as χ moves away from critical nonlinearity χ BS . However, this increase in frequency (ω = 1/T) exhibits distinct behaviors when comparing systems with and without saturable nonlinearity. Figure 3b shows saturable nonlinear media exhibiting a faster frequency increase, which suggests the saturation also influences the intermediate regimes between the uniformly stable and the localized solutions.
Based on previous results, which show the participation function characterizing the existing dynamical regimes, we computed the local minima of the normalized participation number for a large set of nonlinear strengths. By considering possible initial transient time, only minima occurring in the interval 80, 000 < t < 100, 000 are written down. Thus, we show in Fig. 4a plots of such minima as a function of χ for distinct saturating parameters (ζ = 0.0, 0.25, 0.50, 0.75). Note the Fig. 3 a Time evolution of the normalized participation function [P(t) = P(t)/N] for representative settings within the breathing regime. The period of the breaths decreases as we move away from the modulational instability threshold. b Breathing frequency for media with saturable nonlinearity exhibits a faster growth as it moves away from the critical point minima appearing just above the modulational instability threshold, displaying smaller and smaller P min as we increase the nonlinear response χ . Consistent with the development of regular breathings described before, each χ exhibits only one P min . However, a further increase in nonlinearity promotes the emergence of distinct breathing modes with different periods. The periodic pattern is destroyed by interference and differentorder breathing modes appear. Now, the fully extended solution is not always visited after a breath and the wave packet dynamics becomes very sensitive to the initial conditions. Such regime, so-called chaotic-like [35], becomes evident when we observe participation minima spreading over a wide range of values, signaling incommensurate breathing modes along dynamics (see Fig. 4a). In Fig. 4b we show the threshold between regimes of regular and chaotic-like breathings as a function of saturable parameter ζ . Now, we observe the critical nonlinearity (χ C L ) decreasing as the saturation increases. Since the frequency of breathing modes increases faster in saturable media, the interference between different breathers is favored, corroborating the narrowing of the regular breathing regime.
As seen in previous results, the chaotic-like regime is characterized by different-order breathing modes with an undefined pattern. Solitary peaks with high amplitude can therefore emerge in such dynamical regime. This framework, combined with the current hot topic of rogue waves and the long-standing discussion of them being fundamentally driven by nonlinear or linear processes [51,52], led us to examine the emer-  Fig. 5a we show an example: a snapshot of spacetime evolution of the wave packet probability |ψ n (t)| 2 for a system with N = 100 sites, χ = 3.50 and ζ = 0.50. We observe a typical rogue wave event, with a large |ψ n (t)| 2 in comparison to the surrounding profiles. Spatial and temporal cross-sections give a more detailed analysis. Unpredictability and short-lived are observed features that corroborate the rogue wave scenario. Its amplitude is approximately 3 times bigger than the average height among one-third of the highest breaths recorded in a complete time-evolution, overcoming the usual requirement (2 times bigger -dashed line) defining rogue waves [53,54]. Here, such rogue wave threshold was computed by recording amplitudes of breaths in space and time for an ensemble of 40 independent numerical experiments. A better characterization is shown in Fig. 5b, in which we display normalized probability density functions (PDFs) of |ψ n (t)| 2 derived from the ensemble. PDFs exhibit a L-shaped long-tail distribution, which is typical of a rogue wave regime [51,53,54]. Despite the representative values of χ and ζ used in Fig. 5, this behavior is characteristic of the chaotic-like regime, both in media without and with saturable nonlinearity. Here, rogue waves have the initial random noise and the nonlinear influence as underlying mechanisms, not a disordered media as explored in recent studies [51,55].
Finally, we investigate the regime of localized solutions, achieved by self-trapping phenomena, whose underlying mechanism is the non-dispersive character of nonlinearity and the emergence of sufficiently high breaths. Such dynamical regime is reported for strong enough nonlinearities and previous results clearly suggest that the threshold between chaotic-like and localized solutions is changed by saturable parameter. Looking at Fig. 4a, the minima of normalized participation function are on the order of a few sites (P min ≈ 0) for strong enough nonlinearities. By evaluating the average of asymptotic participation function P(t → ∞), we identified the threshold to self-trapping regime increasing with increasing ζ (see Fig. 6a). This result indicates saturation as a relevant parameter, which postpones the emergence of spatially localized solutions as we increase third-order nonlinearity. In Fig. 6b we show spatial profiles for two representative settings of χ and ζ that meet the localized regime. Dashed lines represent the cross-section. These localized profiles exhibit the well-known soliton shape [1,56]. Dashed lines represent fittings using |ψ n | 2 ∝ sech 2 [γ (n − n max )], characterizing both profiles as bright soliton-like solutions. Spatial profiles of such localized structures in diverse other configurations of ζ and χ exhibit the same appearance.

Final considerations
In summary, we bring results that contribute to a better understanding of the crossover between modulational stability and soliton formation in discrete systems. We studied stable uniform solutions evolving to different dynamical regimes as we increase nonlinearity. The system is described by a set of thirdorder saturable nonlinear Schrödinger equations, which can represent either optical pulse propagation in 1D nonlinear waveguide arrays or electron dynamics in one-dimensional lattices. By considering a small noise on the initial uniform wave packet, we describe how saturation changes the thresholds of existing dynamical regimes: uniform, regular breather, chaotic-like breather, and localized solutions. Numerical and analytical studies show the nonlinear strength above which uniform solutions become breather solutions increasing with the saturation parameter. Thus, we observe the saturable parameter contributing to the stability of the uniform solution. By exploring the crossover between regular and chaotic-like breathing regimes, we report the inhibition of regular breathing regime and the widening of the chaotic-like regime. The wave packet dynamics in the latter one exhibit clear signatures of rogue waves, with a greater number of high-amplitude modes than would be obtained from Gaussian statistics.
The regime of localized solutions becomes less accessible with the saturation of nonlinear response since the critical nonlinearity increases linearly with saturable parameter. Such localized solutions are well described by bright soliton shape. Finally, it would be interesting to have future contributions exploring other scenarios, such as higher-order nonlinearities or the presence of time-and space-dependent dissipation. Such results within an analytical framework would also be interesting and would certainly bring valuable insights into the overall wave dynamics in discrete nonlinear lattices. The dynamic regimes and thresholds revealed here can still be combined with machine learning and artificial neural networks, contributing to the development of heterogeneous deep computational systems [57].
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.