In one of his foundational studies, Stokes established the existence of traveling nonlinear periodic wave trains in deep water [1]. The stability of these waves was resolved when Benjamin and Feir proved that in sufficiently deep water the Stokes wave is modulationally unstable. Small perturbations of Stokes waves were found to lead to exponential growth of the side bands [2, 3]. More recently, modulational instability (MI) of the background state is considered to play a prominent role in the development of rogue waves in oceanic sea states, nonlinear optics, and plasmas [4–9].

The nonlinear Schrödinger (NLS) equation (when ε , γ = 0 in Eq. 1) is one of the simplest models for studying phenomena related to MI; as such, special solutions of the NLS equation are regarded as prototypes of rogue waves. Among the more tractable “rogue wave” solutions of the NLS equation are the rational solutions (with the Peregrine breather being the lowest order) and the spatially periodic breathers (SPBs) which are constructed as heteroclinic orbits of modulationally unstable Stokes waves [10–12]. In the case of the Stokes waves with N unstable modes (UMs), the associated SPBs can be of dimension M ≤ N and are referred to as M-mode SPBs; the single mode SPB is the Akhmediev breather [13]. For more realistic sea states with non-uniform backgrounds, heteroclinic orbits of unstable N-phase solutions have been used to describe rogue waves [14–16].

For theoretical and practical purposes it is important to understand the stability of the SPBs with respect to small variations in initial data and small perturbations of the NLS equation. Using the squared eigenfunction connection between the Floquet spectrum of the NLS equation and the linear stability problem, the SPBs were shown to be typically unstable [17]. The effects of damping on deep water wave dynamics, even when weak, can be significant and in many instances must be included in models to enable accurate predictions of laboratory and field data [18–22].

In a recent study the authors examined the stabilization of symmetric SPBs using the linear damped NLS equation (a near-integrable system that preserves even symmetry of solutions) [23]. The route to stability for these damped SPBs was determined by appealing to the Floquet spectral theory of the NLS equation. Degenerate complex elements of the periodic spectrum (referred to as complex double points) are associated with instabilities of the solution and may split under perturbation to the system. In the restricted subspace of even solutions complex double points can reform as time evolves. The damped solutions were found to be unstable as long as complex double points were present in the spectral decomposition of the data (either by persisting or reforming). A key issue in analyzing the route to stability is determining which complex double points in the spectrum of the SPB are split by damping. For an initial SPB with a given mode structure, perturbation analysis showed that only the complex double points associated with resonant modes split under damping while those associated with nonresonant modes remained closed [23].

The evolution of deep water waves is described only to leading order by the NLS equation. A more accurate description of the wave dynamics is provided by the Dysthe equation, obtained by extending the asymptotic analysis used to derived the NLS equation to fourth order. The Dysthe equation has been shown to accurately predict laboratory data for a wider range of wave parameters than the NLS equation [24, 25, 26]. Gramstad and Trulsen brought the Dysthe equation into Hamiltonian form obtaining a new higher order NLS (HONLS) equation (Eq. 1 with γ = 0 ) [27]. Damped versions of the HONLS equation have successfully described ocean swell and frequency downshift of wave trains on deep water [28, 29].

In this paper we examine the competing effects of dissipation and higher order nonlinearities on the routes to stability of the N-mode SPBs in the framework of the linear damped HONLS equation over a spatially periodic domain: i u t + u x x + 2 | u | 2 u + i ε ( 1 2 u x x x − 8 | u | 2 u x − 2 i u [ ℋ ( | u | 2 ) ] x ) + i γ u = 0 (1) where u ( x , t ) is the complex envelope of the wave train, ℋ { f ( x ) } = 1 π ∫ − ∞ ∞ f ( ξ ) x − ξ d ξ is the Hilbert transform of f, and 0 < ε , γ ≪ 1 . The initial data used in the numerical experiments is generated using exact SPB solutions of the integrable NLS equation. The SPBs are over Stokes waves with N unstable modes (referred to as the N-UM regime) for 1 ≤ N ≤ 3 . We interpret the damped HONLS (near-integrable) dynamics by appealing to the NLS Floquet spectral theory.

The higher order nonlinearities in Eq. 1 break the even symmetry of both the initial data and the equation. This raises several interesting questions regarding the damped HONLS equation. Which integrable instabilities are excited by the damped HONLS flow and which elements of the Floquet spectrum are associated with these instabilities? What are the routes to stability under damping; i.e., what remnants of integrable NLS structures are detected in the damped HONLS evolution?

