Complex Solitary Wave Dynamics, Pattern Formation, and Chaos in the Gain-Loss Nonlinear Schr\"odinger Equation

A numerical exploration of a gain-loss nonlinear Schr\"odinger equation was carried out utilizing over 180000 core hours to conduct more than 10000 unique simulations in an effort to characterize the model's six dimensional parameter space. The study treated the problem in full generality, spanning a minimum of eight orders of magnitude for each of three linear and nonlinear gain terms and five orders of magnitude for higher order nonlinearities. The gain-loss nonlinear Schr\"odinger equation is of interest as a model for spin wave envelopes in magnetic thin film active feedback rings and analogous driven damped nonlinear physical systems. Bright soliton trains were spontaneously driven out of equilibrium and behaviors stable for tens of thousands of round trip times were numerically identified. Nine distinct complex dynamical behaviors with lifetimes on the order of ms were isolated as part of six identified solution classes. Numerically located dynamical behaviors include: (1) Low dimensional chaotic modulations of bright soliton trains; (2) spatially symmetric/asymmetric interactions of solitary wave peaks; (3) dynamical pattern formation and recurrence; (4) steady state solutions and (5) intermittency. Simulations exhibiting chaotically modulating bright soliton trains were found to qualitatively match previous experimental observations. Ten new dynamical behaviors, eight demonstrating long lifetimes, are predicted to be observable in future experiments.


Introduction
Related forms of the nonlinear Schrödinger equation (NLS) are used to explore nonlinear phenomena in many distinct physical systems: Ginzburg-Landau equations describe the envelope evolution of mode-locked lasers and superconductivity [1]; the cubic NLS treats deep water waves [2] and the dynamics of spin-wave envelopes in magnetic thin films [3,4]; a driven damped NLS models exciton-polariton and magnon Bose-Einstein condensates (BECs) [5]; and the Gross-Pitaevskii equation models the mean field of atomic and molecular BECs [6,7]. With the increasing supply of cheap computing power these systems have become the subject of extensive and sometimes rigorous numerical study.
A majority of these phenomena were observed on active feedback rings; such feedback structures are ubiquitous within science and physics in general. Rings in particular are commonly used to study wave dynamics when one seeks a quantized wavenumber, periodic pumping, self-generation or other resonant phenomena. Active rings, so called for the presence of periodic linear amplification, allow for the direct compensation of the major loss mechanisms present within a system. This permits one to drive the system into quasi-conservative regimes, enabling the observation of dynamics on scales of several to tens of thousands of round trip times. Dynamics with lifetimes of this order would otherwise be prohibited by the presence of dissipation.
Yet, within the context of spin waves in magnetic thin films, little work has been carried out to develop an adequate theory for describing the rich range of behaviors evident within these recent experimental works. The integrable cubic NLS, while successful in quantitatively describing both dark and bright soliton trains, is unable to reproduce more complex phenomena such as the chaotic oscillation of soliton envelopes. However, there has been significant effort within the mode-locked laser community to study analogous driven and damped systems. Works on dissipation terms and saturation [38][39][40][41], the study of dissipative solitons dynamics [42][43][44][45][46] and other numerical studies of the cubic quintic complex Ginzburg-Landau equation [47,48] are highly relevant to the development of a driven damped model for spin waves in magnetic thin films. These works explore the dynamics of solitary waves and their associated wave equations under the influences of gain and loss. For example, the dynamics of near steadystate dissipative solitons have been considered in detail; such studies include rigorous mappings of stable and unstable regions of parameter space [44][45][46]. There have also been significant investigations made into complex soliton dynamics in optical lattices [49,50]. Similarly, initial transient behaviors have been the subject of significant research efforts by the mode-locked laser community. Transients are of interest for potential applications in signal processing and communication. To date there have been no efforts toward the characterization of long lifetime (>1 ms) dynamics of soliton trains driven from equilibrium within active feedback rings. The work presented here demonstrates such an effort for a generalized NLS with a focus on applications to nonlinear spin waves propagating in an active feedback ring.
Our paper is outlined as follows. Section 2 introduces the model to be studied along with the associated operating limits; here the methodology and scope of the simulations are explicitly defined. Experimental contexts for the work are also considered in section 2. Results in the form of 11 unique complex dynamical behaviors are presented and categorized in sections 3-7. Section 3 contains simulations of chaotically modulating soliton trains. Spatially symmetric and asymmetric solitary wave interactions are presented in section 4. Four examples of dynamical pattern formation are given in section 5. Two cases of steady state evolution are reported in section 6. Intermittent solutions are discussed in section 7. Animations of the dynamical behaviours presented in sections 3-7 are available in the supplementary data. Finally, a discussion of numerical convergence and solution robustness is given in section 8. The work is summarized in section 9.

