Numerical investigation of stability of breather-type solutions of the nonlinear Schrodinger equation

. In this article we conduct a broad numerical investigation of stability of breather-type solutions of the nonlinear Schrödinger (NLS) equation, a widely used model of rogue wave generation and dynamics in deep water. NLS breathers rising over an unstable background state are fre-quently used to model rogue waves. However, the issue of whether these solutions are robust with respect to the kind of random perturbations occurring in physical settings and laboratory experiments has just recently begun to be addressed. Numerical experiments for spatially periodic breathers with one or two modes involving large ensembles of perturbed initial data for six typical random perturbations suggest interesting conclusions. Breathers over an unstable background with N unstable modes are generally unstable to small perturbations in the initial data unless they are “maximal breathers” (i.e., they have N spatial modes). Additionally, among the maximal breathers with two spatial modes, the one of highest amplitude due to coalescence of the modes appears to be the most robust. The numerical observations support and extend to more realistic settings the results of our previous stability analysis, which we hope will provide a useful tool for identifying physically realizable wave forms in experimental and observational studies of rogue waves.


Introduction
Interest in understanding rogue wave phenomena has been steadily growing for the past decade, especially with current concerns over potential climate changes and their effect on the likelihood and height of rogue waves.The focusing nonlinear Schrödinger (NLS) equation often appears in studies of rogue wave formation in deep water when wave amplification is assumed to be primarily due to nonlinear focusing and modulational instability.As a result, several classes of solutions of the NLS equation are considered to be prototypes of rogue waves.For periodic boundary conditions, u(x + L, t) = u(x, t), one such class is the family of homoclinic orbits of unstable plane waves with N unstable modes (Dysthe and Trulsen, 1999;Osborne et al., 2000;Calini and Schober, 2002;Akhmediev et al., 2009a).We will refer to these homoclinic orbits, which can have M ≤ N modes excited, as M mode spatially periodic breather (SPB) solutions (see Figs. 1 and 2).Time-periodic breather-type solutions as well as rational solutions, which arise as singular limits of breather-type solutions and which decay polynomially in space and time, have also been studied (Ankiewicz et al., 2010;Akhmediev et al., 2009b;Ohta and Yang, 2012).
For modeling purposes, the issue of robustness of these families of solutions is important.To successfully observe or reproduce rogue waves in a setting where noise and small higher order nonlinear effects are inherent requires solutions: (i) to remain close to unperturbed ones in the presence of small random variations of initial conditions and (ii) to persist in perturbations of the NLS equation.
In this article we examine the first requirement by investigating the stability with respect to perturbation of initial data of the one-mode SPBs over a plane wave with one or two unstable modes (UMs) and the two-mode SPBs over a plane wave with two UMs.In Sect. 2 we recall the basic elements of the associated Floquet theory, which allows for an exploration of the structure and properties of the SPB solutions.Section 3, the focus of this paper, provides the results of numerical investigations of stability of SPBs with respect to a wide range of initial perturbations f i (x).We consider (i) random shifts in the initial phase, (ii) random spatial perturbations in the height of the wave, (iii) random noise, (iv) localized random Gaussian perturbations, and (v-vi) random high-and low-frequency perturbations.For each type of SPB and for each f i (x) an ensemble of 100 numerical experiments is carried out varying the random component in the initial data.
To study reproducibility/stability numerically, we first find the "closest" element of the SPB family to the perturbed solution.Varying the parameters of the family and using the H 1 norm to measure distances, the closest element is found by minimizing the maximum distance between the perturbed solution and the members of the family of SPBs.Contour plots provide another diagnostic, since they are visually intuitive and show when solutions stay structurally close to each other in "shape".The ensemble estimates of closeness, measured by A(t), indicate that the only neutrally stable SPBs are those for which all the instabilities of the underlying plane wave are saturated (e.g., the two-mode SPB over a plane wave with two UMs).In the numerical simulations the perturbed SPBs may develop a small spatial asymmetry due to the random perturbations.Interestingly, when considering the family of two-mode SPBs, A(t) is smallest for the coalesced two-mode SPB since the spatial asymmetry is minimized.The coalesced SPB was shown (in numerical simulations and by means of perturbation analysis) to be the persistent waveform in various perturbed NLS models on a periodic domain (Calini andSchober, 2002, 2009).This result together with the new observations presented in this article suggests that the coalesced case may be the most robust two-mode SPB in a laboratory setting.Conversely, SPBs that are not fully saturated are sensitive to noisy environments and are unstable.Finally, in Sect. 4 we outline our linear stability analysis of the one-and two-mode SPBs, which support the results of the numerical investigation.

