The local wave phenomenon in the quintic nonlinear Schrödinger equation by numerical methods

The nonlinear Schrödinger hierarchy has a wide range of applications in modeling the propagation of light pulses in optical fibers. In this paper, we focus on the integrable nonlinear Schrödinger (NLS) equation with quintic terms, which play a prominent role when the pulse duration is very short. First, we investigate the spectral signatures of the spatial Lax pair with distinct analytical solutions and their periodized wavetrains by Fourier oscillatory method. Then, we numerically simulate the wave evolution of the quintic NLS equation from different initial conditions through the symmetrical split-step Fourier method. We find many localized high-peak structures whose profiles are very similar to the analytical solutions, and we analyze the formation of rouge waves (RWs) in different cases. These results can be helpful to understand the excitation of nonlinear waves in some nonlinear fields, such as the Heisenberg ferromagnetic spin system in condensed matter physics, ultrashort pulses in nonlinear optical fibers, and so on.

Recently, the specific question of the formation of RWs in the framework of integrable turbulence has been widely studied. In 2009, Zakharov [17] proposed the concept of integrable turbulence to describe a kind of chaotic wave field evolved from nonlinear integrable system when the initial conditions are disturbed by random factors. In fact, the integrability of the governing equations makes it possible to obtain some exact solutions as the basic nonlinear modes of these systems. For example, solitons and breathers are the main modes with significant amplitudes, while radiation waves can be ignored because of their low amplitudes, although they are also the modes that make up the chaotic background field. However, despite knowing the explicit forms of these modes, the complicated chaotic motion that involves the diversity of their interaction in a generally irregular way is still an unsolved problem. Therefore, integrable turbulence is mainly analyzed by numerical simulation [18,19]. The simplest chaotic mode will be formed when multiple solitons or breathers with different amplitudes move in all possible directions. The collision between them will produce a high amplitude local structure, that is, the socalled rouge wave. This kind of chaotic wave field will evolve into soliton turbulence, breather turbulence or mixed turbulence with different values of parameters in the initial conditions [11,20]. Some studies have investigated the physical mechanisms of RWs in the focusing NLS equation from the different viewpoints [20][21][22][23][24][25]. Particularly, the generation mechanism of RWs in the focusing NLS equation and Hirota equation has been studied by inverse scattering analysis method from the viewpoint of their spectral problems [20,26].
In this paper, we focus on the integrable quintic NLS equation [27], which can be used to describe the Heisenberg ferromagnetic spin system in condensed matter physics and ultrashort pulses in nonlinear optical fibers [28,29].
Here t is the propagation variable, x is the transverse variable and the modulus of complex valued function u(x, t) represents the envelope of the nonlinear wave. The real constant δ stands for the fifth-order dispersion. Some scholars have investigated this equation, but they only consider all kinds of exact solutions and conservation laws [28,[30][31][32][33][34]. In this paper, we will study the numerical evolutions under random initial conditions and analyze the dynamical behaviors in the framework of integrable turbulence. In addition, we verify this phenomenon qualitatively by algebraic method.
Eq. (1) is a compatibility condition, xt = t x , of the system of two linear partial differential equations with variable coefficients (also known as the Lax pair) [35], where λ is a complex spectral parameter, is a column vector given by with 1 (x, t) and 2 (x, t) being complex valued functions, and matrices σ, Q, T j ( j = 0, · · · , 5) are given by where C 0 and C 1 are arbitrary real constants, and the superscript "*" represents the complex conjugate. The rest of this paper is arranged as follows. In Sect. 2, we numerically study the spectral signatures of the spatial Lax pair (2a) with distinct potentials of the quintic NLS equation through the Fourier collocation method. In Sect. 3, we discuss the formation of rogue waves in the quintic NLS equation. Distinct nonzero initial conditions are used to analyze the generation mechanism of rogue waves. In Sect. 4, we summarize the paper.

Fundamental analytical solutions
According to the previous researches, Eq. (1) has been proved to possess the following analytical solutions [35,36]: (i) The bright soliton [see Fig. 1a]: for purely imaginary eigenvalues, λ = bi, where b determines the amplitude of the soliton. In this case, the soliton has nonzero velocity, v = −16δb 4 . Additionally, the phase term e 2ib 2 t evolves along the t axis and does not depend on x. When δ = 0, it degenerates to the soliton solution of classical NLS equation.
(ii) The Akhmediev breather [see Fig. 1b]: where is the growth rate of periodic modulation, and (iv) The Rouge wave (RW) [Localized in space and time, see Fig. 1d]: which can be obtained from the limit a → 1 in the AB (6) or c → 1 in the KM breather (7), and it is a regular rational solution with one maximal point and two minimal points.

Numerical computation of spectral portraits of the spatial Lax pair with distinct potentials
As far as we know, there is no research on the spectral signatures of spatial Lax pair (2a) with distinct potentials of Eq. (1). Usually, it is difficult to analytically find the spectral portrait for a general potential u(x, t). Therefore, some numerical methods have been proposed to study the IST spectra [37,38]. Eq. (2b) can be rewritten as a standard linear eigenvalue problem [38], which can be subsequently solved by using the Fourier collocation method [39]. The x axis is truncated into a finite interval [−L/2, L/2], where L is the length of the interval. Then, for a fixed t = t 0 , the eigenvector = ( 1 (x), 2 (x)) T as well as the potential u(x, t) is expanded into Fourier series with 2N + 1 modes, where k 0 = 2π L and t = t 0 are given. Substituting these expansions into Eq. (9) and equating the coefficients of the same Fourier mode, one can obtain the eigenvalue spectrum. Figure 2 shows the eigenvalue spectra of Eq. (9) with distinct potentials u(x, t) (BS, AB, KM, RW). It can be seen that the BS (5) with b = 0.5 has discrete eigenvalues of ±0.5i, while the eigenvalue spectra of the other three potentials are all over the (−1, 1) interval of the imaginary axis. In addition, the AB (6) with a = 0.75 has discrete eigenvalues around ±0.75i, the KM (7) with c = 1.2 has discrete eigenvalues of ±1.2i, and the RW (8) has discrete eigenvalues of ±i, respectively. Furthermore, it is worth noting that in our numerical results, although the potential u(x, t) varies with the change of time t, the corresponding eigenvalue spectrum is always invariant.   Fig. 1 at t = 0, 0.5, 1; a2-d2 the corresponding spectral portraits. A numerical interval of length L = 100 discretized by 1000 points has been used to compute the eigenvalue spectra

Spectral portraits of periodized structures
In order to overcome the fact that a satisfactory numerical IST analysis of RWs cannot be generally achieved in a local way from a single isolated object, we usually need to use a local periodization procedure [22]. In this way, the spectral portraits can still be accurately determined from the numerical IST analysis of the periodized structures obtained from proper periodization of the isolated object.
Similar to [20,26], we consider spectral portraits of periodized wavetrains in the Eq. (1). We extend the central part of the analytical solutions periodically at t = 0, including the local peaks and the surroundings. In this way, the nonlinear interaction between the central peak and the surrounding structure can be effectively guaranteed, and the effective intensity of the interaction is determined by the spatial period T that is used to produce the periodized wavetrains.
From Fig. 3, we can find that after the periodized extension of the central parts of the BS and KM breather solutions, although their potentials become more complex, the corresponding eigenvalue spectra do not change significantly, except that small bands are now found in Fig. 3a3 and b3 instead of single points. The discrete eigenvalues with a maximal Fig. 3 The numerical results of periodized structures. a1, b1 Profiles of |u(x, t)| 2 in Fig. 1a and c at t = 0; a2, b2 the corresponding periodized wavetrains (T = 6); a3, b3 the spectral portraits of the periodized wavetrains showing that the periodization procedures produce same eigenvalue bands norm in Fig. 2a2 and c2 are ±0.5i and ±1.2i, respectively, while the discrete eigenvalues are ±0.4871i, ±0.4916i, ±0.4974i, ±0.5034i, ±0.5083i, ±0.5116i in Fig. 3a3, and ±1.197i, ±1.198i, ±1.199i, ±1.2i, ±1.2i, ±1.201i in Fig. 3b3, respectively.

Numerical simulations via symmetrical split-step Fourier method: the formation of rouge waves
In this section, we simulate the evolution of Eq. (1) from distinct initial conditions by using the symmetrical split-step Fourier method [38]. We take δ = 10 −6 in the following numerical simulations, because Eq. (1) is very sensitive to the coefficient of the quintic terms.

The symmetrical split-step Fourier method
The symmetrical split-step Fourier method simulates the interaction between the dispersion term and the nonlinear term by assuming that they work independently. Therefore, Eq. (1) can be rewritten as follows: where M and N represent dispersion and nonlinear operators, respectively. The idea of symmetrical split-step Fourier method is shown in Fig. 4. The propagation process from t to t +h is divided into three parts, in which the dispersion and nonlinear effects are assumed to work independently. First, the envelope u(x, t) propagates half a step h 2 , in which only the dispersion effect is considered, i.e., N (u) = 0, and it can be transformed into ordinary differential equation to solve by Fourier transform; then, in the plane of t + h 2 , only the nonlinear effect works, i.e., M(u) = 0, and it can be solved by combining Fourier transform and four-order Runge-Kutta method; next, envelope propagates the remaining half step like the first part. In this way, we can obtain u(x, t + h) representing the envelope after the evolution of one step.

Dam-break problem
In this section, we first consider the dam-break problem [26,40] as the initial condition, to numerically simulate the wave propagation in Eq.
(1). The dam-break problem describes the evolution and interaction of two counter-propagating nonlinear wave trains generated in the quintic NLS box problem. We show that the interaction dynamics results in the emergence of modulated large-amplitude quasiperiodic breather lattices whose amplitude profiles are closely approximated by the Akhmediev breathers (AB) within certain space-time domain. [see Fig. 5]. Then, we take the localized structures observed at t = 7.51, 9.04, 10.70 and 12.74 as examples, and extend the central part of them periodically to study their eigenvalue spectra by the Fourier oscillatory method and numerical IST analysis. Note that the criterion for choosing the period T is that any sufficiently small change in T can only results in enough small quantitative change in the spectrum. From Fig. 6a3-d3, it can be found that these spectra are obviously different from Fig. 2(b2), which suggests that they are not identical to the analytical solutions. In fact, from the viewpoint of finite-gap theory, all the IST spectra possess 3 main spectral bands, which indicates that they are non-degenerate genus 2 solutions, while the AB, KM and RW are all degenerate genus 2 solutions [22].

Noise-driven modulation instability and the formation of rogue waves
For the nonlinear wave equations, the modulation instability (MI) may be regarded as a necessary condition for generating the extreme waves. In this section, we consider the plane wave plus a random noise as the initial condition, where η is a random noise uniformly distributed in the (−1, 1) interval and is the intensity coefficient of the noise, to numerically simulate the wave evolution in Eq. (1). We observe that compared with Fig. 5, the noise-driven MI results in a series of richer high-peak structures with random intensity in chaotic field, whose amplitude profiles are closely approximated by the AB, KM and RWs [see Fig. 7]. In the evolution process, a series of irregular localized peak structures appear around t = 3. Then, more complex periodic growth and decay behaviors emerge along the t axis. Notice that the average spatial period here is around x = √ 2π , which corresponds to the reciprocal of the maximum MI gain frequency. Moreover, we compare the profiles of some high intensity peak structures with the corresponding analytical solutions in Eqs. (6)- (8) and find that they agree very well.  Although it is impossible to observe the ideal analytic structures under the given random initial conditions, it is still worth noting how the analytic solutions can be closely mapped to the structures generated by the noise.

Complex random initial condition and the formation of integrable turbulence fields
However, the random noise η in Eq. (13) is not related to the x. In order to make the curve smoother, we take the one-dimensional random rough surface as the initial condition, where x ∈ [−L/2, L/2], L = 100, with a uniform grid of N = 1000 nodes, k j = 2π j/L( j = −N /2 + 1, · · · , N /2), k = 2π/L, F(k j ) is defined as the Fourier transform of f (x n ) as follows: N (0, 1), where S(k j ) is the power spectrum function of random rough surface obeying Gaussian distribution, and its expression is where μ and l represent standard deviation and correlation length, respectively. In this way, we construct a complex random function f (x), whose real and imaginary parts are independent random variables, each obeying the Gaussian distribution with standard deviation μ and correlation length l. The use of two variable parameters enables us to reveal more complex integrable chaotic dynamics than ever before. Figure 8 shows the effect of these two variable parameters on the initial conditions. The random deviations from the cw with amplitude 1 increase with the increase of standard deviation μ [see Fig. 8a]. However, when the correlation length l increases, the deviations don't change significantly, but the period of the random waves becomes larger [see Fig. 8b]. Figure 9 exhibits the numerical results for the initial conditions in Fig. 8a. We numerically solve the eigenvalue problem Eq. (9) using these initial conditions and the sets of complex eigenvalues λ are shown in Fig. 9a1-c1. We can see that when μ = 0.1, eigenvalues are almost located on the imaginary axis, which correspond to the excitation of ABs. With the increase of the value μ, the eigenvalues gradually deviate from the coordinate axis and turn to solitons, meaning that the chaotic pattern gradually changes from the breather turbulence to the soliton turbulence, which is consistent with the actual simulation results in Fig. 9a2-c2. Figure 9a3-c3 displays part of the three-dimensional diagrams in the evolution process, which clearly illustrate the collision process of breathers and solitons. It should be emphasized that the superposition of the surrounding small radiation waves leads to the random irregular variation of the amplitudes in the collision process.
Another factor that plays a significant role in this chaotic evolution is the correlation length l of the initial random field. Similarly, we exhibit the numerical results in Fig. 10 for the initial conditions in Fig. 8b. Compared with Figs. 9a1-c1, the eigenvalues in Figs. 10a1-c1 still remain near the coordinate axis with the increase of l, but the range on the imaginary axis gradually expands, which means that the amplitudes of solitons in the initial field gradually increases. Figure 11 illustrates the results of the simulations with the above initial conditions. It shows the peak amplitude (absolute maximum of |u| over the whole x interval of simulation) as a function of t. As expected from the theory of MI, the maximum increases exponentially at the initial stages of evolution, and then gradually enters into chaotic pattern. Moreover, with the increase of value μ or l, the convergence to the chaotic state is quicker and the amplitudes of RWs in the chaotic field are higher.

Conclusions
In this paper, we focus on the quintic nonlinear Schrödinger (NLS) equation. Firstly, we numerically study the spectral signatures of the spatial Lax pair   Distinct nonzero initial conditions are used to analyze the generation mechanism of rogue waves: (i) the one case is based on a decaying initial condition (dam-break problem). We observe the modulated breather lattices approximated by the AB, and further compare their spectral signatures using periodization procedure. (ii) in another case, we use plane wave plus a random noise as the initial condition such that the modulation instability (MI) seeded from the noise results in a series of richer high-peak structures of random intensity, and they agree well with the analytical solutions. (iii) in the last case, we take the one-dimensional random rough surface as the initial condition and analyze the formation of rogue waves. We find that two parameters of the initial random functions are essential to fix the number of rogue waves in the resulting chaotic wave field: μ influences the relative number of breathers and solitons in the resulting chaotic wave field while l mainly influences the amplitudes of the solitons. They all lead to an increase in the amplitude of rogue waves. The results obtained in this paper will be useful to understand the formation of rogue waves in the quintic NLS equation and other related models, which can provide ideas for the controllable excitation of rogue waves in nonlinear fields such as the Heisenberg ferromagnetic spin sys-