Model overview
To study the long lifetime dynamics of soliton trains driven from equilibrium, motivated by the works discussed in section 1, we propose a generalized governing equation for spin waves in magnetic thin film active feedback rings: the gain-loss NLS (GLNLS): ∂ 2 ∂ x 2 + iL + (N + iC)|u| 2 + (S + iQ)|u| 4 where u = u(x, t) is a dimensionless complex magnetization amplitude defined as |u(x, t)| 2 = m(x, t) 2 /2M s 2 ; here m(x, t) is the dynamic magnetization while M 2 s is the saturation magnetization; D is the dispersion coefficient; N and S are the cubic and quintic nonlinearity coefficients, respectively; t is the 'temporal' evolution coordinate; x is the 'spatial' coordinate of propagation boosted to the group velocity of the envelope; and L, C and Q are the linear, cubic and quintic gains (losses) if positive (negative). All parameters are taken to be real as the complex nature of the coefficients is explicitly accounted for in (1). The local intensity of the magnetization amplitude is given by |u(x, t)| 2 . The norm and energy at a given time, t, are defined as and respectively, where the integrals are taken over the length or circumference, 2π R, of the feedback ring. All norms, intensities and energies given within figures and animations are scaled by u(0) 2 , max[|u(x, 0)| 2 ] and abs[E(0)], respectively, where t = 0 corresponds to the initial condition used during numerical simulation. Numerical values given within the text are not scaled. The specific choice of initial condition is discussed later in this section. The gain and loss present within the GLNLS may be viewed as an expansion of saturable loss expressions studied separately by Ablowitz et al [39], Ablowitz and Horikis [40], Akhmediev et al [44] and Akhmediev and Ankiewicz [45]. The higher order nonlinearity, S, may be used either as a saturation of cubic nonlinearity or an additional self-steepening; both cases are studied in the literature [2]. The GLNLS omits other terms commonly included in cubic quintic complex Ginzburg-Landau equations such as spectral filtering, periodic pumping and integral mean terms, as they are not needed in this physical context. We are likewise compelled by Occam's razor to choose the simplest possible model which nevertheless reproduces measurements in magnetic thin film active feedback rings. NLS-like equations may be derived in magnetic thin films by use of a slowly varying envelope approximation [51], more rigorously through conservation considerations and a Hamiltonian formalism [9,52], or directly from Maxwell's equations using multi-scale methods [53]. The operating limits of the GLNLS are motivated principally by experimental work on the excitation of chaotic solitons in YIG strip-based active feedback rings [36]. A block diagram of the active feedback ring experiment is shown in figure 1. The ring is comprised of a nonlinear propagation medium, in this case a magnetically saturated crystalline YIG thin film, connected via two transducers to an electronic feedback loop. The electronics loop is constructed of a directional coupler, allowing real time observation at an oscilloscope and/or spectrum analyzer, and an amplifier/attenuator pair for real time adjustment of ring gain. The GLNLS demonstrated qualitative agreement with the low dimensional chaotic modulation of a bright soliton train, as will be discussed in section 3. Detailed experimental results are discussed in Wang et al [36]. These experiments indicate that nonlinearity and dispersion are the dominant sources of envelope shaping for chaotic spin wave solitons and that the losses present in the ring are fully compensated for by the amplifier. This imposed two constraints on modeling. (i) The coefficients N and D must be orders of magnitude larger than L, C and Q. (ii) The linear amplifier must compensate both the linear and nonlinear losses present in the film, requiring a net averaged (over many round trips) linear gain, L > 0. Likewise, the dissipative terms represent the net gain and loss processes occurring in the ring averaged over several round trip times. One expects the use of this approximation, and therefore the model, to be valid when the time scale of envelope modulation is much greater than the soliton round trip time.
All simulations of the GLNLS were performed using adaptive time step Cash-Karp Runge-Kutta for temporal evolution and pseudospectral techniques for spatial propagation [54,55]. Periodic boundary conditions modeled propagation around a ring. Note, a detailed discussion on numerical convergence for complex dynamics may be found in section 8. Every simulation began as a bright soliton initial condition obtained via imaginary time relaxation [56], the ground state solution to the GLNLS with S, L, C and Q set as zero, with |u(x, t)| 2 < 1. Experimentally this corresponds to a stable bright soliton circling within a YIG strip-based active feedback ring, a solution analogous to a soliton train. Dynamical results were generated by driving this bright soliton initial condition out of equilibrium via numerical propagation under the influences of non-zero gain and loss terms. This is a process analogous to active feedback ring experiments where gains are increased beyond those which support stable bright soliton trains [19,28,30,36,57]. As previously motivated in section 1 we are interested in identifying dynamical behaviors which are stable on timescales observable in magnetic thin film active feedback ring experiments, (>1 ms) or 7000+ round trips, and analogous physical systems. For this reason any transient dynamics which occurred as a result of driving the bright soliton train out of equilibrium were not explicitly studied within this work, so long as such dynamics were numerically converged. Finally, a rigorous study of transient dynamics is more appropriately done via an iterative simulation method, rather than the slowly varying envelope approximation used here.
Physical GLNLS parameter values are obtained by fitting this initial condition to experimentally observed bright soliton train conditions. This choice of units also fixes the ratio of N /D used in simulations, while the amplitude of N and D dictate the simulation timescale. We assumed that the dimensionless spin-wave intensity is directly proportional to the spin-wave power, |u(x, t)| 2 ∝ P out , since experimental measurements of voltage are taken across a diode with quadratic behavior and are generally taken to be proportional to power. Values typical for a chaotic soliton experiment are T = 165 ns, the round trip time; d = 0.55 cm, the transducer separation; T e = 10 ns, the electronic loop propagation time; V g = d/(T − T e ) = 3.5 × 10 6 cm s −1 , the group velocity; N = −9.24 × 10 9 rad s −1 , the cubic nonlinearity; and D = 510 cm 2 s −1 , the dispersion. Using these parameters one finds [t] ≈ 25 ns where t is the scaled temporal unit used in simulations. This relation may be used to immediately transform code values for L, C, Q and S, which share units of inverse time, to physical values. For example the largest studied linear gain is L = t −1 ≈ 0.05 ns −1 which matches the order of experimentally approximated linear losses for magnetostatic backward volume spin waves in YIG thin films [57]. Experimentally a time series is recorded at the detection transducer with the full waveform being captured once a round trip after the signal has propagated a length d between the transducers and passed through the electronics loop. The length of the ring, , is taken to be the transducer separation, d, as the propagation delay is orders of magnitude smaller than the round trip time, T e < T . Simulations explicitly model the entire feedback loop at the group velocity of the waveform. A time series may be reconstructed from numerical data by concatenating the simulated waveform after a temporal evolution of T or a spatial evolution of d = . In this work we adopt the former convention to ease the direct comparison of simulations to the power versus time data often observed experimentally for spin waves in magnetic thin films. Such a reconstructed time series is labeled u ts (t) throughout the paper. A time series of solitary wave peak intensity at successive round trips is useful in studying modulating single solitary wave trains and is defined by where T is the round trip time and N RT is the total number of round trips. Parameter space explorations were explicitly chosen to encompass the GLNLS operating regime for magnetic thin film systems, while extending into other limits that could be of interest to other systems where the GLNLS is a useful model. Along with the previously mentioned restriction on the sign of L only cases with cubic losses, C 0, were considered. Both instances of saturating 4 quintic gains, Q 0, and supplemental quintic losses, Q 0, were studied. No sign restrictions were placed on quintic nonlinearity, S. The terms were explored in a decadal fashion across the GLNLS scaled values listed here: Ignoring cases with solely gains present we performed 5470 unique simulations. An additional 1530 simulations were undertaken with random parameters. The value for any single parameter in these simulations was generated by multiplying a pseudo random number between zero and one, from the uniform distribution, by an order of magnitude and sign chosen at random, again with uniform weight, from a parameter's allowed values, as defined above. To avoid ambiguity all statements in this paper concerning the relative size of GLNLS parameters refer to the order of magnitude and not the sign. Over 180 000 core hours were utilized to conduct more than 10 000 unique simulations and convergence studies. An initial study of 3500 simulations was undertaken to explore the extent of transient effects and the numerical convergence behaviors of the GLNLS. A summary and analysis of the subsequent 7000 simulations, corresponding to over 3 TB of data, are presented in sections 3-7. Approximately 1500 simulations were evaluated in detail; the remaining simulations were spot checked for consistency. Dozens of complex dynamical behaviors were identified during the course of simulation. We call this system complex because it displays a rich variety of dynamical behaviors, including chaos, robust emergent solitary-wave features and generally multiple scales in both space and time. Solution types were divided into three stability cases, with each case corresponding to roughly 30% of observed dynamics. The three cases are temporally stable, intermittent and unstable. Temporally stable solutions demonstrated substantial observable lifetimes, greater than 1 ms or 7000+ round trips, and robustness to variations in initial conditions of at least 10%. Evolution was found to be least sensitive to changes in S and Q and most sensitive to perturbations in L. In general the effect of changes in initial conditions tended to degrade the lifetime of dynamical behaviors and push solutions toward the intermittent case. Nine temporally stable distinct dynamics and two separate cases of intermittency are discussed below. A summary of the GLNLS parameter regimes which support these identified dynamics is given in table 1. 4 A typical expression for saturable gain is given by iS g (1 + |u(x,t)| 2 I s ) −1 where I s and S g are control parameters.
Expanding denominator to third order yields iS g (1 − |u(x, t)| 2 I −1 s + |u(x, t)| 4 I −2 s + · · · ), hence positive quintic gain being named saturating. Table 1. Overview of identified long lifetime dynamical behaviors solution types and the range of GLNLS parameters which support them. The signs of quintic nonlinearity, S, and gain, Q, values are identified explicitly while the sign of linear gain is taken to be positive and cubic gain is taken to be negative as discussed in section 2. The order of the quintic nonlinearity, S, is often directly compared to that of the cubic nonlinearity, N . The statement O(x) ± y is meant to be read order of x plus or minus y orders of magnitude. For example the symmetric interaction dynamic is observed when the cubic loss, C, is within plus or minus one order of magnitude of the linear gain term, L.