Analytical background
In this section we describe some elements of Floquet spectral theory that are relevant to the stability analysis of the SPBs.The NLS equation is equivalent to the consistency of the Zakharov-Shabat linear system (Z-S) (Zakharov & Shabat, 1972): (2) where λ is the spectral parameter and u(x, t) is a solution of the NLS equation.
For periodic boundary conditions u(x+L, t) = u(x, t), the spectrum of can be described in terms of the Floquet discriminant (u, λ) := tr(M(L; u, λ)), where M(x; u, λ) is the monodromy matrix of 2, as follows: Of particular interest are the following discrete subsets of the periodic spectrum: 1.The simple spectrum, 2. The set of double points, The spectrum of L (x) is invariant under the NLS flow, and each periodic eigenvalue determines the structure and dynamical stability of the corresponding nonlinear mode.In particular, there are no instabilities associated with λ s j or real λ d j , whereas linear instabilities arise when the λ d j 's are complex.
To illustrate the relation between the complex λ d j and the linear instabilities, consider the plane wave solution u a (t) = ae i(2a 2 t+φ) .For small perturbations u(x, t) = u a (t)(1 + (x, t)), | | 1, the quantity is a solution of the linearized equation (5) Thus ∝ e iµ j x+σ j t , where µ j = 2πj/L and σ 2 j = µ 2 j 4|a| 2 − µ 2 j .Then, the plane wave solution is unstable if 0 < (j π/L) 2 < |a| 2 , where the number of unstable modes (UMs) is the largest M such that 0 < M < |a|L/π.On the other hand, one computes the discriminant of the plane wave to be (a; λ) = 2 cos( √ a 2 + λ 2 L), and the discrete spectrum to be λ s 0 = ±ia and (λ d j ) 2 = j π L 2 − a 2 , j = 0. Notice that the λ d j 's are complex if 0 < (j π/L) 2 < |a| 2 , which is the same condition for a mode to be linearly unstable.

SPBs over an unstable plane wave
Explicit representations for the SPBs can be obtained using the Bäcklund-gauge transformation for the NLS equation (see Sect. 4).For an unstable plane wave with N UMs, a single Bäcklund transformation at a complex λ d j generates the one-mode SPB family corresponding to the j th unstable mode, The parameter ρ governs the time at which the mode becomes excited, β j is related to spatial shifts in the solution, µ j = 2πj/L, and p j = arccos πj/aL.The one-mode SPB limits to a phase translation of the plane wave as t → ±∞ with the decay rate σ j .For example, Fig. 1a-b show the amplitudes of the two different one-mode SPBs, U (1) (x, t; ρ) and U (2) (x, t; ρ), over an unstable plane wave with two UMs, for a = 0.5, L = 4 √ 2π, ρ = φ = β = 0, and x ∈ [−L/2, L/2], t ∈ [−10, 10].The one-mode SPB over an unstable plane wave with one UM has the same structure as in Fig. 1a; L is simply adjusted to allow for only one UM.In the next sections we show that the one-mode SPB is neutrally stable and reproducible only when the underlying plane wave has one UM.
Finding a higher dimensional M mode SPB (1 < M ≤ N) requires M iterations of the Bäcklund-gauge transformation Eq. ( 17), where each iteration introduces an additional parameter in the resulting solution.Applying the Bäcklundgauge transformation successively at complex λ d 1 and λ d 2 generates a two-mode SPB family of the form The parameters ρ and τ determine the time at which the first and second mode, respectively, become excited.Ultimately, ρ and τ govern shape, amplitude, and steepness of the SPB, and can be adjusted to excite the modes at the same time.In fact, selecting ρ = −2 and τ = −3 in Eq. ( 7), we obtain what we refer to as the "coalesced" two-mode SPB, whose amplitude is shown in Fig. 2b.Surprisingly, as we will see in the next section, even though the coalesced twomode SPB has steeper gradients, it can be more robust to random perturbations of the initial data than a generic two-mode SPB.