In the present study we observe the onset of novel instabilities as a result of symmetry breaking and the development of critical states in the damped HONLS flow which were nonexistent in the previously examined damped NLS system with even symmetry. Significantly, we determine these instabilities are associated with degenerate complex elements of both the periodic and continuous spectrum, i.e., with both complex “double points” and complex “critical points”, respectively. This association was not previously recognized. With regard to teminology, although double points are among the critical points of Δ , in this paper we exclusively call degenerate complex periodic spectrum where Δ = ± 2 “double points” and reserve the term “critical points” for degenerate complex spectrum where Δ ≠ ± 2 .

The paper is organized as follows. In Section 2 we present elements of the NLS Floquet spectral theory which we use to distinguish instabilities in the numerical experiments. Whether higher phase solutions, such as the even 3-phase solution given in Eq. 8, are unstable with respect to general noneven perturbations and what the Floquet “signature” is of the possible instabilities, has been an open question. The closest stability results we are aware of are for the elliptic solutions of the focusing NLS equation [30]. We numerically show an even 3-phase solution of the NLS equation is unstable with respect to generic perturbations of initial data and find the relevant element of the Floquet spectrum associated with the instability in order to develop a broadened Floquet characterization of instabilities of the NLS equation.

A brief overview of the SPB solutions of the NLS is provided at the end of Section 2 before numerically examining their stabilization under the damped HONLS flow in Section 3. The Floquet decompositions of the numerical solutions are computed for 0 ≤ t ≤ 100 . Complex double points are initially present in the spectrum. If one of the complex double points present initially splits due to the damped HONLS perturbation, the subsequent evolution involves repeated formation and splitting of complex critical points (not double points) which we correlate with the observed instabilities. The Floquet spectral analysis is complemented by an examination of the growth of small perturbations in the SPB initial data under the damped HONLS flow. We determine that the instabilities saturate and the solutions stabilize once all complex double points and complex critical points vanish in the spectral decomposition of the perturbed flow. Variations in the spectrum under the HONLS flow are correlated with deformations of certain NLS solutions to determine the routes to stability for the damped HONLS SPBs.

In Section 4, via perturbation analysis, we examine splitting of the complex double points, present in the SPB initial data, under the damped HONLS flow. We find that for short time, the complex double points associated with modes that resonate with the SPB structure split producing disjoint asymmetric bands, while the complex double points associated with nonresonant modes remain closed, substantiating the initial spectral evolutions observed in the numerical experiments. The nonresonant double points are observed to remain closed for the duration of the experiments, beyond the time-frame of the short time analysis, even though the solution evolves as a damped asymmetric multi-phase state. In this study resonances have a stabilizing effect; the instabilities of nonresonant modes persist on a longer time scale than the instabilities associated with resonant modes.

The nonlinear Schrödinger equation (when ε , γ = 0 in Eq. 1) arises as the solvability condition of the Zakharov-Shabat (Z-S) pair of linear systems [31]: ℒ ( u ) ϕ = λ ϕ ,     ℒ ( u ) = ( i ∂ x u − u ∗ − i ∂ x ) (2) ϕ t = ( − 2 i λ 2 + i | u | 2 2 i λ u − u x 2 i λ u ∗ + u x ∗ 2 i λ 2 − i | u | 2 ) ϕ (3) where λ is the spectral parameter, ϕ is a complex vector valued eigenfunction, and u ( x , t ) is a solution of the NLS equation itself. Associated with an L − periodic NLS solution is it’s Floquet spectrum σ ( u ) : = { λ ∈ ℂ | ℒ ϕ = λ ϕ , | ϕ |   b o u n d e d   ∀ x } (4)

Given a fundamental matrix solution of the Z-S system, Φ , one defines the Floquet discriminant as the trace of the transfer matrix across one period L, Δ ( u , λ ) = T r a c e [ Φ ( x + L , t ; λ ) Φ − 1 ( x , t ; λ ) ] . The Floquet spectrum has an explicit representation in terms of the discriminant: σ ( u ) : = { λ ∈ ℂ | Δ ( u , λ ) ∈ ℝ , − 2 ≤ Δ ( u , λ ) ≤ 2 } (5)

The Floquet discriminant Δ ( λ ) is analytic and is a conserved functional of the NLS equation. As such, the spectrum σ ( u ) of an NLS solution is invariant under the time evolution.