Chaotic modulation
The chaotic modulation of stable solitary wave trains was observed for solutions containing strongly saturated cubic nonlinearity, S 10 −2 , and the lowest studied ring gains, L = 10 −7 , with matching orders of cubic and/or quintic losses. A single bright soliton is observed to circulate within an active feedback ring while exhibiting complex modulations in peak intensity. Low ring losses are anticipated for this solution type, as experimentally observed chaotically modulating soliton trains have lifetimes measured in seconds. The presence of a single stable bright soliton suggests that nonlinearity and dispersion are the dominate forces in peak shaping. These are two conditions used during the derivation of the GLNLS, (1), as discussed previously in section 2. The chaotic nature of measured time series was verified by using standard phase space reconstruction techniques available in the open source nonlinear time series (TISEAN) package to arrive at a stable correlation dimension, D 2 [58]. The correlation dimension, a phase space invariant, was estimated via computation of the correlation sum for increasing embedding dimensions of the time series [59]. The standard embedding procedure of Taken and Sauer was followed using time-delayed reconstruction of the time series [60,61]. The time delay was chosen as the first minimum of autocorrelation to maximize the linear independence of the time delayed vectors. As the phase space was reconstructed from a single time series, a Theiler window of ten times the single round trip time was used to avoid the misinterpretation of temporal correlation as geometrical structure on the attractor [62]. If the correlation dimension was observed to saturate with increasing embedding dimension the time series was said to have a stable correlation dimension. If the stable correlation dimension was not an integer then the system was said to be chaotic.
We further required the correlation dimension to be stable across a wide range of embedding parameters as one expects the reconstructed attractor to be invariant under smooth transformations. This requirement was extremely conservative as it was computationally onerous and sensitive to noise. However, such a requirement forbids the optimization of phase space invariants by the tuning of embedding parameters, and the requirement of saturation across embedding dimension eliminates any assumptions required to study a single reconstruction. Additional indicators of chaos include broadband spectra and positive Lyapunov exponents [59]; note both these properties are shared with noise so a finite correlation dimension is necessary to demonstrate chaotic, rather than random, motion. The principle challenge to finding a stable correlation dimension was isolating a stationary solution.
Two examples, with peak variations of 2.0 and 5.1% about their mean, are shown in figures 2 and 3, respectively. Per cent peak variation is defined as where u peak was previously defined in equation (4) and var is the sample variance. These values are chosen to match the peak variation of two low ring gain chaotic solitary wave trains observed experimentally by Wang et al [36]. In both figures panel (a) shows the intensity of the single soliton peak for 2500 consecutive round trips while panel (b) shows 100 round trips as would be observed experimentally, as in figure 4(d) below, and each vertical line is in fact a bright soliton of finite width. The single soliton peak intensity is immensely complex on inspection and is at worst random and at best chaotic or quasi-periodic. Animations of the data shown in figures 2(b) and 3(b) are available in the supplementary data as animation 1 and animation 2, respectively (see stacks.iop.org/NJP/16/023025/mmedia). Figures 2(c) and 3(c) show the phase space reconstruction for an embedding dimension of 2 and a time delay of 1, also known as a return plot, for 100 000 round trips of the full time series. The finite width and structure of the reconstructed attractor is one indication of chaotic, as opposed to random, motion. In figures 2(d) and 3(d) is shown correlation dimension versus embedding dimension for each variation case. Both cases saturate above an embedding dimension of 15 to a correlation dimension of 1.26 ± 0.03 and 1.66 ± 0.07, respectively. Error estimates are 95% confidence intervals given by two times the standard deviation for values of D 2 for embedding dimensions above saturation. This low dimensional chaos closely matches the low ring gain experimental observations by Wang et al [36] where 2.0% variation yields a correlation dimension of D 2 = 1.27 ± 0.12. However the numerically generated 5.1% peak variation does not reproduce The cause of D 2 collapse at embedding dimensions 6, 16 and 26 for the 5.1% modulation case has not been rigorously determined but is robust against reasonable perturbations in embedding parameters. The periodicity of the effect suggests the cause is related to sensitivities in the correlation sum to temporal correlations and finite time series. The embedding procedure is also sensitive to time series periodicity, which is present in these low dimensional examples [60,61]. Low dimensional chaos often presents as widened Fourier peaks rather than pure broadband spectra. The oscillation of D 2 for low embedding dimension is a common phenomenon as the embedding procedure is not an accurate reconstruction of phase space unless the embedding dimension is at least twice the box counting dimension of the system's attractor [59]. We find numerically that amplitude of peak modulation and the dimensionality of the chaos are principally dependent on the magnitude of the saturating quintic nonlinearity, Q. The presence of both a linear gain and nonlinear loss term is necessary for a stable correlation dimension to be determined. Chaotic modulations of the train envelope are the most complex examples of a more general modulation behavior. Parameter space explorations yielded examples of bright soliton trains with no, periodic, multi-periodic or quasi-periodic modulations. We note these types of deepening modulations were experimentally observed as the first generations of soliton fractals [30]. Dashed curve is provided as a guide to the eye; points represent actual data. Reproduced from Wang et al [36]. An animation of (b) is available in the supplementary data as 'animation 2'.