Numerical evidence of stability
To integrate the NLS Eq. ( 1) with periodic boundary conditions, we use a highly accurate and efficient exponential integrator that uses Padé rational-function approximations to the exponential, a Fourier-mode decomposition in space, and a fourth-order Runge-Kutta discretization in time (Khaliq et al., 2009).This scheme has been extensively tested with a variety of known analytical solutions and provides, for refined meshes, sufficient accuracy to simulate solutions of the NLS equation on the time frame under investigation.Chaotic behavior does not develop within the framework of the integrable NLS equation.On a longer time frame, chaotic behavior may develop due to perturbations to the NLS equation arising from the numerical scheme (or, in an experimental setup, from higher order effects).For example, using N = 256 Fourier modes in space and a time step t = 10 −3 , we find that the H 1 norm of the difference between the analytical and numerical solutions is at most O(10 −12 ).On the other hand, the error in the global invariants -the norm, the momentum and the Hamiltonian -is at most O(10 −9 ).
For simplicity, we examine the stability of the one-and two-mode SPB solutions with respect to perturbations in the initial data.The results are generalizable to the case of an M mode breather over an unstable plane wave with N ≥ M unstable modes.We begin by letting or where 0 < 1 should be chosen on the order of experimental error.The parameters ρ and τ are selected so that the difference of U (j ) (x, 0; ρ) and the plane wave is O(10 −3 ), and the difference of U (1,2) (x, 0; ρ, τ ) and the plane wave is O(10 −2 ), in order to avoid exciting any of the instabilities of the plane wave.In all the numerical experiments, the perturbation parameter is = 10 −4 and the time frame is t ∈ [0, 30].There is an inherent limitation to the time frame considered, since eventually the solution will enter a neighborhood of the plane wave and the associated instability becomes manifested due to the numerical error.
We consider the following cases: a one-mode SPB over the plane wave with (i) N = 1 UM or (ii) N = 2 UMs; and (iii) a two-mode SPB over the plane wave with N = 2 UMs.If U (x, t) remains close (in an appropriate sense described below) to an element of the respective family, U (j ) (x, t; ρ) or U (1,2) (x, t; ρ, τ ), then this indicates that the SPB is neutrally stable; otherwise the SPB will be classified as unstable.
In each of the three cases (i-iii) and for each of the initial perturbations f i (x), i = 1, . . ., 6, described below and shown in Fig. 3, an ensemble of 100 numerical experiments was carried out by varying the random component in the initial data: is a spatially random perturbation in the height of the wave. 3.
, where r j (x) ∈ [0, 1] indicates random fields.This represents a set of localized Gaussian perturbations about the points x j . 5. f 5 (x) = K k=−K r k (x)e i2π kx/L for small K, where r k (x) ∈ [0, 1] are random fields.This gives a lowfrequency perturbation.
x)e i2π kx/L for large K, where r k (x) ∈ [0, 1] are random fields.This gives a high-frequency perturbation.Remark.For perturbed initial data, the resulting NLS solution U (x, t) no longer possesses the simple structure of an SPB, as an infinite number of modes become excited.Although over a long time its dynamics may deviate significantly from that of the initially close SPB, numerical investigations of the short-to-moderate-time evolution provide information about the robustness of the SPB within the integrable NLS model, and lay the groundwork for a stability analysis of these solutions.
To study reproducibility/stability numerically, we track the evolution of the norm of the difference of the perturbed solution and the closest element of the unperturbed family.For example, in the case of the one-mode SPB 6, in order to determine the closest element of U (j ) (x, t; ρ) to the perturbed solution, we introduce the quantity compute and then determine the parameter value ρ * , which minimizes As such, U (j ) (x, t; ρ * ) is the closest element, and the evolution of H (j ) (t; ρ * ) provides a measurement of how close the perturbed solution is to an element of the one-mode SPB family.For each f i , we estimate an ensemble measure of "closeness" using the average of H (j ) (t; ρ * ) over all 100 simulations, denoted by A (j ) i (t).(Note that ρ * is different for each simulation.) We also use contour plots as a reproducibility/stability diagnostic tool, since they are visually intuitive and show when solutions remain structurally close to each other in "shape", a feature that cannot be determined by examination of A (j ) i (t) alone.In the contour plots we superimpose the contour of the amplitude obtained from the numerically generated solution U (x, t) onto that of the respective unperturbed analytical solution, U (j ) (x, t; ρ * ) or U (1,2) (x, t; ρ * , τ * ).While only sample contour plots for the different cases are presented, the graphs of A (j ) i (t) provide information obtained from the entire ensemble for each perturbation f i .The numerical results consistently indicate that only the SPBs whose instabilities are saturated are neutrally stable.
Case one.We consider the one-mode SPB over a plane wave with one UM, in particular Eq. ( 6) with j = 1, a = 0.5, ρ 0 = 5.0, and L = 2 √ 2π. Figure 4a shows H (1) mm occurs at ρ * ∼ 5.04.The contours of |U (x, t)| and of |U (1) (x, t; ρ * )|, the nearest one-mode SPB found by minimizing H (1) max (ρ), are given in Fig. 4b.Here, U (x, t) and the nearest SPB are visually identical.Figure 4c shows the evolution of A (1) i (t) for each f i .The small growth in A (1) i (t) to 10 −3 at t ≈ 11 for all f i is due to a small spatial asymmetry that develops in the perturbed solution due to the random nature of the f i .This growth is not significant, as compared, for example, to the growth in A (2) i (t) in Figs. 5 or 6 when the underlying plane wave has two UMs.These results show that the perturbed solution stays close to U (1) (x, t; ρ * ) for a substantial period of time, an indication of the neutral stability of the one-mode SPBs when the underlying plane wave has only one unstable mode.
Case two.Next we consider the one-mode SPB over a plane wave with two UMs, namely formula 6 with j = 1, 2, a = 0.5, ρ 0 = 0.0, and L = 4 and of |U (1) (x, t; ρ * )| are given in Fig. 5a.The closest onemode SPB found by minimizing H (1) max (ρ) matches only the first mode of the perturbed solution.A second mode is excited by the perturbation of the initial data at t ≈ 20, which does not develop in any element of |U (1) (x, t; ρ)|.In fact, small perturbations in the initial data generate quasi-periodic solutions of the NLS equation whose amplitudes resemble a superposition of |U (1)  | and |U (2)  | on this time frame.In Fig. 5b, A (1) i (t), the ensemble measure of closeness, undergoes a rapid growth to O(10) as the second mode develops.This second mode is excited in U (x, t) for all random f i 's, and in fact the maximum of A (1) i (t) is larger for the other perturbations.Figure 6a shows the corresponding contours when U (x, 0) = U (2) (x, 0; ρ 0 ) + f 1 (x) (for k = 1 in f 1 (x)).Similar rapid growth in A (2) i (t) is observed (see Fig. 6b), indicating that the one-mode SPBs over plane waves with N ≥ 2 UMs are unstable.
Case three.Finally, we consider the two-mode SPB over a plane wave with two UMs, given by Eq. ( 7) with i = 1, j = 2, a = 0.5, ρ 0 = −2.0,τ 0 = −10.0,and L = 4 √ 2π.The parameters ρ and τ determine the time when the first and second modes of the SPB become excited.In this case we need to find the element of the family U (1,2) (x, t; ρ, τ ) closest to U (x, t).We find first ρ * and then τ * minimizing the differences between the first and second modes of the perturbed and unperturbed solutions.Namely, to determine the closest element of U (1,2) (x, t; ρ * , τ ) to the perturbed solution, we consider and find the unique τ * , which minimizes H (1,2) max (ρ * , τ ); that is  As before, the ensemble measure of closeness, A (1,2) i (t), is the average of H (1,2) (t; ρ * , τ * ) over all 100 simulations for each f i .