The spectrum consists of the entire real axis and curves or “bands of spectrum” in the complex λ plane ( ℒ ( u ) is not self-adjoint). The periodic/antiperiodic points (abbreviated here as periodic points) of the Floquet spectrum are those at which Δ = ± 2 . The endpoints of the bands of spectrum are given by the simple points of the periodic spectrum σ s ( u ) = { λ j s | Δ ( λ j ) = ± 2 , ∂ Δ / ∂ λ ≠ 0 } . Located within the bands of spectrum are two important spectral elements:

1. Critical points of spectrum, λ j c , determined by the condition ∂ Δ / ∂ λ = 0 .

Double points are among the critical points of Δ . However, in this paper, we exclusively call the degenerate periodic spectrum where Δ = ± 2 “double points” and reserve the term “critical points” for degenerate elements of the spectrum where Δ ≠ ± 2 .

The Floquet spectrum can be used to represent a solution in terms of a set of nonlinear modes where the structure and stability of the modes are determined by the band-gap structure of the spectrum. Simple periodic points are associated with stable active modes. The location of the double points is particularly important. Real double points correspond to zero amplitude inactive nonlinear modes. On the other hand, complex double points are associated with degenerate, potentially unstable, nonlinear modes with either positive or zero growth rate. When restricted to the subpace of even solutions, exponential instabilities of a solution are associated with complex double points in the spectrum [32].

A concrete example illustrating the correspondence between complex double points in the spectrum and linear instabilities is the Stokes wave solution u a ( t ) = a e i ( 2 a 2 t + ϕ ) . For small perturbations of the form u ( x , t ) = u a ( t ) [ 1 + ε ( x , t ) ] , | ε | ≪ 1 , one finds ϵ satisfies the linearized NLS equation i ε t + ε x x + 2 | a | 2 ( ε + ε ∗ ) = 0 (6)

Representing ϵ as a Fourier series with modes ε j ∝ e i μ j x + σ j t , μ j = 2 π j / L , gives σ j 2 = μ j 2 ( 4 | a | 2 − μ j 2 ) . As a result, the jth mode is unstable if 0 < ( j π / L ) 2 < | a | 2 . The number of UMs is the largest integer M such that 0 < M < | a | L / π .

The Floquet discriminant for the Stokes wave is Δ = 2 ⁡ cos ( a 2 + λ 2 L ) . The Floquet spectrum consists of continuous bands ℝ ∪ [ − i a , i a ] and a discrete part containing λ 0 s = ± i | a | and the infinte number of double points ( λ j d ) 2 = ( j π L ) 2 − a 2 ,   j ∈ ℤ ,   j ≠ 0 (7) as shown in Figure 1. Note that the condition for λ j d to be complex is precisely the condition for the jth Fourier mode ε j to be unstable. The remaining λ j d for | j | > M are real double points.

As the NLS spectrum is symmetric under complex conjugation, we subsequently only display the spectrum in the upper half λ-plane.

In our examination of the even 3-phase solution in Section 2 we found that when evenness is relaxed novel instabilities arise which are associated with complex critical points. Armed with this result, in this section we return to the questions posed at the outset of our study: 1) Which integrable instabilities are excited by the damped HONLS flow and what is the Floquet criteria for their saturation? 2) What remnants of integrable NLS structures are detected in the damped HONLS evolution?

The notation used in this section is as follows: 1) The “N UM regime” refers to the range of parameters a and L for which the underlying Stokes wave initially has N unstable modes. 2) The initial data used in the numerical experiments is generated using exact SPB solutions of the integrable NLS equation. The perturbed SPBs are indicated with subscripts: U ε , γ ( j ) ( x , t ) refers to the solution of the damped HONLS Eq. 1 for one-mode SPB initial data U ( j ) ( x , 0 ) . Likewise U ε , γ ( i , j ) ( x , t ) is the solution to Eq. 1 for iterated SPB initial data.

The damped HONLS equation is solved numerically using a high-order spectral method due to Trefethen [37]. The integrator uses a Fourier-mode decomposition in space with a fourth-order Runge-Kutta discretization in time. The number of Fourier modes and the time step used depends on the complexity of the solution. For example, for initial data in the three UM regime, N = 1024 Fourier modes are used with time step Δ t = 7.5 × 10 − 5 . As a benchmark the first three global invariants of the HONLS equation, the energy E = ∫ 0 L | u | 2 d x , momentum P = i ∫ 0 L ( u ∗ u x − u u x ∗ ) d x , and Hamiltonian H = ∫ 0 L { − i | u x | 2 + i | u | 4 − ε 4 ( u x u x x ∗ − u x ∗ u x x ) + 2 ε | u | 2 ( u ∗ u x − u u x ∗ ) + i ε | u | 2 [ ℋ ( | u | 2 ) ] x } d x are preserved with an accuracy of O ( 10 − 12 ) for 0 ≤ t ≤ 100 . The invariant for the damped HONLS system, the spectral center k m = − P / 2 E , is preserved with an accuracy of at least O ( 10 − 12 ) in the experiments.