Symmetric and asymmetric interacting solitary waves
When more than one solitary wave propagate with differing group velocities, enabling dynamics such as collisions, we say these waves interact. Two distinct cases where the spatial features of solitary wave interactions are symmetric or asymmetric under rotation are discussed below.

Symmetric interaction
Symmetric interaction solutions are highly complex, but ordered, gain driven interactions between a number of intensity peaks varying from 2 to more than 20. These solutions evolve in intricate and complicated patterns but maintain symmetry in space under a rotation of π rads. The solution intensity exhibits a constrained modulation about a stable mean, but is energetically unstable. The energy of the system grows approximately linearly in time and is closely correlated, with a correlation coefficient of r > 0.95, to the time-averaged number of peaks present in the system. The sample correlation coefficient is a measure of the linear correlation between two variables and is defined as  where σ x is the standard deviation of x; P i and E i are the number of peaks and system energy at the ith round trip. This relationship suggests every intensity peak present in the system has similar energy. Peaks undergoing symmetric interactions also demonstrate persistence in time under collisions and have linear or constant phases, both characteristics of bright solitons. Further, individual intensity peaks may also be fit to a sech 2 profile when they are spatially isolated from other peaks circulating the ring. A typical example is illustrated in figure 4(a) by a spatiotemporal plot of intensity across 800 round trips, each vertical slice shows the waveform on the ring at a specific round trip. There exists a stark symmetry in dynamics with respect to a rotation by π rads. An animation of the dynamics shown in figure 4(b) is available in the supplementary data as 'animation 3', see stacks.iop.org/NJP/16/023025/mmedia. Figures 4(b) and (c) show the scaled norm and energy, respectively, for the same time frame. Over these 800 round trips we note the norm varies about a stable mean by ±1% while the system energy increases by 8%.  Symmetric interactions are observed to evolve in systems with linear gains between 10 −5 and 10 −3 and cubic losses of the same, or ±1, orders of magnitude. This long time stable evolution requires a near-balanced system where linear gain is the dominant force and peak growth is meaningfully restricted by the presence of nonlinear losses. No solutions had initial peak growth above 400% prior to the initial splitting event.
The dynamics within this regime demonstrate a characteristic splitting process, diagrammed in figures 5(a)-(d). The initial bright soliton modulates and grows until the domination of linear gain over nonlinear loss in low-intensity regions yields a non-zero intensity floor. Energy enters the system until these low lying excitations reach intensities where attractive nonlinearity and dispersion may shape the excitation into a stable solitary wave close in form to the well-known hyperbolic secant. The new peak then begins to interact with its neighbors. This same procedure results in the generation of a second, then third and so on, intensity peak. Thus, in contrast to more typical nonlinear partial differential equations which give rise to fixed soliton dynamics for all times, the GLNLS here displays a particular soliton dynamics on long but not infinitely long time-scales. This gives rise to the possibility of a new form of integrability which is relevant on long but not infinite times, and may require the development of new mathematical formalisms, in particular a multiscale approach in time. The timing of the initial splitting event varies from 100 µs to 1 ms where t = 0 is defined as the moment gain and losses are turned on. The effect of quintic loss/gain is superficial to this solution category until orders above 10 −2 when it begins to dominate the dynamics. Quintic losses (gains) result in slower (faster) rates of initial splitting, but do not have any meaningful impact on the rate of energy gain. This splitting process is stabilized (weakened) by the addition of an attractive (repulsive) quintic nonlinearity term of the same order as the cubic present in the system. Higher orders of quintic nonlinearity destroy the stability, driving the dynamics into the intermittent regimes described later in section 7. This solution type demonstrates a high sensitivity to initial conditions, which is discussed in section 8. A single round trip of a symmetric interaction solution closely resembles the multi-peaked solitons previously reported by Wu et al [35].