H
(1,2) mm exhibits a larger variance, as can be seen in Table 1, displaying the minimum, mean, median and maximum of H (1,2) mm over the entire ensemble of experiments for each f i .We find H (1,2) mm is at most O 10 −1 (obtained with the random phase f 1 ), with all other f i yielding smaller asymmetries and f 6 , the random high-frequency perturbation, yielding the smallest.One may ask whether the observed spatial asymmetry can be captured explicitly by finding the solutions of Eq. ( 5) since, for random variations in the data, the squared eigenfunctions will no longer be centered around the origin.Since A (1,2) i (t) grows to at most O 10 −1 , U (x, t) stays near U (1,2) (x, t; ρ * , τ * ) for a substantial period of time; that is, the two-mode SPB over a plane wave with two UMs is neutrally stable.
Finally, we consider the special case of the coalesced two-mode SPB over a plane wave with two UMs (recall Fig. 2b).Here U (x, 0) = U (1,2) (x, t; ρ 0 , τ 0 ) + f 5 (x) with a = 0.5, ρ 0 = −2.0,τ 0 = −3.0,and L = 4 √ 2π. Figure 11a shows the contours of |U (x, t)| and of the two-mode SPB 11b) is significantly smaller than in the generic two-mode SPB case (compare with Fig. 10), and is mm is ∼ 0.0091 and ∼ 0.2068, respectively. of the order of A (1) i (t).In this case U (x, t) remains closer to U (1,2) (x, t; ρ * , τ * ), since the coalesced modes appear together earlier in time, and as such U (x, t) is not as susceptible to growth of spatial asymmetries.Vice versa, assuming initial data for an SPB with distinct modes, but with ρ 0 and τ 0 chosen close to the parameter values for the coalesced SPB, it is possible to observe the coalesced SPB due to the shifts in the parameters.
Remarkably, the coalesced two-mode SPB appears to also be more robust under certain types of perturbations of the NLS equation (Calini and Schober, 2002).These two observations indicate that the coalesced case may be the most robust two-mode SPB in a laboratory setting.

