{'fraction_non_alphanumeric': 0.039342425613513456, 'fraction_numerical': 0.012474953138129404, 'mean_word_length': 4.4482920530578705, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coupled complex Ginzburg-Landau systems with saturable nonlinearity and asymmetric cross-phase modulation

We formulate and study dynamics from a complex Ginzburg-Landau system with saturable nonlinearity, including asymmetric cross-phase modulation (XPM) parameters. Such equations can model phenomena described by complex Ginzburg-Landau systems under the added assumption of saturable media. When the saturation parameter is set to zero, we recover a general complex cubic Ginzburg-Landau system with XPM. We first derive conditions for the existence of bounded dynamics, approximating the absorbing set for solutions. We use this to then determine conditions for amplitude death of a single wavefunction. We also construct exact plane wave solutions, and determine conditions for their modulational instability. In a degenerate limit where dispersion and nonlinearity balance, we reduce our system to a saturable nonlinear Schr\"odinger system with XPM parameters, and we demonstrate the existence and behavior of spatially heterogeneous stationary solutions in this limit. Using numerical simulations we verify the aforementioned analytical results, while also demonstrating other interesting emergent features of the dynamics, such as spatiotemporal chaos in the presence of modulational instability. In other regimes, coherent patterns including uniform states or banded structures arise, corresponding to certain stable stationary states. For sufficiently large yet equal XPM parameters, we observe a segregation of wavefunctions into different regions of the spatial domain, while when XPM parameters are large and take different values, one wavefunction may decay to zero in finite time over the spatial domain (in agreement with the amplitude death predicted analytically). While saturation will often regularize the dynamics, such transient dynamics can still be observed - and in some cases even prolonged - as the saturability of the media is increased, as the saturation may act to slow the timescale.

Spatio-temporal chaos has been observed in the complex Ginzburg-Landau equation [11,12]. Such dynamics can then be related to soft or hard turbulence [13]. Transition to spatio-temporal turbulence from traveling waves due to instabilities was discussed in [14]. Spatial patterning is also possible, with solutions such as spiral waves being found [15]. These dynamics, while irregular, appear in a bounded region within C, and are distinct from parameter regimes giving blow-up solutions. Bounds on the attractor for complex Ginzburg-Landau equations have been considered in [16,17]. There have also been studies on the control of such turbulence [18][19][20][21].
Coupled or vector complex cubic GL equations have attracted interest. For example, the effect of slow real modes in reaction-diffusion systems close to a supercritical Hopf bifurcation was studied in [44] via a system of two coupled GL equations, where it was pointed out that a single GL equation cannot reproduce the dynamics from the reactiondiffusion system, suggesting the necessity of considering such coupled systems in some contexts; see also [45]. Such systems arise in various applications in nonlinear optics [46]. Linearly coupled systems of GL equations were considered in [47,48], where solitary pulses were studied. In the case where the squared-modulus of the waves feature into each equation (in contrast, say, to the case of linear coupling [48]), the cross-phase modulations parameters may vary. This is modeled by assuming different coefficients for the self-and cross-interaction terms, as was done in [49]. Phase instabilities in vector GL equations have been studied [50]. A coupled GL system formulation of the weak electrolyte model for electroconvection in planarly aligned nematic liquid crystals has been considered [51]. A coupled complex GL system for the amplitudes of waves along a gas-less combustion front was derived in [52]. There have also been studies on discrete coupled GL systems [53,54]. These replace the Laplacian operator with a discrete analogue, and can be viewed as an oscillator system on a network.
A variety of solutions to systems of coupled cubic GL equations have been studied in the literature. Stability or instability of solutions of two coupled complex GL equations was considered in [49,55]. Instability and transition from chaotic traveling waves to stable states such as soliton lattices was studied in [56,57]. The effect of spatial frequency forcing on standing-wave solutions of coupled complex GL equations was considered in [58]. Regarding spatial patterning and spatio-temporal chaos, [59] show that the generalized synchronization in spatially extended chaotic GL systems is determined by the additional dissipation introduced into the response system, and that a coupling parameter increase is equivalent to a simultaneous growth of the dissipation and the amplitude of the external signal. Spatio-temporal chaos in coupled complex GL equations was studied in [60]. Creation and annihilation processes of different kinds of vector defects in coupled GL equations are considered in [61], as are transitions between different regimes of spatio-temporal dynamics. Phase singularities in such systems are studied in [62].
Although far less studied than their cubic counterparts, coupled cubic-quintic GL equations have also been studied in the literature. Modulational instability of linearly coupled equations was considered in [63]. The collision of pulses or of solitary waves has also attracted attention [64][65][66]. The spatio-temporal structure of pulsating solitons was approximated using the variational approximation [67]. Modulational instability of plane waves in two coupled cubicquintic GL equations coupled with a XPM term (in the cubic nonlinearity alone) was considered in [68] and also [69], while the dynamics of stationary pulse solutions for such equations were studied in [70]. Synchronization of solutions in such systems was studied in [71].
Polynomial nonlinearities are not the only possibility for self-or cross-interaction terms, however cubic or cubicquintic nonlinearity is what is standard in the literature. In the case of nonlinear Schrödinger equations, a variety of other self-interaction terms are often studied. Saturable nonlinearities consist of a rational or other bounded function of the complex modulus of the macroscopic wave function [72,73]. These find use in nonlinear optics, with application to nematic liquid crystals [74,75], optically induced nonlinear photonic lattices [76], and have been used to study lefthanded properties of metamaterials with saturable nonlinearity [77]. Despite this interest in saturable nonlinearties, there do not appear to be vector or even scalar extensions of complex GL equations which employ such nonlinearity. This extension shall be the focus of the present paper.
In the present paper we shall study a coupled complex GL system with saturable nonlinearity and XPM parameters in the nonlinear terms, taking the form where , a, b ∈ R, µ > 0, n = 1, 2, . . . , and α 1 , α 2 are XPM parameters. To the best of our knowledge, saturable nonlinearity has not been considered in complex GL dynamics, although for the simpler NLS dynamics, the saturable n = 1 case was considered in [72,73,76,77], while the n = 2 case was considered in [74,75], all without arbitrary XPM terms. When µ = 0, we recover the cubic GL system and hence both µ and n determine respectively the strength and type of saturable term. Regarding XPM in cubic GL systems, note that the literature consists of papers which take α 1 = α 2 , and often both are set to unity. The case of different XPM parameters has apparently not been considered under cubic GL dynamics (coupled complex cubic [46] or cubic-quintic [68][69][70] GL equations coupled with the two XPM parameters always equal have been considered but with very specific solutions sought and generic dynamics not explored, while coupled complex GL equations with two distinct XPM parameters do not appear to have been written down let alone explored). As we shall see, interesting dynamics will arise from (1)-(2) when α 1 = α 2 , even when µ = 0 (corresponding to the cubic GL system (3)-(4)). This suggests that even in the µ = 0 limit, this system may give new dynamics for some combinations of XPM parameters. The remainder of this paper is organized as follows. In Section 2, we study the system (1)-(2) analytically, approximating the absorbing set for the solution trajectories and demonstrating generic modulational instability of the solutions. One interesting phenomena observed is extinction of one of the wavefunctions. This depends strongly on the XPM parameters -and we arrive at local conditions for this to occur for either u or v. We give analogous results for the scalar saturable GL equation, as well. In Section 3, we verify the aforementioned analytical results numerically, while also demonstrating other features of the system dynamics. As suggested analytically, there exists an absorbing set for sufficiently small , otherwise the solutions blow-up for large time. Spatio-temporal chaos is common, although there are some parameter regimes where coherent patterns -such as uniform states or banded structures -appear to persist. A partitioning of the support of wave functions, with the wave functions tending to separate into distinct regions over time, is observed in many cases. In the extreme case, one wavefunction will decay to zero with the other persisting, as the local analytical results suggested. In addition to these global behaviors, we observe a variety of highly local structures, and while spatio-temporal chaos appears quite common, we note for some parameter regimes the appearance of transient defects and what appear to be rogue waves at intermediate time intervals. The primary role of the saturation appears to be regularization of the cubic GL dynamics and a modification of the timescale of the dynamics, depending on the parameter regime, and the kinds of dynamics we observe can all be observed in cubic GL systems (save for blow-up or high-intensity solutions which are particular to the saturable model) provided XPM parameters are taken to differ in value from self-interaction parameters. We summarize the interesting findings and discuss our results in Section 4.

II. ASYMPTOTIC DYNAMICS OF SATURABLE COMPLEX GL SYSTEMS
In the present section we shall study the asymptotic dynamics of (1)-(2) under reasonable assumptions. We first construct an approximation to the absorbing set for solution trajectories in (u, v) ∈ C 2 , and later numerical simulations validate the results we obtain despite the fact that these results are asymptotic in nature. Simulations suggest that the dynamics of (1)-(2) are often modulationally unstable, and we are able to demonstrate that there always exist unstable wavenumbers for generic parameter values which induce modulational instability. One interesting feature observed in numerical simulations is the dissipation or extinction of one wavefunction under certain restrictions on the XPM parameters, and we are able to determine these conditions analytically. Finally, all of our results for the system (1)-(2) can be deduced for the scalar saturable GL equation, and we list those results for completeness.
A. Approximating the absorbing set for the system (1)- (2) In order to construct an approximation to the absorbing set for solutions in C 2 , we shall make the simplifying assumption that there is one dominant spatial mode for sake of analytical tractability. Assuming that arbitrarily small modes are excited on the plane, we shall then take the wavenumbers corresponding to dominant modes to be arbitrarily small. While this is not generically true, we note that asymptotic results obtained under this assumption show agreement with results from numerical simulations (which shall be considered later in Section 4).
Consider a solution of the form Such a solution will give the dynamics corresponding to a single wavenumber. The transformation (5) puts (1)-(2) into the form Consider U (t) = ρ u (t) exp(iθ u (t)) and V (t) = ρ v (t) exp(iθ v (t)), which puts (6)-(7) into the form Observe that we may write and solving for θ u we obtain Similarly, This then leaves us with two coupled ODEs for ρ u and ρ v , (8) and (10). Define the intensity We shall observe that it behaves as a Lyapunov function. Note that I(0, 0) = 0, while I |u| 2 , |v| 2 > 0 for all |u| 2 , |v| 2 = (0, 0), and it is radially unbounded. Note that As I is equivalent to the Euclidean norm on R 4 , we have four-volume contraction when dI dt < 0 and expansion when dI dt > 0.
To better understand the structure of (16), write ρ u = r cos(φ) and ρ v = r sin(φ) so that the ρ u -ρ v plane is put into polar coordinates. This puts (16) into the form where Here F has domain s ≥ 0 and χ ∈ [0, 1]. We shall only be interested in solutions for which I(t) remains bounded -that is, solutions with bounded modulus for all time. Therefore, we require that F < 0 as s → ∞. From (18), we have Since n ≥ 1, we have that F → −∞ as s → ∞ for Now, we desire solutions which have bounded modulus for arbitrary wavenumbers, so we set |k u | = |k v | = 0 as these will include the least stable wavenumbers. With this, we have There are two cases to consider. First, if α 1 + α 2 ≥ 2, then as we must maintain stability independent of the angle taken toward s → ∞ in the plane. On the other hand, when α 1 + α 2 < 2, we have as for this case the least stable trajectories as s → ∞ correspond to χ = 1 2 . Therefore, we have obtained the condition One may show that > (α 1 , α 2 , µ, n) will give F → ∞ as s → ∞ for at least some χ ∈ [0, 1]. Let us next note that F (0, χ) = 0. Meanwhile, Then, when > 0 there exist arbitrarily small wavenumbers for which F > 0 in a neighborhood of the origin. This, in turn, implies instability near the origin (instability of initial solutions with arbitrarily small modulus). On the other hand, if < 0, then the origin is stable. As we shall be interested in bounded yet nontrivial dynamics, we consider the case where > 0.
The above analysis suggests that in order to obtain bounded yet non-uniform dynamics, we should consider 0 < < (α 1 , α 2 , µ, n). Now that we have this condition, we shall attempt to approximate the attractor. Note that our analysis presupposes that a fixed wavenumber will be excited, while in general multiple wavenumbers may be excited simultaneously. We shall therefore consider the case where the wavenumbers are zero, which is the least stable case to consider. This will then provide an approximation to the absorbing set, which can be viewed as being a bound for the attractor [78].
One may verify that for µ > 0, n ≥ 1, 0 < < (α 1 , α 2 , µ, n), the function ∂F ∂s can have either zero or one turning points; either it decreases monotonically toward a fixed limit given by lim s→∞ F (s, χ)/s or it deceases, then increases monotonically toward a fixed limit given by lim s→∞ F (s, χ)/s. By our choice of bounds on , the limit will always be negative. Therefore, ∂F ∂s > 0 at s = 0, lim s→∞ ∂F ∂s < 0 , and since there is at most one turning point, for each χ ∈ [0, 1, there will be a unique s = s * * (χ) at which ∂F ∂s = 0. This, in turn, suggests that along any fixed χ ∈ [0, 1], there will exist a unique s = s * (χ) at which F is maximal (and, positive). Hence, along any fixed χ ∈ [0, 1], we have that F increases from zero to F (s * (χ), χ), and then decreases without bound (F → −∞ as s → ∞). From this and continuity of F , there is a unique curve s = ζ(χ) > 0 such that F (ζ(χ), χ) = 0.
To find the curve ζ(χ), we set F (ζ, χ) = 0 and, using ζ > 0, we have which gives the relation Let us again assume arbitrarily small wavenumbers are excited. Then, we have Returning to coordinates (r, φ), we have that r 2 = ζ sin 2 (φ) . Returning to original coordinates, we have the curve where φ = tan −1 (|V |/|U |) provided |U | = 0 and φ = π 2 if |U | = 0. From the above analysis, the curve is closed and bounded away from the origin and away from infinity, given µ > 0, n ≥ 1, 0 < < (α 1 , α 2 , µ, n). Putting this all together, let use define the sets and We have that dI dt > 0 on B i and dI dt < 0 on B o . As I is simply a scaling of the Euclidean norm, we therefore have that trajectories grow in norm on B i while they contract in norm on B o , suggesting that as t → ∞ we have that trajectories (U, V ) converge to the attractor A.
While the curve given in terms of ζ may seem a bit abstract, note that we may bound ζ in the |U |-|V | plane. Define constants Then, the surface A ∈ C 2 is enclosed between a 3-sphere of radius min {ζ 1 , ζ 2 } and a 3-sphere of radius max {ζ 1 , ζ 2 }. Which of ζ 1 or ζ 2 is greatest will depend on α 1 + α 2 ≷ 2. This suggests that, when arbitrary small wavenumbers are excited, the large-time dynamics will be found within the region min {ζ 1 , ζ 2 } < |U | 2 + |V | 2 < max {ζ 1 , ζ 2 }. This should be more numerically feasible to detect than the explicit form of the attractor A.
Note that, in the case where there is a minimal non-zero wavenumber excited (bounded away from zero), the general idea is the same, with the form of ζ including |k u | and/or |k v |, as needed.

B. Generic modulational instability of plane wave solutions
Having constructed an approximate absorbing set for the dynamics of (1)-(2), we note that simulations (to be discussed in Section 4) generically show that solutions do not tend to globally uniform steady dynamics. This motivates us to study the modulational stability of plane wave solutions. As our simulations suggest, we are able show that for generic parameter combinations there will exist a family of unstable wavenumbers leading to modulational instability. As we shall show in Section 4, this gives the appearance of spatio-temporal chaos in the dynamics for many parameter regimes.
We begin by looking for the plane wave solutions where A, B, ω u , ω v ∈ R and k u , k v ∈ R n . Substituting these solutions into (1)-(2) and separating the real and imaginary parts we find from which we find The conditions for existence of plane wave solutions (35) are therefore and Under these conditions the plane wave solutions will exist in the absorbing set A.
In order to determine the stability or instability of the plane wave solutions (35), we follow the procedure described in [69]. We first introduce perturbations to the plane wave solutions of the form where δ 1. Substituting the perturbed solutions (43) into (1)- (2), we find that the O(δ) perturbations are given by where u * is the complex conjugate, and the real-valued constants F u , F v , G u , and G v are given by We take the ansatzû which puts (44) into the form of the matrix system where the elements of M are In order to obtain a non-trivial solution for the amplitudes we require that the determinant of the matrix M to be zero, giving a characteristic polynomial for Ω: Ω 4 + Λ 3 Ω 3 + Λ 2 Ω 2 + Λ 1 Ω + Λ 0 = 0, where the Λ j 's are functions of the entries in M . We find that Λ 3 = −i (M 11 + M 22 + M 33 + M 44 ), while the other components will not be of use to us. From the form of (46), the plane waves (35) will be unstable if Im(Ω) > 0 for at least one of the four complex roots.
Note that if Im(Λ 3 ) < 0, then at least one of the roots has positive imaginary part. To show this, assume we have a fourth order polynomial with roots Ω j with j = 1, . . . , 4. Then, we can show (from Vieta's formula) that Therefore, if Im(Λ 3 ) < 0, then 4 j=1 Im(Ω j ) > 0, which means that at least one of the roots Ω j must have positive imaginary part. Therefore, Im(Λ 3 ) < 0 is a sufficient condition for instability.
From the definition of M we have that so the instability condition is given by Im(Λ 3 ) = − 4|K| 2 + 2 (F u + G v ) + 4 < 0, which is satisfied for any K such that As we can choose K arbitrarily, we conclude that the system is always modulationally unstable on an unbounded domain.

C. Local stability conditions for single wavefunction extinction
In numerical simulations we occasionally observe one of the wavefunctions u or v tends to zero as t → ∞. We shall now derive conditions under which one of the (u, v) = (0, v) or (u, v) = (u, 0) states are locally stable. These conditions, in turn, will give parameter restrictions for which one of the two extinction states is stable. The XPM parameters will play a strong role in these conditions, with the results suggesting that situations where α 1 = α 2 = 1, which are commonly taken in the literature when studying cubic GL systems, do not result in this wavefunction extinction, while asymmetric parameter conditions can allow for dissipation of one wavefunction as t → ∞.
For large time we were able to show that solutions tend to the absorbing set A, so let us consider a solution on the absorbing set, again assuming that µ > 0, n ≥ 1, 0 < < (α 1 , α 2 , µ, n). Then, parameterizing the moduli of the solutions to (8) and (10) as ρ u = r cos(φ), ρ v = r sin(φ), we obtain Now, on the absorbing set, the r equation is satisfied for r = ζ sin 2 (φ) , while the φ equation becomes We have that and We arrive at the stability criteria for extinction of one of the wavefunctions. Given µ > 0, n ≥ 1, 0 < < (α 1 , α 2 , µ, n), the state (u, v) = (u, 0) (corresponding to φ = 0) is locally stable on the absorbing set A provided that while the state (u, v) = (0, v) (corresponding to φ = π 2 ) is locally stable on the absorbing set A provided that D. Asymptotic dynamics from a single scalar saturable complex GL equation While the system (1)-(2) is novel, note that a complex scalar GL equation with saturable nonlinearity has not been considered in the literature, with only NLS equations with saturable nonlinearity attracting attention thus far. Therefore, the scalar equation is of interest in its own right. However, note that obtaining results for such a complex scalar equation is not quite as simple as setting the XPM terms to zero (i.e. α 1 = α 2 = 0), as the quantity |u| 2 + |v| 2 appears in the denominator of the system (1)-(2), so one should start with the scalar equation and directly obtain the needed results. We sketch the obtained results here, for sake of completeness. To begin with, the scalar complex saturable GL equation reads A plane wave solution of the form u(x, t) = Ae i(k·x−ωt) with exists for (58) provided 0 < − |k| 2 < µ −1/n . As we have done in the previous section, this plane wave solution may be shown to modulationally unstable against arbitrary perturbations with wavenumber K satisfying |K| 2 > − − |k| 2 1 − µ − |k| 2 n . To approximate the absorbing set, we again write u = e ik·x U (t) = ρ(t)e iθ(t) e ik·x , where θ = b − (a + b)|k| 2 t − b log(|ρ|) and the dynamics of ρ are given by Define the intensity for the scalar equation (58) as I = |u| 2 . Note I(0) = 0, I(|u| 2 ) > 0 for non-zero |u|, and I is unbounded as |u| grows. Further, Now, F (ρ) < 0 if < |k| 2 + ρ 2 ((1+µρ 2n ) 1/n , and for asymptotically large ρ this gives the condition < µ −1/n , which is the first of the two possible parameter restrictions we obtained in the two-species case. Further, we have that F (ρ) > 0 for ρ ∈ (0, ρ 0 ) and F (ρ) < 0 for ρ ∈ (ρ 0 , ∞), where Therefore, the norm I will expand for ρ < ρ 0 and contract for ρ > ρ 0 , making ρ = ρ 0 the absorbing set. These results are in line with the two-species case, only here we obtain a circle of radius ρ 0 in C as the absorbing set rather than a more complicated submanifold of C 2 . Note that the solution is only bounded for <¯ = 0.8964. We take random initial data in (0, 1) for both plots. Panel (a) demonstrates a solution which remains bounded for all time, while in panel (b) we find that the solution grows without bound for large enough time, owing to the fact that the value corresponding to this solution is larger than¯ .

III. NUMERICAL SIMULATIONS AND EMERGENT DYNAMICS
We numerically solve equations (1)-(2) using finite-differences in Matlab. We write four separate PDE for the real and imaginary parts of each wavefunction, and use centered differences to discretize the Laplacians, and solve the resulting large system of ODEs using the Runge-Kutta solver 'ODE45'. We use absolute and relative tolerances of 10 −9 and constrain the time step to confirm convergence of solutions. We pick suitably large domains of the form [0, L] and [0, L] 2 with periodic boundary conditions to approximate R and R 2 respectively. We take L = 10 3 in all simulations. Unless stated otherwise, we always use random initial conditions in the interval (0, 0.1). In the 1-D setting we used m = 10 3 grid points, and in the 2-D setting we used m = 400 × 400 grid points.

A. Blow-up and boundedness
We numerically confirm the boundedness of solutions consistent with the results approximating the absorbing set in Section II A. For all values of α 1 + α 2 ≥ 2 we find that solutions remain bounded for large time for all parameter ranges we explored if <¯ . This bound is sharp, as for >¯ , random small initial data grow arbitrarily large and hence appear to become unbounded. This is a very rapid and spatially uniform process, in that within 1000 units of time t the intensity is larger than 10 300 . In the other case, where α 1 + α 2 < 2, we again confirm boundedness for <¯ , and unbounded solutions for ¯ , but there are more complicated behaviours for ¯ . We show an example of this in Figure 1, where bounded initial data remain so for sufficiently small. When this bound is exceeded, however, small pockets of large intensity appear as time progresses, and these appear to become larger over time, growing without bound. We confirm this by increasing the magnitude of the initial data, our simulations showing an increase in the magnitude of the spikes later in time.

B. Separation of intensities
In this section, we demonstrate an interesting effect for large XPM (that is, when α i > 1) whereby the intensity of each wavefunction, defined as I u = |u| 2 and I = |v| 2 respectively, become spatially separated. That is, for large times and for α i 1.4, we always observe that each respective intensity is nonzero only when the other is approximately zero. In Figure 2 we plot these corresponding intensities for three levels of α 1 = α 2 in each plot to demonstrate this effect. For α 1 = α 2 = 1, plots instead appear like that in Figure 1(a).
Within these regions of nonzero intensity (or in the case of equal intensities where each wavefunction is mixed) the generic behave of the intensity is a kind of spatiotemporal chaos. In some cases, this induces a meandering boundary of the region where one is defined. We note that some regions that have nonzero intensity for one wavefunction can disappear due to this chaotic motion changing the boundary of this region; see for instance, Figure 3. In each case, we observe a separation of intensity of the two wave functions. As the XPM parametersα1 and α2 increase in value, we find that the solutions tend to segregate within bands. Within a band, the dynamics appear to give spatio-temporal chaos, which verifies the fact that modulational instability of wavefunctions appears generic to the system.
If α 1 = α 2 , and one of the conditions (56) or (57) holds, then the equation containing the larger XPM parameter will be driven to zero. We give examples of this in Figure 3. In the first example (corresponding to α 1 = 1.5, α 2 = 1.4), the XPM parameters are not large enough to fully separate the wavefunctions, and hence I u persists in sporadic pockets over time. In the second example (corresponding to α 1 = 5, α 2 = 4), I u is driven to zero as I v tends to a spatiotemporal-chaos state throughout the domain. We observe that some 'tendrils' of u extending in time, corresponding to small regions of nonzero I u that persist in time, but these are eventually extinguished by fluctuations of the larger wavefunctions.

C. Comparison to the cubic complex Ginzburg-Laundau system
We now vary µ to compare the behavior of the system with saturable kinetics to that of the cubic complex GL system. We note that the boundedness result for < 1/µ 1/n (assuming α 1 + α 2 > 2), implies that for a given , solutions will be bounded for µ < 1/ n , and we exploit this to consider large values of µ in the saturable model.
We use increasing values of µ and plot intensities I u and I v in Figure 4. Generally, the differences between the saturable dynamics and the µ = 0 cubic case depends on how close, in a relative sense, the parameters are to the blowup boundary of¯ . We note that the timescales are taken larger with µ, in order to observe the behavior. Similarly, the maximum intensity of each plot increases as µ increases. We note qualitatively different transient behaviors occur for changes in the saturation parameter µ.
We note that for large values of the XPM, for some parameter regimes, we observe one of the wavefunctions becoming extinct (as in the case with unequal α i ). In Figure 4, the top left panel exhibits such behavior, as does a panel in Figure 5. We note that this is not unique to saturable kinetics; the top left panel of Figure 4 corresponds to the purely cubic kinetics. However, since the complex GL system with cubic nonlinearity is not usually studied with general XPM parameters, we are not aware of this being pointed out previously in the literature. We observe that larger values of µ tend to smooth and coarsen the dynamics, or in some cases regularize them as in Figure 4 where for µ = 10, both intensities tend to steady but simply-patterned states, whereas for smaller values of µ the intensities behave as more generic spatio-temporal fluctuations. A complete classification of when the dynamics coalesce to form steady intensity profiles is not tractable, but we do note that this behavior is extremely rare if a and b have the same sign (with spatio-temporal chaos appearing far more often in simulations), but common if they have opposing signs as in Figure 5. Note also that, for a and b of opposing signs, we observe that while there is still a coarsening of the dynamics for large saturability, there may be loss or extinction of one of the wavefunctions, as seen in Figure 5(e). Therefore, the coexistence of smooth, temporally steady patterns depends strongly on not only the saturability of the media, but on the parameters a and b, as well.

D. Localized structures
In addition to spatio-temporal chaos or regular steady state patterns such as bands or even uniform states, there are other emergent dynamics possible from (1)-(2) in restricted parameter regimes, at least at intermediate timescales.
In Figure 6, we observe transient defects in the intensity plots that eventually deteriorate into spatio-temporal chaos. However, as we increase µ, we see these transient defects persist for longer periods of time, again emphasizing the kind-of regularizing nature of being closer to the blow-up parameter boundary. The triangular shaped transient defects have been observed previously in the literature on vector cubic GL systems [60]. For µ = 0, we recover the cubic complex GL system, and with the choice of different XPM parameters, we have the extinction of the u wavefunction for large enough time. As µ increases, the solutions are in a sense regularized, and the two wavefunctions are seen to both persist for µ = 1. The dynamics of both µ = 0 and µ = 1 appear to exhibit spatio-temporal chaos. Interestingly, when µ becomes quite large (here we set µ = 10), the dynamics segregate into discrete bands for large time, with these structures appearing to become stable. The amplitude of these solutions is large (relative to the other cases) indicating that the solutions grow in amplitude and become self-reinforcing against perturbations when saturability of the medium is strong, for our choice of a and b. For these parameters (in particular, for a < 0), we find that small values of µ are somewhat stabilizing, resulting in banding patters akin to those found in Figure 4. However, these bands are fairly irregular in appearance, with some bands eventually merging at large time. This is in contrast to the discrete nature of the bands observed in (e)-(f) of Figure 4. Increasing µ from µ = 0 to µ = 50, we find that the overall structure of the solutions is the same, with the only qualitative change being an increase in the amplitude of the solutions. However, when we increase µ further, there is a bifurcation into a new regime, with initial attempts at banding instead dying out (e) or merging into a single region (f). This provides an example of a parameter regime where increasing the saturability in the model actually results in a loss of regularity in the discrete banded structures, with one of the wavefunctions dissipating completely for large time. Note that the transient defects are more prevalent when µ is large. So, not only do the defects persist for longer, they also are a more dominant proportion of the dynamics. Form small µ the transient defects will give way to spatio-temporal chaos, while when µ is larger they may tend to uniform or existence states. For larger µ we also observe a resurgence of new transient defects (see Figure 6(d)), with new defects forming after the initial defects have dissipated. This is in contrast to the small-µ limit (which includes the cubic complex GL system), where the initial transient defects will give way to spatio-temporal chaos. This again suggests one role of µ is to regularize the dynamics of the saturable GL system.
Finally, we note the existence of what appear to be rogue waves -isolated spatial maxima of amplitude greater than three times the surrounding wave height -that occur within pockets of spatio-temporal chaos. We note examples of these in Figure 7, where these highly localized waves occur as bursts in similar spatial regions to where they occurred previously. While they share some features of breathers, they exhibit far less regularity and appear to always die out after sufficient time. It is interesting to see such dynamics even when saturation is very large, as an increase in µ often results in a stabilization of the dynamics. Optical rogue waves were previously reported in higher-order (fourth-order, in particular) complex scalar GL equations [10], and so the emergence of rogue waves in certain parameter regimes for our model may not be surprising.

E. Dynamics and pattern formation in two spatial dimensions
Finally we demonstrate patterns observed in two-dimensional spatial domains with specific initial data. Generic behavior for small random initial conditions corresponds to random fluctuations as we have seen in the 1-D case, and we also remark that for small µ we can recover many kinds of patterns found in the cubic complex GL system in the literature. Here we give an example of spontaneous patterning arising from an initial seed of u(0, x, y) = Within each band there appear to be highly localized waves with peak amplitude much larger than that of surrounding waves. The localization is in space and in time, with the waves pulsing several times before dissipating. As the amplitude of these waves is up to and greater than three times the mean amplitude of surrounding waves, these localized structures may be candidates for transient rogue waves. (1 + i) cos(x) cos(y) and v(0, x, y) = (1 + i) sin(x) sin(y). In Figure 8 we demonstrate the evolution of this initial data for α 1 = α 2 = 2, and in Figure 9 for α 1 = α 2 = 5. Note that the color bars have been scaled to emphasize regions of high-intensity.
As we might anticipate from the one-dimensional simulations, for α 1 = α 2 = 1, patterns form and mingle between the two wavefunctions, without any clear separation between them (although we do notice that they pattern differently, likely corresponding to interactions between them). For larger XPM parameters, we see separation between the two wavefunction intensities as they form intricate patterns on the boundaries between high intensity areas. Over long periods of time, these patterns break down into spatio-temporal fluctuations as expected, but they seem to persist in a dynamic way (e.g. continually changing and shifting) for long periods of time depending on parameters chosen.

IV. DISCUSSION
We studied the dynamics emergent from a complex Ginzburg-Landau system with saturable nonlinearity including cross-phase modulation (XPM) parameters. In the case where the saturation parameter is set to zero, we recover a general complex cubic Ginzburg-Landau system with XPM. First studying the system analytically, we obtained conditions for the existence of bounded dynamics. Within the parameter regime where bounded dynamics exist, we were able to approximate the absorbing set for solutions analytically, showing that it is bounded between two 3spheres in C 2 . We were also able to demonstrate generic modulational instability of plane wave solutions for arbitrary values of the model parameters. This suggested that transient rather than steady dynamics should prevail when bounded solutions exist, as was later verified with numerical simulations. For XPM parameters much larger than self-interaction parameters (the latter normalized to unity), we found from numerical simulations that segregation of wavefunctions into different regions of the spatial domain was a common occurrence. Within each respective region, the underlying dynamics for XPM parameters set to unity were found, with behaviors ranging from steady uniform dynamics or patterning to spatio-temporal chaos depending on the values of other model parameters. Furthermore, when XPM parameters were large and yet of distinct values, one wavefunction was seen in simulations to decay to zero in finite time over the spatial domain. In order to better understand these numerical results, we obtained analytic conditions for amplitude death of one of the wavefunctions in the presence of unequal XPM parameters, and we found amplitude death of u for large enough α 1 and amplitude death of v for large enough α 2 .
As the generic modulational instability suggests, spatio-temporal chaos is common in the bounded dynamics regime; however, there are certain parameter regimes where coherent patterns including uniform states or banded structures persist even for large time. Again, depending on the XPM parameters, these steady state solutions may separate into different regions of the spatial domain, or may only exist for one of the wavefunctions (due to amplitude death of the other wavefunction). While numerical simulations helped in the verification of large-time asymptotic results obtained in Section 3, they also were essential for the detection of other interesting transient features of the dynamics not apparent from the asymptotic results. For instance, we found a variety of transient dynamics which persist for intermediate times, including transient defects (such as those of triangular form shown in Figure 6) and spatial patterns in two-dimensional spatial domains (as shown in Figures 8-9). Additionally, in some simulations we found what appear to be rogue waves, with single large amplitude waves surrounded by regions of far smaller amplitude background waves.
When the XPM parameters exceed the self-interaction parameters by a sufficient amount, there is a separation of wavefunctions. When all four parameters are the same (or roughly of similar value) then both wavefunctions frequently exhibit spatio-temporal chaos over the spatial domain. However, when the XPM are taken to be larger than the self-interaction parameters (we always scale the latter to unity), then we observe a separation of wavefunctions into distinct regions of the spatial domain. The size and extent of the regions changes with time, and the support of one wavefunction may occupy a larger region of the spatial domain at one time, only for the support of the other wave function to dominate at a later time. Within the support of each wavefunction, spatio-temporal chaos is often found (see, for instance, Figure 2 for an example), although for some specific parameter regimes localized structures or stable patterns may emerge. One example is the development of banded patterns, such as those seen in Figure  4. For large enough XPM parameters, these bands will occur for each wavefunction over mutually exclusive subsets of the spatial domain. This has not been frequently observed in the literature for cubic complex GL systems, as self-interaction and cross-interaction terms are often take to be the same. However, in the limit where µ = 0 (where saturation is negligible and we recover the cubic nonlinear kinetics), we observe the same behaviors for large enough XPM parameters, suggesting that these dynamics occur in appropriate complex GL systems generically, independent of the precise form of nonlinearity so long as the nonlinearity includes self-and cross-interactions, and provided that cross-interactions (XPM parameters in our model) are sufficiently dominant.
We also observe extinction, or amplitude death, of one wavefunction when there is sufficient asymmetry in the coefficients in the nonlinear terms. Indeed, while large enough α 1 = α 2 > 1 promotes a segregation of the individual wavefunctions, we find amplitude death of one wavefunction when α 1 and α 2 are large and also take different values. In particular, if α 1 is dominant, then we observe amplitude death of u, while if α 2 is dominant, then we observe amplitude death of v. One such case is observed in Figure 3. We note that amplitude death due to delayed coupling in GL systems has very recently been observed [79], where it was shown that amplitude death can occur in a pair of one-dimensional CGL systems coupled by diffusive connections. The analytical results of [79] reveal that amplitude death never occurs in a pair of identical complex cubic GL systems coupled by a "no-delay connection", but can occur in the case of the delay connection considered. In our work, using XPM parameters different from unity and obeying certain parameter restrictions, we observe such amplitude death in perhaps a simpler or more natural way than introducing delayed couplings. The results of [79] stating that coupled GL equations with equivalent kinetics do not have amplitude death hold true for our saturable model as well, as this corresponds to the case of α 1 = α 2 = 1 in our model, which never satisfies the parameter restrictions (56)- (57). Amplitude death or partial amplitude death is more commonly seen in literature on network or lattice equations such as oscillator systems [80][81][82][83][84], and is seemingly novel in reaction-diffusion PDE systems.
In many cases, we observe that saturability (as measured by increasing the parameter µ) can regularize or stabilize dynamics. This is apparent in Figure 4 where increasing µ takes the system from amplitude death (µ = 0, the cubic GL system limit) to coexistence of wavefunctions in a segregated manner with each giving spatio-temporal chaos (µ = 1) to finally steady banded patterns (µ = 10). We have carried out many more simulations than shown here, observing that µ can also change the timescale of the dynamics. While we have not commented much on the role of n in the saturable terms, we note that an increase in n will result in a stabilization of the dynamics and a slowing of the timescale. Recall that a condition < Cµ −1/n derived in Section 2 determines boundedness of dynamics (with the most interesting bounded dynamics occurring near the boundary of this condition in parameter space), so n can be seen as a scaling of this condition, with increasing n regularizing the dynamics. In our numerical simulations of Section 3 we focused on n = 1 and n = 2 as these are the two values commonly taken in the literature. However, the tend of increasing n regularizing the dynamics is seen as we take n ≥ 3, as well.
These results suggest that one might modify the saturability of the media in order to develop a novel control mechanism for the spatio-temporal chaos emergent from complex GL equations and systems, which has not apparently been considered in the earlier mentioned control papers. Additionally, one could attempt to modify the XPM parameters in order to better control the placement of wavefunctions or force amplitude death if that is desirable. Used in conjunction, saturability and XPM parameters could be used to position and reinforce the intensity of specific patterns. We are unaware of this approach being used in the literature.
While often regularizing, we observe that saturable terms can prolong the existence of transient defects, such as those found in [60]; see Figure 6. Furthermore, there are regimes where saturability permits blow-up of solutions yet no such blow-up occurs in the cubic nonlinearity case of µ = 0. This suggests that the dynamics are highly sensitive to the model parameters, even in the presence of strong saturation in the media. It was interesting to observe rogue waves even when saturation was very large, as seen in Figure 7. Unlike for cubic nonlinearity in the NLS, calculating the rogue wave envelope analytically is made more challenging due to the saturable nonlinearity, and so we only have numerical indications of the emergence of rogue waves. We are not aware of the emergence of rogue waves in saturable models, although there have been experimental suggestions of rogue waves in saturable media. Optical rogue waves were previously reported in complex cubic scalar GL equations [10] with higher-order derivatives, and the emergence of rogue waves in certain parameter regimes for our model for weak saturability is therefore not surprising. Although not much has been done in this direction involving saturable nonlinearity, we point to a very recent study [85] which found candidate rogue waves numerically in a nonlinear Schrödinger equation with saturable nonlinearity. Like what we have shown for our solutions in Figure 7, [85] found isolated large-amplitude waves surrounded by background waves of much smaller mean wave height. A more systematic study on how saturability of a medium might influence the stability of highly localized structures, such as rogue waves or breather-like solutions, would be an interesting topic for future work.