Asymmetric interaction
Asymmetric interaction solutions are loss driven solutions which behave similarly to the symmetrical case discussed in section 4.1 but do not maintain a spatial symmetry with respect to rotations around the ring. The number of interacting peaks was observed to vary from five to twenty depending on the parameters of the simulation. The total number of peaks is conserved, in an average sense, after spatial symmetry about the feedback ring center breaks and is closely correlated, r > 0.98, to the system's energy. Here r is the sample correlation coefficient defined by (6). An example is shown in figure 6(a) by a spatiotemporal plot of intensity over 10 000 round trips. Symmetry about the feedback ring center, x = 2.75 mm, can be seen breaking near round trip 4300. An animation of this symmetry breaking is available as animation 4 in the supplementary data. A scaled intensity plot of a single symmetric (asymmetric) round trip is shown in figure 6(b) ( figure 6(c)). The interacting peaks are seen to be node-less, and evolve about a non-zero, |u(x, t)| 2 > 0, intensity floor. The stability of the asymmetric interaction solution type is demonstrated in figures 6(d) and (e) showing scaled norm and energy, respectively, over the same 10 000 round trips. Normalization varies about a stable mean by ±0.001% while energy modulates by ±3%; this stands in contrast to the energetic instability inherent to symmetric interactions, section 4.1.
Asymmetric interactions are observed in systems with linear gains, L, between 10 −4 and 1 and cubic loss, C, of the same order of magnitude. Quintic gains and losses, Q, are stable up to this same order of magnitude. The number of peaks and peak height increased (decreased) with the presence of attractive (repulsive) quintic nonlinearity of the same order as the cubic. Higher orders of quintic nonlinearity push the solution into intermittency, a temporally unstable class of solutions discussed in section 7. The solution intensity floor varies with parameter choice, including nonlinearity, but trends towards the constant intensity which satisfies the energy balance of the GLNLS. The balance is given explicitly by the expression where L, C, Q and u(x, t) are the same terms as in (1)  This regime demonstrates a characteristic splitting process, diagrammed in figures 5(e)-(h). An initial bright soliton initial condition grows and flattens into a plateau under the influence of a strong linear gain and saturating nonlinear losses. Once the non-zero plateau expands to fill the feedback ring the central peak undergoes a splitting procedure similar to that observed for symmetric interactions, diagrammed in figures 5(a) and (b). The domination of linear gain over nonlinear losses in low amplitude regions produces small peaks. These smaller excitations grow until the system's attractive nonlinearity and dispersion shape them into solitary wave intensity peaks. Unlike the process for symmetric interactions, section 4.1, this splitting process occurs within the first 10 µs of evolution, where t = 0 is defined as the moment gains and losses are turned on, and saturates within the first 1 ms yielding an energetically stable excitation. The amplitude of intensity peaks relative to the plateau intensity varies from 1 to 10%, but the peak heights measured from the plateau mean are of the same order as those observed in symmetric interactions.
This solution type demonstrates a high sensitivity to initial conditions, which is discussed in section 8. This sensitivity and the highly complex nature of the evolution are hallmarks of chaotic dynamics. However, attempts to arrive at a converged correlation dimension using the methods discussed in section 3 were inconclusive. Such sensitivity is typically characterized by a positive Lyapunov exponent [59]. While a careful determination of the largest Lyapunov exponent requires a rigorous reconstruction of phase space we may estimate the exponent numerically by evolving nearby trajectories in time. Direct measurement suggests a Lyapunov exponent between λ = 2 × 10 4 and 1 × 10 5 s −1 . This rate of trajectory separation is of the same order as that observed experimentally (λ = 1.9 ± 0.2 × 10 5 s −1 [36]) for the 5.1% modulating soliton train discussed previously in section 3.

Dynamical pattern formation
Four distinct robust dynamical patterns which demonstrate lifetimes of at least 1 ms or 7000 round trips were located during GLNLS parameter space exploration. Solutions of this group differ from previously discussed solution behaviors in that they exhibit a periodic recurrence of their characteristic dynamic. Self organization of this kind is common in open nonlinear systems [63]. These examples are discussed to demonstrate the breadth of pattern formation supported by the GLNLS under fixed choice of N and D. The regions of parameter space supporting dynamical pattern formation violates the assumptions underlying the derivation of the GLNLS in the context of magnetic spin waves, as discussed in section 2, owing to the high order of quintic nonlinearity and losses which drive evolution. However, the GLNLS is a useful model in a variety of systems including laser cavities, as discussed in section 1, and these dynamics may appear in such contexts.