Squared eigenfunctions and linear stability
To support the results of the numerical investigation, we outline the linear stability analysis of the one-and two-mode SPB solutions carried out in Calini and Schober (2013).The key observation is that, for a given solution u(x, t) of the NLS equation (e.g., one of the SPBs), its associated "squared eigenfunctions" satisfy the linearized equation about u (i.e., Eq. ( 5) with u a replaced by u).In particular, for a one-mode SPB, if φ and ψ satisfy the Z-S system 2 at U (j ) (x, t), then f (x, t) = φ 1 ψ 1 + φ2 ψ2 and g(x, t) = i(φ 1 ψ 1 − φ2 ψ2 ) solve the linearized NLS equation.Thus, determining stability becomes simply a question of examining the behavior in time of f (x, t) and g(x, t).
The Bäcklund-gauge transformation (Sattinger and Zurkowski, 1987) allows one to transform both the "seed" solution u(x, t) and its eigenfunctions while preserving spatial periodicity, as follows: let φ := α + φ + + α − φ − , α ± ∈ C, where φ + and φ − are linearly independent solutions of the Z-S system at (u, λ j ), with λ j one of the complex λ d j .Construct the following gauge matrix: Then, φ (j ) (x, t, λ; λ j ) = G(λ; λ j , φ)φ(x, t, λ) solves the Z-S system 2 at (U (j ) (x, t), λ), where is the new NLS solution.We use the following notation: the value of superscript j indicates the λ j used in Eq. ((17)), while the number of superscripts is the number of iterations of the Bäcklund-gauge formula.When the seed solution is an unstable plane wave with N unstable modes, for each complex λ d j , the new solution U (j ) (x, t) is the one-mode SPB associated with the j th UM.One iteration of the Bäcklund-gauge transformation produces a two-mode SPB (e.g., U (1,2) (x, t)) as well as its associated eigenfunctions.Since we are interested in the stability of the SPBs, we need the explicit time dependence of the transformed eigenfunctions.A pair of linear independent eigenfunctions of the plane wave is given by where k(λ) = √ λ 2 + a 2 .If the plane wave has only one UM associated with complex λ d 1 , the entries of G(λ; λ j , φ) are bounded in time, since where µ 1 = 2k(λ 1 ) = 2π/L .It follows that the only possible source of exponential-in-time growth of f (x, t) and g(x, t) comes from the eigenfunctions Since χ ± becomes linearly dependent at λ 1 , it suffices to examine χ + (x, t; λ 1 ), which turns out to have no exponential time dependence.In fact, .
Likewise, one finds that the first component is also bounded in time.Therefore, f (x, t) and g(x, t), the solutions of the linearized NLS equation, are bounded in time; in other words, Bäcklund-gauge transformation at λ 1 saturates the associated UM of the plane wave.We conclude that, when the underlying plane wave solution has only one unstable mode, the one-mode SPB is neutrally stable.Similarly, if the plane wave solution has two unstable modes, then applying the Bäcklund transformation successively at λ 1 and λ 2 saturates the associated UMs.Therefore, the two-mode SPBs over a plane wave with two UMs are neutrally stable.
On the other hand, a one-mode SPB over the plane wave with two unstable modes is linearly unstable.In this case, the eigenfunctions χ ± (x, t; λ) of U (1) (x, t; ρ), obtained by implementing the Bäcklund-gauge formula at λ 1 , can be shown to be linearly independent at λ = λ 2 and to exhibit exponential growth in time.In particular, their first component is of the form where B ± (x, t) are bounded and σ = −4iλ 2 k(λ 2 ) is real.In this case f (x, t) and g(x, t) grow exponentially in time; thus U (1) (x, t; ρ) is linearly unstable.These results suggest that only those SPBs for which all UMs are saturated are neutrally stable.