Nonlinear mode decomposition of the damped HONLS flow: At each time t we compute the spectral decomposition of the damped HONLS data using the numerical procedure developed by Overman et. al. [34]. After solving system Eq. (2), the discriminant Δ is constructed. The zeros of Δ ± 2 are determined using a root solver based on Müller’s method and then the curves of spectrum filled in. The spectrum is calculated with an accuracy of O ( 10 − 6 ) which is sufficient given the perturbation parameters ϵ and γ used in the numerical experiments are O ( 10 − 2 ) .

Notation used in the spectral plots: The periodic spectrum is indicated with a large × when Δ = − 2 and a large box when Δ = 2 . The continuous spectrum is indicated with small × when the Δ is negative and a small box when Δ is positive.

Interpreting the damped HONLS flow via the NLS spectral theory: A tractable example which illustrates the use of the Floquet spectrum to interpret near integrable dynamics is the spatially uniform solution (there is no depenence on ϵ) of the damped HONLS Eq. 1, u a , γ ( t ) = a e − γ t e i ( | a | 2 ( 1 − e − 2 γ t ) γ ) (12)

At a given time t = t ∗ the nonlinear spectral decomposition of Eq. 12 can be explicitly determined by substituting u a , γ ( t ∗ ) into ℒ ( u ) ϕ = λ ϕ . We find the periodic Floquet spectrum consists of λ 0 s = ± i a e − γ t ∗ and infinitely many double points ( λ j d ) 2 = ( j π L ) 2 − e − 2 γ t ∗ a 2 ,   j ∈ ℤ ,   j ≠ 0 (13) where λ j d is complex if ( j π L ) 2 < e − 2 γ t ∗ a 2 . Under the damped HONLS evolution the endpoint of the band of spectrum λ 0 s and the complex double points λ j d move down the imaginary axis and then onto the real axis. Similarly, at t = t ∗ , a linearized stability analysis about u a , γ ( t ∗ ) shows the growth rate of the jth mode σ j 2 = μ j 2 ( 4 e − 2 γ t ∗ a 2 − μ j 2 ) . Thus the number of complex double points and the number of unstable modes diminishes in time due to damping. As a result, u a , γ ( t ) stabilizes when the growth rate σ 1 = 0 , i.e., when λ 1 d = 0 , giving t ∗ = ln ( a L / π ) / γ .

Saturation time of the instabilities: Since the association of complex critical points with instabilities is a new result, we supplement the spectral analysis with an examination of the saturation time of the instabilities for the damped SPBs U ε , γ ( j ) ( x , t ) and U ε , γ ( j , k ) ( x , t ) as follows: we examine the growth of small asymmetric perturbations in the initial data of the following form, U ε , γ , δ ( j ) ( x , 0 ) = U ( j ) ( x , 0 ) + δ f k ( x ) and U ε , γ , δ ( i , j ) ( x , 0 ) = U ( i , j ) ( x , 0 ) + δ f k ( x ) where.

i. f k ( x ) = cos ⁡ μ k x + r k e i ϕ k ⁡ sin ⁡ μ k x , μ k = 2 π k / L , 1 ≤ k ≤ 3 ,

To determine the closeness of U ε , γ ( j ) ( x , t ) and U ε , γ , δ ( j ) ( x , t ) as time evolves we monitor the evolution of η ( t ) , as given by Eq. 9. We consider the solution to have stabilized under the damped HONLS flow once η ( t ) saturates.

In the damped HONLS numerical experiments we obtain a new criteria for the saturation of instabilities: η ( t ) saturates and the SPB stabilizes once damping eliminates all complex double points and complex critical points in the spectrum.

While examining the route to stability of the SPBs under the damped HONLS several novel results arose. One feature was that the instabilities of nonresonant modes persist longer than the instabilities of the resonant modes. We are interested in the fate of complex double points under noneven perturbations induced by HONLS as they characterize the SPB. Following the perturbation analysis in [39] used to determine the O ( ε ) splitting of double points for single mode perturbations, we carry the analysis to higher order for noneven multi mode perturbations of the SPBs. We find 1) additional modes resonate with the perturbation and 2) complex double points associated with nonresonant modes remain closed.