Central peak recombination
Central peak recombinations exhibit a complex five peak solitary wave recombination pattern with a periodicity of 180-250 round trips, depending on parameter choices. This behavior is driven by a strongly attractive quintic loss, Q = −1 and a linear gain of L = 10 −2±1 with cubic loss of C = −10 −3±1 . The presence of quintic nonlinearity has a severely negative impact on the behavior lifetime. The median wave height of central peak recombination solutions satisfies the energy balance equation, (7). For the example shown in figure 7 we predict an average intensity of |u(x, t)| 2 = 0.0995, corresponding to the parameters have L = 0.01, C = −0.001 and Q = −1, which closely matches the observed numerical average intensity, |u(x, t)| 2 = 0.0934 ± 6 × 10 −3 . The error estimate is defined as in section 4.2. An example of central peak recombination is shown in figure 7. Panel (a) shows a spatiotemporal plot of scaled intensity over 3500 round trips (see animation 5 in the supplementary data) and (b) shows the scaled energy over the same round trips. The periodicity of the recombination is evident in the spatiotemporal plot, which contains 16 periods. The dynamics are most readily described starting when the central peak collapses. Immediately following the collapse, the peaks on either side of the ring center propagate toward the middle of the ring and recombine into a new central peak matching the original peak's amplitude. At the same time the outlying peaks split into two. The innermost of these new peaks grows until one finds three central peaks of equal amplitude. At this point the central peak undergoes collapse and the process repeats. A single bright solitary wave propagates unperturbed along the edge of the ring.

Complex co-propagation
The complex co-propagation solution was so named as it resembles the steady state copropagation solution (see section 6.2 below) and is likewise energetically stable. It differs primarily in that the waveform undergoes complex, but periodic, modulation. The dynamics also occur on a non-zero density floor satisfying the GLNLS energy balance equation, (7). The

Spatial shifting
Spatial shifting solutions are simulations which exhibit energetically stable evolution with a well-defined and periodic shifting of the spatial location of the dynamical behaviors. In all observed cases the underlying energetically stable dynamics are evenly distributed bright solitary waves co-propagating on an intensity floor which satisfies the GLNLS energy balance given by equation (7). For the example shown in figures 8(a) and (b) we have L = 0.028 99, C = −0.062 19 and Q = 0.000 648 corresponding to |u(x, t)| 2 = 0.0438 which closely matches the numerically observed average intensity of |u(x, t)| 2 = 0.0432 ± 2 × 10 −4 . The solitary waves spontaneously split at a constant periodicity and reform into an identical set of co-propagating peaks with a spatial shift defined as L 2N s where L is the feedback ring length and N s is the number of peaks present in the simulation. All peak properties as well as the splitting dynamics remain consistent through multiple periods. A strong attractive quintic nonlinearity is required to support this dynamical behavior, as seen previously with central peak recombinations and complex co-propagation in sections 5.2 and 5.1. Spatial shifting is seen in simulations with quintic nonlinearities of S −0.8 and moderate linear gains of L = 10 −2±1 . Cubic losses near C = 10 −1 support this behavior, while quintic losses were found to be unimportant until above values of Q = ±10 −1 where they dominated the dynamics.
An example of temporal shifting is illustrated in figures 8(a) and (b). Panel (a) shows two bright solitary waves co-propagating while undergoing a spatial shift of 5.5 4 mm every 22 000 round trips (see animation 7 in the supplementary data). The shifting event occurs over 1500 round trips. Panel (b) shows the solution's scaled energy over these same round trips; the energetic stability of the co-propagation regimes is demonstrated. The energy profile of each shifting event was found to be identical.

Breathers
Solitary wave breathers on a ring are characterized by a single solitary wave which undergoes a periodic disappearance of the peak and reappearance at the other side of the ring. The frequency of breathing increases with system energy. The solution is not energetically stable and breathing frequency increases until the system reaches a new dynamical behavior. Numerically observed lifetimes were never less than 20 000 round trips, or 3 ms. The wave breathing is driven by a strong quintic loss, Q = 10 −1 , with comparatively weak linear gain, L = 10 −5±1 and cubic loss, C = 10 −4±1 , terms. The solution type is sensitive to the presence of quintic nonlinearity with any magnitude above 10 −2 , whether attractive or repulsive, pushing the dynamics into the intermittent regime, discussed later in section 7. Linear gain dominates during low intensity periods of the breathing behavior, resulting in a non-zero intensity floor which ultimately drives the collapse of stable breathing.

Steady state solutions
Simulations which evolved into energetically stable static wave forms were named steady state solutions. Two distinct steady state solutions were isolated from the parameter space exploration: multi-peaked solitary waves and co-propagating solitons.

Multi-peaked solitary waves
Multi-peaked solitary waves were characterized by energetically stable, to machine precision, nodeless complex waveforms that evolve without exhibiting any time dependence in their intensity. The shape of the wave and the number of principle peaks varies from two to eight in studied cases, depending on parameter choice. Symmetric and asymmetric waveforms were observed. Multi-peaked solitary waves were observed for any linear gain, L, below 1 and cubic losses, C, of ±3 orders of magnitude. The impact of quintic losses and gains principally affected the median wave height according to the GLNLS energy balance equation, (7). For the multipeaked solitary wave shown in figure 9 we have L = 0.1, C = −1 and Q = −1 corresponding to an estimated average intensity of |u(x, t)| 2 = 0.0916 which closely matches the observed value, |u(x, t)| 2 = 0.0916 ± 7 × 10 −7 . The error was previously defined in section 4.2. As with previous examples, high values of Q relative to L lead to the term dominating dynamics and the solution leaving the steady state solution class. Positive, or saturating, values of quintic nonlinearity lead to reductions in secondary peak heights while attractive values leads to the presence of additional principal peaks via further shaping of secondary peaks. The overall shape of the multi-peaked solitary wave, including the number of principal and secondary peaks, is dependent on the choice of parameters. A typical example is shown in figure 9 of a symmetric multi-peaked solitary wave with two principle and two secondary peaks (see animation 9 in the supplementary data). Figure 9(a) shows a spatiotemporal plot of scaled intensity over 12 000 round trips with each vertical slice showing the intensity across a single round trip. Panel (b) is the scaled intensity plot of the final round trip and panel (c) shows the static solution energy over the same evolution period. Not all multi-peaked solitary waves travel at the group velocity as the example in figure 9. This solution type is the most commonly observed long time behavior in studied simulations and was one of the behaviors present in a majority of the intermittent cases, discussed further in section 7.