Conclusions
In this article we present the results of extensive numerical investigations of stability of spatially periodic NLS breathers for a variety of random perturbations commonly encountered in experimental settings.It is our hope that such results will provide insight on how to generate reproducible rogue waves in the laboratory.In fact, for a waveform to be reproducible, it should in particular be robust with respect to small perturbations of the initial condition.For the specific case of an SPB over an unstable background with N unstable modes, we observe that only the maximal breathers are stable, in the sense that a small initial perturbation will not grow exponentially in time.This kind of stability should in practice ensure that experimental noise introduced when initializing the wave stays controlled and that the given SPB (or rogue wave) is thus physically realizable or reproducible.More precisely, initializing an experiment with either U (x, 0) = U (1) (x, 0; ρ) or U (x, 0) = U (1,2) (x, 0; ρ, τ ) (for an unstable background with either one or two unstable modes, respectively), and allowing for noise, the generated wave will remain close to an element of the unperturbed family.
A particularly interesting outcome of the numerics is that, among the maximal (and thus stable) two-mode SPBs, the one whose spatial modes have coalesced appears to be the most robust, and therefore it may be the most appropriate candidate for laboratory experiments.Furthermore, in order to facilitate post-processing of the data in a lab setting, our results suggest that an a priori estimate of the shifts in parameters for a given SPB and a prescribed level of noise should be useful in identifying the nearest SPB to the generated wave (i.e., the perturbed solution).Concretely, in every numerical simulation initialized with an order O(10 −4 ) random perturbation of an SPB, the following shifts in parameters were obtained: for one-mode neutrally stable SPBs, the parameter shift is h = ρ * − ρ 0 ≈ O(10 −2 ); and for twomode SPBs, the shifts in the parameters are h = ρ * − ρ 0 ≈ O(10 −2 ) and k = τ * − τ 0 ≈ O(10 −1 ).(These orders of magnitude are consistent with computing h and k by equating a Taylor expansion of U (1) (x, t; ρ 0 + h) with U (1) (x, t) and U (1,2) (x, t; ρ 0 + h, τ 0 + k) with U (1,2) (x, t).)While we utilized the H (1) norm as a measure of closeness, it may be more feasible to compare the maximum amplitudes of the physical and analytical solutions.

Figure 2 .
Figure 2. Amplitude plots of the two-mode SPB over a plane wave with two UMs when the modes are (a) distinct and (b) coalesced.
Figure 5. (a) Contours of |U (x, t)| (dashed line) and the one-mode SPB |U (1) (x, t; ρ * )| (solid line) over a plane wave with two UMs.(b) Evolution of A for each f i .

Figure
Figure 7. H

Table 1 .
The minimum, mean, median and maximum of H