To obtain linearized initial conditions for the one and two mode SPBs we use the Hirota formulation of the SPBs [40]. For example for the one mode SPB one obtains, u ( j ) ( x , t ) = a e 2 i a 2 t 1 + 2 e 2 i θ j + Ω j t + γ ⁡ cos ⁡ μ j x + A 12 e 2 ( 2 i θ j + Ω j t + γ ) 1 + 2 e Ω j t + γ ⁡ cos ⁡ μ j x + A 12 e 2 ( Ω j t + γ ) (14) where μ j = 2 π j / L , Ω j = μ j 4 a 2 − μ j 2 , sin ⁡ θ j = μ j / 2 a , A 12 = sec 2 θ j , and γ is an arbitrary phase.

The appropriate linearized initial conditions for the one and two mode SPBs, u ( i ) ( x , 0 ) and u ( i , j ) ( x , 0 ) respectively, are obtained by choosing t and γ such that ε ˜ s = 4 i ⁡ sin ⁡ θ s   e Ω s t + γ , s = i , j , are small. After neglecting second-order terms we obtain: u ( j ) ( x , 0 ) = a ( 1 + ε ˜ j e i θ j ⁡ cos ⁡ μ j x ) (15) u ( i , j ) ( x , 0 ) = a ( 1 + ε ˜ i   e i θ i ⁡ cos ⁡ μ i x + ε ˜ j   e i θ j ⁡ cos ⁡ μ j x ) (16)

The damped HONLS yields the following noneven first order approximation for small time, u ( j ) ( x , h ) = a [ 1 + ε ˜ j ( e i θ j ⁡ cos ⁡ μ j x + r j e i ϕ j ⁡ sin ⁡ μ j x ) ] (17) u ( i , j ) ( x , h ) = a [ 1 + ε ˜ i ( e i θ i ⁡ cos ⁡ μ i x + r i   e i ϕ i ⁡ sin ⁡ μ i x ) + ε ˜ j ( e i θ j ⁡ cos ⁡ μ j x + r j e i ϕ j ⁡ sin ⁡ μ j x ) ] (18) where θ s ≠ ϕ s and a , ε ˜ s , r s are functions of h and the damped HONLS parameters ϵ and γ, for s = i , j . For simplicity we set ε = ε ˜ s and suppress their explicit dependence on ε , γ : u = a + ε [ e i θ i ⁡ cos ⁡ μ i x + r i e i ϕ i ⁡ sin ⁡ μ i x + Q ( e i θ j ⁡ cos ⁡ μ j x + r j e i ϕ j ⁡ sin ⁡ μ j x ) ] = a + ε u ( 1 ) (19) where r s ≠ 0 and Q can be 0 or 1, depending on whether a one or two mode SPB is under consideration.

Since Δ ( λ , u ) and the eigenfunctions v n = [ v n 1 v n 2 ] are analytic functions of their arguments, at the double points λ n we assume the following expansions: v n = v n ( 0 ) + ε v n ( 1 ) + ε 2 v n ( 2 ) + ⋯ (20) λ n = λ n ( 0 ) + ε λ n ( 1 ) + ε 2 λ n ( 2 ) + ⋯ (21) Substituting these expansions into Eq. 25 we obtain the following: O ( ε 0 ) :   ℒ v n ( 0 ) = 0 (22) O ( ε 1 ) :   ℒ v n ( 1 ) = [ − i λ n ( 1 ) v n 1 ( 0 ) + u ( 1 ) v n 2 ( 0 ) − i λ n ( 1 ) v n 2 ( 0 ) + u ( 1 ) ∗ v n 1 ( 0 ) ] ≡ F (23) O ( ε 2 ) :   ℒ v n ( 2 ) = [ − i λ n ( 1 ) v n 1 ( 1 ) − i λ n ( 2 ) v n 1 ( 0 ) + u ( 1 ) v n 2 ( 1 ) − i λ n ( 1 ) v n 2 ( 1 ) − i λ n ( 2 ) v n 2 ( 0 ) + u ( 1 ) ∗ v n 1 ( 1 ) ] ≡ G (24) where ℒ = [ ∂ / ∂ x + i λ n ( 0 ) − a − a − ∂ / ∂ x + i λ n ( 0 ) ] (25)