Co-propagating solitary waves
Co-propagating solitary waves are the second steady state isolated during parameter space exploration. Co-propagating solitary waves are time independent solutions where N s identical bright solitary waves propagate alongside one another without interacting. Similar N s soliton solutions of the simple cubic NLS are well studied and the number of solitons is found to be proportional to the power of the initial condition relative to the value of N /D [2,23]. Periodic boundary conditions require solutions have an even number of nodes. Within studied solutions N s was observed to vary from two to eight. Figure 10 shows the same physical quantities as plotted in figure 9, with panel (c) again demonstrating the energetic stability of the solution type. The peak shape, shown in panel (b), is not consistent with either bright or dark solitons. The example plotted in panel (a) exhibits a modulation in peak heights with a variance of 10 −3 % about the mean. While the variation is not visible in panel (a) it can be observed in animation 10 of the evolution in the supplementary data (stacks.iop.org/NJP/16/023025/mmedia). Stable co-propagation was observed only in an isolated region of parameter solutions with L = 10 −4±1 and C = 10 −2 . Quintic gains, Q, of orders higher than the cubic present or quintic nonlinearities, S, with magnitude higher than 0.01 (the lowest order studied) drove the solution out of the steady state and in general pushed solutions into the intermittent class, discussed in section 7. Lower orders of quintic gain did not have any meaningful effect on stability, the number of peaks or peak height.

Intermittent solutions
Intermittent solutions demonstrate numerous distinct dynamical behaviors as the waveform evolves in time. The lifetime of these behaviors ranges from hundreds of round trips to hundreds of thousands. This corresponds to up to 1 ms before the dynamics transitions from one behavior to another. These solutions are robust to at least 10% variation of initial conditions in the sense that they do not degrade to noise or experience blow-up. Such perturbations do have significant effects on the relative lifetime of each dynamical behavior and even the types of behaviors a simulation exhibits. Quantitative matching of the intermittent dynamics to experiment will offer a challenge due to their highly transient nature; however, qualitative behaviors should be observable experimentally. In general, intermittent solutions spend a majority of their time in aperiodic evolution between distinct dynamical behaviors. Intermittent solutions can exhibit all of the behaviors previously described as temporally stable for a finite numbers of round trips. Intermittency is the typical dynamic exhibited when stable solutions are perturbed and is therefore not observed only in isolated regions of parameter space. Stable solutions are robust to variations in initial conditions, as previously stated in section 2. Intermittency is observed when perturbations exceed 10%, however it bears mention that the necessary value is ultimately highly dependent on both solution type and the parameter being perturbed. Hundreds of intermittent simulations were identified during the study.
Two illustrative examples are shown in figure 11. The same physical quantities are shown as in figure 8. Panel (a) shows a typical simulation with three distinct multi-peaked solitary wave regimes separated by two aperiodic regimes exhibiting splitting, modulation and copropagation behaviors (see animation 11 in the supplementary data). The energy is shown in plot (b) and was relatively constant during each of the multi-solitary wave regimes. The aperiodic regimes demonstrate significantly lower energy than the finite lifetime multi-solitary wave excitations. Panel (c) shows a simulation which exhibits periods of complex four solitary wave co-propagation interspersed with periods of aperiodic dynamics (see animation 12 in the supplementary data). The lengths of successive periods of dynamical behavior are highly variable and sensitive to both parameter choice and initial condition. This sensitivity makes numerical convergence difficult to demonstrate, as discussed below in section 8.

Numerical convergence and quantitative versus qualitative robustness
Simulations used well established and understood algorithms, fifth order adaptive Cash-Karp Runge-Kutta in time and pseudospectral methods in space [54,55]. Detailed evaluation of psuedospectral methods for similar nonlinear equations may be found in [54]. Simulations were run with a spatial grid of 256 and a single step truncation error of 10 −12 . The maximum number of time steps performed in a single simulation was 10 8 . Initial conditions were generated using imaginary time propagation and a single step truncation error of 10 −18 . Stability of initial conditions was confirmed via real time propagation. To machine precision all initial states exhibit zero change in energy when propagated in real time for 10 9 steps. Initial conditions for different spatial resolutions have fixed differences in energy owning to discretization. This discretization error decreases exponentially for increasing spatial resolution: |E 256 − E 512 | = 10 −7 , |E 512 − E 1024 | = 10 −8 and |E 1024 − E 2046 | = 10 −9 . Results for each solution class were compared across grid sizes of 512 and 1024 and a single step truncation error of 10 −18 to verify numerical convergence for algorithms used for both space and time propagation. We present convergence data for two distinct groupings of solution classes, those which exhibit extreme sensitivity to initial conditions and those which do not. The former demonstrate a qualitative robustness, while the latter are quantitatively robust.
Convergence can be demonstrated by relative difference. Given two sets, {x n } and {y n }, of data consisting of N directly comparable observations then the relative difference at the ith entry is defined as This quantity offers a simple, unitless measure of the relative difference between two quantities.

Quantitative robustness
Solution classes which did not demonstrate a marked sensitivity to initial conditions were numerically converged in a traditional manner. A distinct measurable, in this case energy, is quantitatively compared across successive time steps under different spatial and temporal resolutions. Convergence data is graphically displayed in figure 12 for the initial condition used during simulations, figure 12(a), as well as four categorical behaviors, panels (b)-(e). In each case the solid blue line compares spatial resolutions of 256 and 512 grid points, while the dashed red line compares the spatial resolutions of 512 and 1024 grid points. Figure 12(a) shows the fixed discretization error discussed previously in section 8 while figures 12(b)-(e) demonstrate the spatial convergence of each dynamical behavior over the entire evolution period is as good or better than that of the initial conditions. The greatest observed single time step relative spatial resolution error was 10 −3 %. The greatest observed single time step relative temporal resolution difference was 10 −8 %. The solution types listed here were quantitatively converged: • complex co-propagation; • spatial shifting; • breathers; • multipeaked solitary waves; and • co-propagating solitary waves.

Qualitative robustness
A subset of observed dynamical behaviors, from both the temporally stable and intermittent categories discussed in sections 3-7, demonstrate an extreme sensitivity to initial conditions. These solutions were robust to variations in initial conditions and parameters of at least 10% in the sense that such perturbations did not yield a shift in their categorization. However, changes in initial energy or in loss parameters of the order 10 −9 and lower resulted in distinct dynamics within that categorical behavior and shifts in the starting and ending times. We note shifts in spatial resolution introduce variations of this order to the relaxed initial condition. Therefore solutions exhibiting this sensitivity may not be converged numerically in the traditional sense. An illustrative example is given by the asymmetrical interaction behavior, discussed in section 4.2, after spatial symmetry about the feedback ring center has broken. Figures 13(a)-(c) depict scaled intensity across the same 2500 round trips for three different choices of spatial resolution: 256, 512 and 1024 spatial grid points, respectively. The behaviors across the three spatial resolutions are qualitatively similar with each exhibiting an asymmetric N s solitary wave interaction which is characteristic of the solution type. However, the detailed dynamical behaviors of each case are quantitatively different. Further, as shown graphically in figure 13(d), the three solutions have energies which vary by less than 0.32%. The behaviors listed below all demonstrated sensitivity similar to that discussed here with a relative energy difference no greater than 0.1%: • symmetric interaction; • asymmetric interaction; • chaotic modulation; and • central peak recombination.
This sensitivity is a typical property of evolution toward and around a strange attractor. Major attempts were made to quantify the dimensionality of the attractor as discussed in section 3 by the authors. No stable correlation dimensions were located for any of the complex dynamical behaviors which demonstrated sensitivity to initial conditions, excluding chaotic envelope modulation. A Lyapunov exponent was estimated for the asymmetric interaction case,

Unstable
Any simulation which evolved into a trivial result, degraded to noise, blew up or decayed to zero is considered to be unstable. Approximately 30% of studied simulations are unstable. This is not surprising as the explored parameter space includes cases when either gain (or loss) dominate the dynamics by orders of magnitude. A subset of simulations are observed to degrade into noise due to spatial resolution issues. Pseudospectral methods rely on discrete fast Fourier transforms (dFFT) which provide excellent convergence for well-behaved curves. However, if the spatial features of u(x, t), the complex spin wave amplitude, approach the length of the numerical lattice spacing singularities may appear and the dFFT algorithms will no longer converge locally. These errors will continue to grow and propagate over time in an L , Q > 0 evolution. Further exploration of these cases with finer spatial resolutions is prohibited by computational resource constraints. In contrast, all results previously presented in sections 3-6 are converged in both space and time, as demonstrated in sections 8.1 and 8.2.

Conclusions
We report the numerical identification of nine distinct long lifetime complex dynamical behaviors as part of six broad solution classes of the GLNLS, (1). Behaviors were located during an extensive numerical exploration of six dimensional parameter space. A minimum of eight decades were examined for each gain term while five decades of higher order nonlinearities were considered at fixed dispersion and cubic nonlinearity. The GLNLS served as a driven damped model of long lifetime spin wave dynamics in magnetic thin film active feedback rings and analogous driven damped nonlinear physical systems. Agreement of GLNLS low dimensional chaotic modulating bright soliton trains with experimental measurements [36] was discussed in detail. We predicted additional GLNLS dynamical behaviors including two distinct steady state solutions, four unique examples of stable dynamical pattern formation and the intricate spatially symmetric/asymmetric interactions of solitary wave peaks. Finally we reported the existence of intermittent regimes within GLNLS parameter space, a phenomena typical of chaotic dynamical systems. Two unique examples of intermittency were presented which demonstrated finite periods of two distinct dynamical behaviors. The variety of presented GLNLS solution types matches the scope of dynamical behaviors observed experimentally in YIG film spin wave systems, as well as predicting new behaviors that can be tested in present experiments. The GLNLS thus presents a simple yet viable and fundamental model for driven, damped nonlinear waves propagating in dispersive mediums.
We neglected the periodic effect of amplification within the feedback ring, so the gain and loss terms presented in this work represented averaged quantities. Highly variable solution types such as the symmetric and asymmetric interactions potentially violate the GLNLS operating regime, with gain and loss driven dynamics occurring on the scale of a single round trip. Future studies of this limit, adiabatically driven soliton trains and transient dynamics are warranted. A fine grained exploration of parameter space may also be justified to identify distinct domains of stability for each observed behavior. In the future a rigorous study of GLNLS phase space would be useful to determine the cause of intermittency and potentially locate chaotic attractors of higher dimension.