Coupled Nonlinear Schr\"odinger (CNLS) Equations for two interacting electrostatic wavepackets in a non-Maxwellian fluid plasma model

The nonlinear dynamics of two co-propagating electrostatic wavepackets, characterized by different wavenumbers and amplitudes, in a 1D non-magnetized plasma fluid model is considered, from first principles. The original plasma model, consisting of \kappa-distributed electrons evolving against a cold ion background, is reduced, by means of a multiple-scale perturbation method to a pair of asymmetric coupled nonlinear Schr\"odinger (CNLS) equations for the dynamics of the wavepacket envelopes. Exact analytical expressions are derived for the dispersion, self-modulation, and cross-modulation coefficients involved in the CNLS equations, as functions of the wavenumbers and the spectral index \kappa characterizing the electron profile. An analytical investigation of the modulational instability (MI) properties of this pair of wavepackets reveals that MI occurs in most parts of the parameter space. The instability windows and the corresponding growth rate are calculated in a number of case studies. Two-wave interaction favors MI by extending its range of occurrence and by enhancing its growth rate. Growth rate patterns obtained for different \kappa suggest that deviation from Maxwellian equilibrium, for low \kappa values, leads to enhanced MI of the interacting wave pair. To the best of our knowledge, the dynamics of two co-propagating wavepackets in a plasma described by a fluid model with \kappa-distributed electrons is investigated thoroughly with respect to their MI properties as a function of \kappa for the first time, in the framework of an asymmetric CNLS system. Although we have focused on electrostatic wavepacket propagation in non-Maxwellian plasma, the results are generic and may be used as basis to model energy localization in nonlinear optics, in hydrodynamics or in dispersive media with Kerr-type nonlinearities where MI is relevant.


Introduction
Research on plasma dynamics has grown significantly during the last decades, because of the prospects for a wide range of applications related to space and fusion, among others.Plasmas in a state of collisional equilibrium exhibit a Maxwellian distribution.Space plasmas, however, are very often in stationary states out of equilibrium that can be described by kappa distributions, which are characterized by a Maxwellianlike "core" and a high-energy tail that follows a power law.Such distributions have been observed in the solar wind and planetary magnetospheres [1,2,3,4,5].The intrinsic dispersion and nonlinearity in plasmas can be balanced to lead to collective dynamics and the emergence of localized excitation in the form of solitons, breathers, rogue waves, or shocks.Plasma fluid models mimicking Navier-Stokes equation in hydrodynamics (see e.g. in [6,7]), can be reduced through standard perturbative methods to the celebrated nonlinear Schrödinger (NLS) equation known to govern the envelope dynamics of a modulated wavepacket in nonlinear and dispersive media.In the context of plasma fluids, the NLS equations has been derived for the first times in the early '70s [8,9], and subsequently for a large variety of electrostatic plasma models.Recently, its localized solutions (envelope solitons, breathers) have in fact been associated with freak or rogue waves in plasmas [10,11,12,13].Formally analogous NLS equations focusing on electromagnetic waves in plasmas modeled by fluid-Maxwell equations have been also derived [14,15,16].In recent years, the kappa distribution of electrons has been used more and more in NLS-reduced plasma fluid models [17,18,19,20,21,22,23].
Different versions of NLS equations have been also employed in other contexts, e.g., in the investigation of the dynamics of blow-up solutions in trapped quantum gases in Gross-Pitaevskii equation [24] which in a sense generalizes the three-dimensional NLS equation.Alternatively, plasma fluid model may be reduced to the Kortewegde Vries (KdV) equation, which has been recently investigated in a generalized form (i.e., with high power nonlinearities) with respect to solitary, compacton, and periodic types of solutions [25].Furthermore, the NLS equation with higher order dispersion and parabolic nonlinearity has been found to admit dark, bright, kink, as well as triangular solitary wave solutions [26].Also, several types of localized excitations have been constructed for extended Kadomtsev-Petviashvili equations of the B type, such as dromion-and lump-lattice structures as well as various folded solitary waves [27,28].
More specifically, pairs of CNLS equations have been employed as the workhorse model for systems that belong to diverse areas of science such in nonlinear optics, hydrodynamics, and dispersive media, for the investigation of their dynamics and their modulational instability properties.Below we refer to some of these work in which formally equivalent forms of CNLS equations to those presented in Section 2.3 or 2.4 were employed.
Systems of CNLS equations have been obtained in Ref. [70] for the description of interacting nonlinear wave envelopes and formation of rogue waves in deep water.A condition for rogue wave formation was introduced which relates the angle of interaction with the group velocities of these waves.The CNLS system exhibit modulational instability (MI), which develops much faster (i.e., it has higher growth rates) than in each of the corresponding decoupled CNLS equations, and constitutes an important aspect of rogue wave dynamics.Further, CNLS equations were employed as a model describing the formation of rogue waves of the Peregrine breather type in water, to confirm corresponding experimental evidence.Indeed, in Ref. [69], it was demonstrated that nonlinear wave focusing resulting in strong localization on the water surface which originates from MI evolving from two counter-propagating water waves.The observed experimental results are in excellent agreement with the corresponding results obtained from the hydrodynamic CNLS equation model.In another relevant work [53], a weakly nonlinear model consisting of a pair of CNLS equations was derived for two water wave systems known as crossing sea states in the context of ocean waves which propagate in two different directions.It was shown that when the first wave is modulationally unstable, then a second wave propagating in different direction results in an increase of the MI growth rates and enlargement of the instability region through their interaction.
Furthermore, a pair of CNLS equations which describe the interaction between the high-frequency Langmuir and low-frequency ion-acoustic waves in a plasma, and breather and rogue wave solutions have been investigated via the Darboux transformation [52].The co-propagation of two circularly polarized strong laser pulses in a magnetized plasma in a weakly relativistic regime and their MI properties were investigated in Ref. [38] within the framework of two CNLS equations.Further works on dispersive media include the investigation of a CNLS equations model resulted from Maxwell's equations with nonlinear permittivity and permeability in a lefthanded metamaterial [57], for which several types of vector solitons were identified, the investigation of the MI profile of the coupled plane-wave solutions in left-handed metamaterials using CNLS equations in Ref. [10], and the investigation of the coupling between backward-and forward-propagating waves, with the same group velocity, in a composite right-and left-handed nonlinear transmission line [71].
Perhaps the largest amount of work on CNLS equations has been performed in optics, in which long ago those equations were shown to govern light propagation in isotropic Kerr materials [54].In that work, the existence of a novel class of vector solitary waves was revealed, which bifurcate from circularly polarized solitons.Also, the generation of vector dark-soliton trains by induced MI in a highly birefringent fiber was experimentally observed in Ref. [55], in excellent agreement with theoretical predictions from numerical simulations of the CNLS equations.More recent works include the emergence of temporal patterns such as bright solitons and rogue waves in optical fibers investigated theoretically via a special version of CNLS equations, i.e., the Manakov system, which are initiated by the MI process [56].The Manakov system has been also recently investigated with respect to the existence of novel, nondegenerate fundamental solitons corresponding to different wavenumbers [58].More optical solitons have been obtained for the Manakov system of CNLS equations in Ref. [59], while analytical one-and two-soliton solutions in birefringent fibers have been obtained through the Hirota bilinear method in Ref. [74].
In the present work, a pair of CNLS equations is derived from a plasma model comprising of a cold inertial ion fluid evolving against a thermalized electron background that follows a kappa distribution, using a standard multiple scale approach (the Newell method).The coefficients of the CNLS equations have been calculated as functions of the wavenumbers of the two interacting waves and the coefficients of a Mc Laurin expansion of the kappa distribution (which depend on the spectral index κ).For arbitrary wavenumbers, the coefficients of the CNLS equations do not exhibit any known symmetry and hence the system is rendered non-ingtegrable.A compatibility condition is derived for that system through a detailed modulational instability analysis, that provides a polynomial relation (fourth order in the perturbation frequency) between the perturbation frequency and wavenumber [60].The MI properties of the CNLS equations are then investigated with respect to the wavenumbers of the two carrier waves and the wavenumber of the perturbation for several values of the spectral index κ of the electron distribution.
In comparison with earlier similar works, which in the context of plasma are rather scarce, we should stress that the obtained CNLS equations for the envelops of the two interacting wavepackets are the most general as far as the values of their coefficient are concerned; those coefficients are derived exactly in terms of k 1 , k 2 , and κ without any other additional assumption (e.g., equal group velocities of the wavepackets).The choice of the fluid model was made on the basis on simplicity; it contains only one parameter, the spectral index κ (besides the externally determined k 1 and k 2 ).This constitutes advantage to other models since it isolates the effect of κ on the investigated MI properties.Moreover, a thorough MI analysis is performed as a function of κ for the first time, which illuminates its role for values that are relevant to space plasmas and have been deduced from observational data.A comparison of the results for small values of κ (in the interval [2,3]) with those using the common Maxwell-Boltzmann electron distribution reveals dramatic differences as far as the strength of the growth rate of MI and the MI window are concerned, which have not been shown before.
In the following Section 2 the plasma fluid model equations are introduced, from which the CNLS equations are derived using the multiple-scale perturbation method.In Section 3, we undertake an analytical study of the modulational instability analysis of plane wave solutions of the CNLS equations.A detailed parametric analysis of MI occurrence is carried out, based on numerical calculation of the growth rate in Section 4. Our main results are briefly discussed and summarized in Section 5.
2 Derivation of a system of CNLS equations for two interacting wavepackets

Plasma fluid model equations and electron distribution
A one-dimensional (1D), non-magnetized plasma model is considered here for electrostatic excitations, that consists of a cold inertial ion fluid of number density n i = n and velocity u i = u, evolving against a thermalized, kappa-distributed electron background of number density n e .The continuity, momentum and Poisson partial differential equation(s) governing the dynamics of that model are given by [61] respectively, with φ being the electrostatic potential.In Eqs. ( 1)-(3) above, the spatial and temporal variables x and t, are normalized to the inverse ion plasma frequency ω −1 pi and to the Debye screening length λ D , respectively, where The number density variable(s) of the ions n and the electrons n e are normalized to their respective equilibrium values n i,0 and n e,0 = z i n i,0 , the ion velocity u is normalized to i.e. essentially the (ion-acoustic) sound speed.Finally, the electrostatic potential φ is normalized to k B T e /e.Adopting a standard notation, the symbols ε 0 , k B , T e , e, m i , z i appearing in the latter expressions denote respectively the dielectric permittivity in vacuum, the Boltzmann's constant, the (absolute) electron temperature, the electron charge, the ion mass, and the degree of ionization (z i = q i /e, where q i is the ion charge).
As announced above, we shall assume the electrons to follow a kappa-distribution [1,62,63] that is given by [64,65] Note that this function is characterized by a real parameter κ (known as the spectral index).The kappa distribution adopted here for the nonthermal electrons is a straightforward replacement for the Maxwell-Boltzmann distribution when dealing with systems in stationary states out of equilibrium [66].Such stationary states are commonly found in space and astrophysical plasmas [2,3], whose nature allows for extracting the spectral index κ in each particular case.Possible values of κ in Eq. ( 7) lie in the interval (3/2, ∞), while for κ = +∞ the kappa distribution reverts to the Maxwell-Boltzmann one.Lower values of κ are associated with a long tail in the particle (velocity) distribution, in account of a strong suprathermal component (i.e., for large velocity) [63].This is a common occurrence e.g. in Space plasmas [1].
Space plasma observations have shown an empirical separation for the spectrum of the values of κ into a near-(3/2 < κ < 2.5) and a far-(2.5 < κ < ∞) equilibrium region.In particular, kappa distributions with 2 < κ < 6 have been found to fit the observations and satellite data in the solar wind [62] (and references therein).In later sections, we are using values of κ evenly distributed in the interval 2 < κ < 6 as well as a large value of κ = 100 for which practically the electron component in the considered plasma model has reached equilibrium (Maxwell-Boltzmann).
Kappa distributions have by now become one of the most widely used tools for characterizing and describing space plasmas and are known to affect the fundamental properties of plasma waves and also of electrostatic solitary waves [67].For low values of the electrostatic potential φ, i.e. for |φ| ≪ 1 (in normalized units), the kappa distribution Eq. ( 6) can be expanded in a Mc Laurin series in powers of φ as from which we shall keep terms up to order φ 3 .In that case, the first three expansion coefficients are given by

Application of the reductive perturbation method
We introduce fast and slow variables so that the derivatives are approximated as where ε is a small quantity, and x k = ε k x, t k = ε k t (for k = 0, 1, 2, 3, ...).The dependent variables n, u, and φ can be expanded around their equilibrium value as The expansions of {n, u, φ} in Eqs. ( 13)- (15) along with those for the derivatives in Eq. ( 12) and the expansion of the electron distribution in Eq. ( 7) are substituted into the model equations Eqs. ( 1)-(3), in order to obtain sets of perturbative equation of different orders in ε up to order ε 3 .The fluid equations obtained at first order (∝ ε 1 ) read where x 0 = x and t 0 = t.In these equations, the following Ansatz will be introduced 1,1 e +iθ1 + S 1,2 e +iθ2 + c.c., where S 1 = {n 1 , u 1 , φ 1 } and the phases are θ j = k j x − ω j t with k j and ω j being the wavevector and the corresponding (angular) frequency of the j−th wave (for j = 1, 2).Note that the derivatives with respect to t = t 0 and x = x 0 that appear in Eqs. ( 16)- (18) do not act on the amplitude(s) S 1,j --(but only on the harmonic carrier functions, i.e. the exponentials, above), since the latter does not depend on these variables.After some straightforward calculations we get the equations (j = 1, 2) +n 1,j = 0.
The equations for j = 1 are not coupled to those with j = 2 and thus the two waves do not interact at this order of approximation.As a result, two uncoupled 3 × 3 systems with unknowns {n 1,j , φ 1,j } for j = 1 and j = 2, respectively, are obtained as For these systems to have a non-trivial solution, it is required that their determinants should vanish, i.e., which provides the linear frequency dispersion relation (ω) and the corresponding group (v g ) and phase (v ph ) velocity for the j−th wave (j = 1 or 2) as where the dependence of the expansion coefficient c 1 on the spectral index κ is explicitly indicated.The systems of Eqs. ( 21) can be solved in terms of the variables The frequency dispersion ω 1 and the group velocity v g,1 which have been calculated Fig. 1 The frequency dispersion ω 1 (a), the corresponding group velocity v g,1 (b) as calculated from Eq. ( 24), and the phase velocity v ph,1 = ω 1 /k 1 (c) as calculated from Eq. ( 25), as a function of the wavenumber k 1 for κ = 2 (black), κ = 3 (red), κ = 4 (green), κ = 5 (blue), κ = 6 (orange), κ = 100 (brown).
from Eqs. ( 23) and ( 24) are shown in Fig. 1 as a function of the wavenumber k 1 for six values of the spectral index κ.The frequency dispersion curves seem to have the same qualitative features for any value of κ.However, the frequency (or fixed k 1 ) increases with increasing κ for any k 1 in the (physically relevant) interval shown, while it saturates for large spectral index κ to its Maxwell-Boltzmann value.Note that the curve obtained for κ = 100, shown as brown curve in Fig. 1(a), is practically identical to a curve that would have been obtained using the Maxwell-Boltzmann distribution for the electrons.
The corresponding group velocity curves in Fig. 1(b), exhibit two distinct behaviors; for low k 1 (i.e., k 1 0.9) the group velocity v g,1 increases with increasing κ just as the frequency does in Fig. 1(a).However, for high k 1 (i.e., k 1 0.9) that tendency is reversed, and the group velocity v g,1 decreases with increasing κ.Note that for very low k 1 (i.e., for k 1 → 0), the group velocity increases from v g,1 ≃ 0.59 for κ = 2 to v g,1 = 1 for κ = 100 (practically the same as in the case of Maxwell-Boltzmann distributed electrons, as expected for large values of κ).This constitutes a relative increase of up to ≃ 70% in group velocity (value) between the extreme values of κ adopted here.We draw the conclusion that lower values of κ, i.e. a stronger deviation from the thermal (Maxwellian) picture result in lower group velocity (i.e.slower wavepackets/envelopes) for electrostatic wavepackets of this kind, at least for long wavelengths (i.e. the Debye length, more or less).The inverse trend is observed for short (subscreening-radius) wavelengths (i.e. for wavenumber exceeding unity, roughly), since lower κ values entail a slightly faster v g (value).
The phase velocity v ph,1 , on the other hand, exhibits a simpler behavior, as illustrated in Fig. 1(c).It turns out that v ph,1 increases with increasing κ for all wavenumbers, while at the same time the phase speed decreases asymptotically towards zero for large k 1 (i.e. for vanishing wavelength), regardless of the value of κ, as expected from the analytical expression (25).For negligible k 1 (i.e. in the infinite wavelength limit), the expressions for both group and phase speeds recover the sound speed lim k1→0 v g,1 = lim k1→0 v ph,1 = 1/ c 1 (κ), as can be inferred from Eqs. ( 24) and (25).
Obviously, all of the observations in the latter two paragraphs hold for the second wave also, upon a simple subscript permutation 1 → 2.
The equations obtained at the second order in ε are where F 1 , F 2 , and F 3 are known functions of Ψ j which are given in Appendix B. To solve Eqs. ( 28), we use the following ansatz 2,1 e 2iθ1 + S (1) 2,+ e i(θ1+θ2) + S (1) where S = n, u or φ; by "c.c" we have denoted the complex conjugate of the expression within the curly brackets.Upon substitution of Eqs. ( 29) into Eqs.(28), we obtain the systems of equations 2,j = µ 3,j , resulting from all terms proportional to e iθ1 and e iθ2 , respectively (to be discussed later; note that µ 1,j , µ 2,j , and µ 3,j , involving temporal and spatial derivatives of Ψ j , for j = 1, 2, are given in Appendix B), for j = 1 and 2, the latter resulting from the terms proportional to e 2iθ1 and e 2iθ2 , respectively.The second harmonic amplitudes are given by solutions in the form where the coefficients u,2,j , and C φ,2,j are given in Appendix A, and finally 2,+ = f and where f 2 , and f are functions of k 1 and k 2 and are given in Appendix B. The last two systems result from the terms proportional to e i(θ1+θ2) and e i(θ2−θ2) , respectively.Their solution reads: where again the coefficients C n,2,± , C u,2,± , and C φ,2,± are given in Appendix A. In the same (second) order (∼ ε 2 ), we also obtain an equation involving the "constant" terms (zeroth order harmonics) in the form which will be used in the next subsection.Before moving on to the equations which are third order in ε, we return to the discussion of the solution of Eqs.(30).The latter, which actually consists of two 3 × 3 uncoupled systems for the variables {n 2,1 } and {n 2,2 }, can be written in matrix form as where µ m,j (m = 1, 2, 3 and j = 1, 2) are provided in Appendix B. Note that the determinants of the two inhomogeneous systems of Eqs. ( 13) are zero, since this was the requirement for obtaining the frequency dispersion relations.Thus, in order to find nontrivial solutions for the two inhomogeneous systems of Eqs.(13), it is imposed that the following determinants should also vanish: These determinants result from the replacement of the first column of the determinant D j with the vector µ j = [µ 1,j µ 2,j µ 3,j ] T .After calculations we obtain that in order for D ′ j to vanish, the following compatibility conditions (for both j = 1 and 2) should hold.The physical implication is that the (two) wavepacket envelopes travel at their respective group velocity.This was qualitatively expected from earlier works [10,68].
Isolating the ℓ = 0 contributions, we obtain two equations for the zeroth harmonic amplitudes which can be combined with Eq. ( 42) to provide a 3 × 3 system to be satisfied by the (3) unknowns n In writing down the latter equations, we have also used the relations (27).To solve the 3 × 3 system thus obtained for the zeroth harmonics, we shall introduce the Ansatz into Eqs.( 48)-( 50) and use Eq. ( 15) to determine the coefficients C φ,2,j , C u,2,j , and C n,2,j whose expressions are given in the Appendix A. Finally, the equations obtained (in third order in ε) for the (3rd-rder) first harmonic mode ℓ = 1 read where the functions R 1,j , R 2,j , and R 3,j are given in Appendix B.
A lengthy algebraic manipulation of Eqs. ( 52) -aimed at eliminating the variables n where the expressions for R m,j (m = 1, 2, 3 and j = 1, 2) are given by Eqs.(B20)-(B22) in Appendix B. Finally, by substituting the expressions for R m,j into Eq.( 53) and using the obtained solutions for the harmonic amplitudes n 2 , and φ (0) 2 , we obtain after a tedious and lengthy calculation with equations where Note that the terms proportional to the variables n 2,j were eliminated with the help of Eqs.(30).
From Eqs. ( 54) and ( 55), it is straightforward to obtain a pair of coupled NLS equations in the familiar form where (for j = 1, 2) are the (2) dispersion coefficients (analogous to the group-velocity-dispersion/GVD terms in nonlinear optics; note that P j = 1 2 ∂ωj ∂kj , and are the self-modulation (self-nonlinearity) coefficient and the (2) cross-modulation (coupling nonlinearity) coefficients, respectively.(The tilded quantities in the latter three expressions were defined above.)In Fig. 2, the dispersion coefficient P 1 and the nonlinearity coefficient Q 11 are shown as a function of the wavenumber of the first carrier wave k 1 .Note that the corresponding quantities for the second carrier wave (i.e., P 2 , Q 22 , etc.) as a function of the wavenumber k 2 have exactly the same form.As evident in Fig. 2(a), for the model considered here, the dispersion coefficients P j (j = 1, 2) are negative for any k j , converging to zero for large values of k j .Also, for the interval of k 1 values shown in the figure, the magnitude of P 1 increases with increasing spectral parameter κ.In the Maxwell-Boltzmann case (i.e., corresponding practically to the value of κ = 100), the minimum value of P 1 is located at k 1 ≃ 0.5, while it shifts to high values of k 1 with increasing κ.
On the other hand, the nonlinearity coefficients Q jj can be either positive or negative, depending on the value of k j , as can be seen in Fig. 2(b).In the Maxwell-Boltzmann case, that corresponds to very high values of the spectral index κ of the kappa distribution of the electrons, the coefficients Q jj change sign at a particular value of the wavenumber of the carrier wave k j , i.e., at k j ≡ k r = 1.47, in agreement with earlier results [8,9,61].
Critical carrier wavenumbers for non-Maxwellian plasmas.Generally, for arbitrary κ, the symbol k r denotes the value of k 1 at which Q jj has a single root, i.e., k r is defined as that k 1 for which Q jj (k j ; κ) = 0.As it can be implied from the definition of Q jj (for either j = 1 or 2), k r is a(n undisclosed) function of κ.
The dependence of the root k r on the spectral index κ is shown in Fig. 3(a) for a large interval of κ values.As it can be observed in Fig. 3(a), the root k r diverges as κ → 3/2.The value of the root k r decreases rapidly as κ increases from the value 3/2, then reaches a minimum around κ ≃ 2.5 close to the border of far-to near-equilibrium, and then, with further increasing κ, it increases slowly to converge to k r = 1.47 from below for very large values of κ (Maxwell-Boltzmann regime).The coupling coefficients Q 12 and Q 21 , on the other hand, depend on both k 1 and k 2 as can be observed from Eqs. (62).It turns out that, on the k 1 − k 2 plane, these coefficients may have either positive or negative values.The boundary between the regions corresponding to opposite signs for Q 12 is shown for several values of κ in Fig. that might be associated with negative Q 12 for Maxwellian electrons (i.e.large κ) now switch to a positive Q 12 in the case of suprathermal electrons, i.e. for smaller κ.The analogous curves (boundaries) for Q 21 are obtained upon reflection with respect to the diagonal k 1 − k 2 , hence in the areas above the boundaries Q 21 is positive (and negative underneath).
It is certainly useful to see not only these boundaries but all the values that Q 12 and Q 21 can take on the k 1 − k 2 plane.The first of them, i.e., the coefficient Q 12 is mapped on the k 1 − k 2 plane in Fig. 4 for several values of κ. (Note that Q 12 and Q 21 exhibit reflection symmetry around k 1 = k 2 .)As it can be observed in Fig. 4, the coefficient Q 12 is mostly negative on the k 1 − k 2 plane for κ = 2 (Fig. 4(a)), while the areas of negative and positive values are more or less of the same size for the other values of κ, i.e., for κ = 3, 4, 5, 6, 100 (Figs.4(b)-(f)).Also, it is clear that Q 12 takes very large positive or negative values in certain areas of the k 1 − k 2 plane.In particular, Q 12 is very large and positive for low k 2 , while it is has very large negative values in the upper left corner of the k 1 − k 2 plane.Actually, the plots of Q 12 are shown only for k 2 > 0.1 to avoid very large and divergent values.These large values of the coupling coefficients Q 12 and Q 21 are probably responsible for the wide occurence of modulational instability in the CNLS system that will be discussed in the next section.Interestingly, only one (i.e.not both) of the Q 12 and Q 21 coefficients can become zero for a (any) particular pair of k 1 and k 2 values.In that case, for, say, Q 21 = 0, the second of the CNLS equations Eq. ( 60) decouples from the first one, Eq. ( 59), although the latter is still affected by the former [61].

Reduced form of the CNLS system
Equations ( 49) and ( 50) can be transformed to a more common form of CNLS equations with a change of the independent variables, followed by a transformation of the wavefunctions Ψ 1 and Ψ 2 .Let us introduce the new independent variables ξ and τ through ξ = x − vt and τ = t, where the subscript in x and t have been dropped, and v = (v g,1 + v g,2 )/2 is the half-sum of the group velocities v g,1 and v g,2 .By applying this change of variables to Eqs. ( 49) and ( 50) we get where δ = (v g,1 − v g,2 )/2 is the half-difference of the group velocities.In writing the Eqs.( 63) and ( 64), it is implicitly understood that the differentiations involved in the dispersive terms are with respect to the variable ξ 1 , while the higher-order differentiation in the beginning of the left-hand-side are with respect to ξ 2 .The subscripts are nonetheless omitted, for simplicity, as they will not affect the final result.Then, by applying to Eqs. ( 63) and ( 64) the transformation Ψ j = Ψj exp i δ 2 4Pj τ − δ 2Pj ξ , where j = 1, 2, we obtain where Ψj are complex functions of the new variables ξ and τ .
It should be noted that CNLS equations analogous (and formally equivalent) to either Eqs.( 59)- (60) or to Eqs. ( 65) and (66) above, have been obtained in various physical contexts earlier, describing nonlinear processes including pulse propagation in water waves [69], deep water wave propagation [70], soliton propagation in left-handed transmission lines [71], pulse propagation in optical nonlinear media [72], emergence of vector solitons in left-handed metamaterials [57], propagation of pulses with different polarization in anisotropic dispersive media [73].Analogous systems of equations have been obtained from various plasma fluid models [29,30,32,33,34,35,37,38,61], but never in the case of a (non-symmetric) pair of electrostatic wavepackets (i.e. with different wavenumbers and amplitude), which is the focus of this study.

Modulational instability analysis -analytical setting
Modulational instability (MI) analysis for two co-propagating pulses in the plasma fluid described by Eqs. ( 1) can be performed on the CNLS equations ( 63) and (64) following the procedure in Refs.[38,60].We shall first seek a plane-wave solution for the system, which turns out to be of the form where Ψ j,0 (j = 1, 2) is a constant amplitude (here assumed to be real) and ωj is a frequency, to be determined.By inserting Eq. ( 67) into Eqs.( 63) and ( 64), we obtain Equations ( 67) with the amplitude-dependent frequencies given in Eqs. ( 68) and ( 69) provide a nonlinear mode of the CNLS system wherein the propagation of waves is governed by Eqs. ( 63) and (64).The stability behavior of these modes against a small perturbation can be addressed through standard MI analysis.Let us perturb the nonlinear modes by adding a small complex quantity ε j to their amplitudes (|ε j | ≪ |Ψ j,0 |), as Ψ j = (Ψ j,0 + ε j ) e iωj τ .(70) By inserting the perturbed Ψ j into Eqs.( 63) and ( 64) we obtain Then, by setting ε j = g j + ih j , where g j and h j are real functions, and subsequently substituting into Eqs.( 71) and ( 72) and eliminating h j from the equations, we get Eventually, by setting g j = g j,0 e i(Kξ−Ωτ ) + c.c, (75) where c.c. denotes the complex conjugate, and then substituting into the earlier equations, we obtain after some algebra the compatibility condition where (j = 1, 2) For modulational instability to set in, Eq. ( 76) must have at least one pair of complex conjugate roots, so that the (positive) imaginary part of these roots can be identified as the growth rate of the modulationally unstable modes Γ.When Eq. ( 76) has two pairs of complex conjugate roots, then the largest of their imaginary parts is identified as the growth rate of the modulationally unstable modes, Γ, i.e., where Ω r denotes one (each) of the four roots of the compatibility condition.
Note that a similar analysis may be applied to incoherent pairs of CNLS equations [74] as well as to higher dimensional, non-autonomous CNLS equations [46].

Decoupled envelopes.
In order to facilitate later discussion, we consider now the special case in which the CNLS equations ( 73) and ( 74) are decoupled, i.e., the case in which either Q 12 or Q 21 (or even both) are equal to zero.In that case, the compatibility conditions Eq. ( 76) simplifies significantly, and the perturbation frequencies Ω as a function of the perturbation wavenumber K are given by for the first and the second NLS equation that results from decoupling Eqs. ( 73) and ( 74), respectively.Correspondingly, the criterion for modulational instability in the first and second decoupled equation becomes, respectively, P 1 Q 11 < 0 and P 2 Q 22 < 0. Note that these products depend on k 1 only and k 2 only, respectively.Since P j (j = 1, 2) in the considered model is always negative, the sign of the product P j Q jj depends here on Q jj .Now, if we recall that Q jj > 0 for k j < k r (κ) while Q jj < 0 for k j > k r (κ), we arrive at the following result: Within the considered model, the first NLS equation (for the 1st wavepacket) is modulationally stable for k 1 < k r (κ), while it is modulationally unstable otherwise.An identical conclusion is drawn and can be formulated for the second wavepacket (if decoupled from the first), upon formally shifting 1 → 2, i.e. it will be stable for k 2 < k r (κ) and unstable otherwise.The conclusions in the latter paragraph recover what has long been known for the modulational stability profile of the NLS equation [75], summarized in that stability is prescribed for P Q < 0 and instability occurs otherwise.63) and ( 64) are either both in the stable regime (for k 2 < kr where P 1 Q 11 < 0 and P 2 Q 22 < 0), or the first is in the stable regime and the second in the unstable (for k 2 > kr where P 1 Q 11 < 0 and P 2 Q 22 > 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig. 3(a).

Modulational instability: numerical parametric analysis
We have calculated the growth rate Γ numerically for the CNLS equations by finding the roots of the polynomial in Ω function resulting from Eq. ( 76), and subsequently  identifying the one with the highest imaginary part.That highest imaginary part is plotted as a function of the wavenumber of the perturbation K and the wavenumber of the second carrier wave k 2 for six values of the spectra index κ, i.e., for κ = 2, 3, 4, 5, 6, 100, and fixed value of k 1 and the amplitudes of the nonlinear wave modes A 1,0 and A 2,0 (the latter are fixed throughout the paper to A 1,0 = A 2,0 = 0.15).In Fig. 5, in particular, the growth rate Γ is mapped on the k 2 − K plane for k 1 = 0.4 < k r (κ) and the six values of the spectral index κ mentioned earlier.For that value of k 1 , the first of the decoupled CNLS equations is modulationally stable (P 1 Q 11 < 0), while the second one is also modulationally stable for k 2 < k r (κ) (P 2 Q 22 < 0) while it is modulationally unstable for k 2 > k r (κ) (P 2 Q 22 > 0).The instability boundary for the second decoupled CNLS equation is indicated in Fig. 5 (as well as in Figs.7  and 9 below) by the black-solid vertical line.All the sub-figures Fig. 5 reveal complex instability patterns, in which modulational instability appears both in areas of the k 2 − K plane where P 1 Q 11 < 0 and P 2 Q 22 < 0 (stable-stable regime), as well as in areas where P 1 Q 11 < 0 and P 2 Q 22 > 0 (stable-unstable regime).The growth rate Γ has not been mapped for k 2 < 0.1 on the k 2 − K plane in all sub-figures for clarity, since it acquires very large and even divergent values in this areas.63) and ( 64) are either both in the stable regime (for k 2 < kr where P 1 Q 11 < 0 and P 2 Q 22 < 0), or the first is in the stable regime and the second in the unstable (for k 2 > kr where P 1 Q 11 < 0 and P 2 Q 22 > 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig. 3(a).
Generally speaking, the decrease of the spectral index κ increases the areas on the k 2 − K plane in which MI occurs.Moreover, the maximum values of the growth rate Γ seem to increase considerably with decreasing κ (note the different scales in the colorbars).Note also the difference in the mudulational instability patterns for values of κ in the far-and near-equilibrium region, i.e., for κ = 2 and κ = 3, 4, 5, 6, 100,  , the stable-stable regime for which P 1 Q 11 < 0 and P 2 Q 22 < 0 (k 2 < kr, left from the vertical line) and the stable-unstable regime for which P 1 Q 11 < 0 and P 2 Q 22 > 0 (k 2 > kr, right from the vertical line).The stability situation in these maps is the same as that in figs.5 and 6.Note that kr depends weakly on the spectral index κ, but the differences cannot be observed in this scale.
respectively, which are separated by a rough empirical boundary at around κ ≃ 2.5.Indeed, for κ = 2, the growth rate Γ exhibits very large values, and the CNLS system is modulationally unstable in the larger part of the k 2 − K plane in which the decoupled CNLS equations are in the stable-stable regime (P 1 Q 11 < 0 and P 2 Q 22 < 0).For the next higher value of κ, i.e., for κ = 3, the CNLS system is modulationally stable in the most part of the area in which the decoupled CNLS equations are in the stablestable regime.With further increasing the value of κ, the areas of the k 2 − K plane in which the CNLS system is modulationally unstable gradually shrink slightly, while the higher values of the growth index Γ also gradually decrease.In Fig. 5, as well as in Figs.7 and 9, it is evident that the most significant variation in the patterns of Γ seems to occur for values of the spectral index κ within the strongly nonthermal region, i.e., for 1.5 < κ < 2.5.To the contrary, the variation of the Γ pattern for κ = 3 and up to infinity is rather slight, at least as long as the size of the modulational instability areas on the k 2 − K plane and the magnitude of Γ in these areas are concerned.In Fig. 6, maps of Γ on the k 2 −K plane are shown for more values of κ within the strongly nonthermal region, i.e., for κ = 1.75, 2.25, and 2.50.In the figure, the map of the already shown map for κ = 2 is also shown along with the other three.Inspection of the four subfigures of Fig. 6 readily reveals the dramatic variations of the modulational instability Γ patterns with the variation of κ.For κ = 1.75 (Fig. 6(a)), in particular, modulational instability accompanied by very high values of Γ is apparent over almost the whole of the k 2 − K plane shown.As we shall see also below, this is true for the other case studies presented here.Empirically we find that this the situation for κ up to 1.90, where where some modulationally stable islands begin to appear on the k 2 − K plane.When κ is increased to 2, as we have also discussed in the previous figure, a substantial modulationally stable appear on the k 2 − K plane that extends to both the stable-stable and stable-unstable regions of the decoupled CNLS equations.At the same time, the mudulationally unstable areas, which remain large, are localized to two basic lobes, the one in the stable-stable and the other in the stable-unstable region of the decoupled CNLS equations.In the centers of these lobes, the highest values of Γ occur.With further increasing κ, the modulationally stable area increases against the modulationally unstable ones.Moreover, the strongly modulationally unstable areas gradually concentrates in the stable-unstable region of the decoupled CNLS equations while another strongly modulationally unstable area remains for low k 2 .Remnants of the latter can be even observed in Fig. 5(b) (for κ = 3).
The results shown in Fig. 7 have been obtained with the same parameters as those used in Fig. 5, except that in this case k 1 = 0.7.Thus, as long as the stability of the decoupled CNLS equations is concerned, the situation is exactly the same as that in Fig. 5, i.e., the decoupled CNLS equations are in the stable-stable regime (P 1 Q 11 < 0 and P 2 Q 22 < 0) for k 2 < k r (κ), while they are in the stable-unstable regime for k 2 > k r (κ) (P 1 Q 11 < 0 and P 2 Q 22 > 0).In comparison with Fig. 5, we observe that for κ = 3, 4, 5, 6, 100, areas of strong instability with relatively high values of the growth rate Γ occur at low k 2 for most of the perturbation wavenumbers shown in these figures.Again, these strongly unstable areas shrink slightly with increasing spectral index from κ = 3 to 100.A similar situation can be described for κ = 2, where high values of the growth rate Γ occur at low k 2 as in the case of higher κ; however, the values of Γ are significantly higher in this case (note again the difference in the scales of colorbars).Also, in a large part of the stable-stable area of the k 2 − K plane modulational instability with moderate values of Γ occurs.Generally speaking, the strength of the modulational instability is weaker and the areas where that instability occurs are smaller in Fig. 7 for k 1 = 0.7 than in Fig. 5 for k 1 = 0.4 (except at low k 2 ).This can be observed more clearly for κ ≥ 3.
As above, we also present maps of Γ for a few more values of κ in the strongly nonthermal region, i.e., in the interval 1.5 < κ < 2.5 in Fig. 8. Roughly, the remarks made for the patterns of Γ in Fig. 6 also hold here.For κ = 1.75 (Fig. 8(a)), almost all of the shown area of the k 2 − K plane is modulationally unstable, with the growth rate achieving very high values of the same order as those in Fig. 6(a).Again, with increasing κ, the modulationally stable areas on the k 2 − K plane grow to the expense of the modulationally unstable ones.The pattern of Γ in Fig. 8(b) is similar to that in in Fig. 6(b) with the two modulationally unstable lobes and the strongly modulationally unstable area for low k 2 .In Fig. 8(c), the two lobes disappear, while a strongly  63) and ( 64) are either both in the unstable regime (for k 2 > kr where P 1 Q 11 > 0 and P 2 Q 22 > 0), or the first is in the unstable regime and the second in the stable (for k 2 < kr where P 1 Q 11 > 0 and P 2 Q 22 < 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig. 3(a).
modulationally unstable area at low k 2 remains, along with another modulationally unstable area which is located in the stable-unstable region of the decoupled CNLS equations, i.e., in the area with k 2 > k r (κ).Peculiarly, the modulationally unstable areas in Fig. 8(d) are larger to those in Fig. 8(c).However, the general tendency of shrinking of the modulationally unstable areas with increasing κ continuous for  κ > 2.5.Note that in this case the strongly modulationally unstable area for low k 2 persist, although becoming much smaller with increasing κ, even for very large values of κ, i.e., κ = 100 (Fig. 7(f)).The results shown in Fig. 9 have been obtained with the same parameters as those used in Figs. 5 and 7, except that the value of k 1 is now k 1 = 1.6 > k r (κ).That means that the decoupled CNLS equations are either in the unstable-stable regime (P 1 Q 11 > 0 and P 2 Q 22 < 0) for k 2 < k r (κ), i.e., to the left from the black-solid vertical line, while they are in the unstable-unstable regime (P 1 Q 11 > 0 and P 2 Q 22 > 0) for k 2 > k r (κ), i.e., to the right of the black-solid vertical line.In this case, we observe a change in the obtained instability patterns, especially for κ ≥ 3.For all values of κ, we again observe strong modulational instability with high values of Γ for low k 2 .For κ = 2, modulational instability occurs in most parts of the areas of the k 2 − K plane in which the decoupled CNLS equations are either in the unstable-stable or the unstableunstable regime.For κ ≥ 3 we again observe some shrinkage of the instability areas with increasing κ.We can also observe for these values of κ that the modulationally unstable areas are mostly in the unstable-stable regime than in the unstable-unstable regime of the decoupled CNLS equations.
In the case shown in Fig. 10, also, the patterns of the growth rate Γ on the k 2 − K plane for a few values of κ in the strongly nonthermal region, i.e. for small values of κ.The evolution of these patterns is similar to the corresponding evolution described in Fig. 10.Thus, from Figs. 5-10 we can safely conclude that the most dramatic variation in the onset (region) of modulational instability on the k 2 − K plane occurs in the interval of κ values within the strongly nonthermal region.
Fig. 11 The growth rate Γ of the CNLS system (black curves) and the growth rates Γ 1 and Γ 2 (red and green curves, respectively) of the first and the second decoupled CNLS Eqs. ( 63) and (64 as a function of the wavenumber of the perturbation K for κ = 3, Ψ 1,0 = Ψ 2,0 = 0.15, and (a) k 1 = 1 and k 2 = 0.6, for which P 1 Q 11 < 0 and Representative examples of growth rate variation.It may be interesting to present some illustrative examples for the growth rate Γ as a function of the perturbation wavenumber K which may be distinctly different from those one may observe in a single NLS equation.In Fig. 11 Γ(K) curves are shown in three cases that differ in respect to the stability regime of the decoupled CNLS equations, for κ = 3.Along with these curves, we also show the corresponding growth rates Γ 1 and Γ 2 obtained for the first and the second decoupled CNLS equation, respectively.In Fig. 11(a), in which k 1 = 1 < k r (κ) and k 2 = 0.6 < k r (κ), i.e., in the stable-stable regime of the decoupled CNLS equations, we obtain a Γ(K) curve with two instability lobes.Interestingly, stronger modulational instability occurs in the second lobe, for K between 3 and 4, than in the first lobe for K between 0 and 0.5 (low K).Since both the decoupled equations are modulationally stable, the corresponding growth rates Γ 1 and Γ 2 are both zero.In Fig. 11(b), in which k 1 = 1.6 > k r (κ) and k 2 = 1.35 > k r (κ), i.e., in the unstable-stable regime of the decoupled CNLS equations, we obtain a Γ(K) curve with two humps.Moreover, we obtain a non-zero Γ 1 (κ) curve which has the commonly observed shape.A comparison between the two indicates that besides the difference in shape, the values of Γ 1 (κ) are always significantly smaller than those of Γ(κ).Moreover, the modulationally unstable interval in K is much smaller in the first decoupled NLS equation than that of the coupled CNLS system.In Fig. 11(c), in which k 1 = 1.6 > k r (κ) and k 2 = 1.45 > k r (κ) = 1.42, i.e., in the unstable-unstable regime of the decoupled CNLS equations, we obtain three curves for the growth rates Γ(K), Γ 1 , and Γ 2 which have the same shape but very different size.As a result, the interval of the perturbation wavenumber K in which the coupled CNLS system is modulationally unstable, is much larger that the corresponding one of either the first or the second decoupled NLS equation.

Discussion and conclusions
A pair of CNLS equations was derived for two co-propagating and interacting electrostatic wavepackets in a plasma fluid, comprising cold inertial ions evolving against a thermalized, kappa-distributed electron background.A multiple-scale technique has been adopted, similar to the well known Newell's method in nonlinear optics, for the reduction of the original plasma fluid model to the pair of CNLS equations which is generally non-integrable since the dispersion, nonlinearity, and coupling coefficients P j and Q ij (j = 1, 2) depend on the carrier wavenumbers k 1 and k 2 (both assumed arbitrary) of the two interacting wavepackets.The exact mathematical expressions for those coefficients have been analytically obtained, and in fact depend parametrically on the plasma parameters too.Although P j for the particular model considered here is negative, while it approaches asymptotically zero for large k j , the coefficients Q jj can be either positive or negative.The value of k j that separates positive from negative Q jj values, k r , depends on the value of the spectral index κ which characterizes the kappa distribution of the electrons, i.e., k r = k r (κ).The spectral index enters into the expressions of Q jj through the coefficients c 1 , c 2 , and c 3 of the Mc Laurin expansion of that distribution.Interestingly, the κ−dependence of k r is not monotonous but it exhibits a minimum around κ ≃ 2.5 which has been set as a rough boundary between electrons far-and near-equilibrium by empirical investigations of Space plasmas.The k r (κ) dependence is weak for most values of κ, while it becomes stronger for κ < 2 where it increases abruptly until its divergence at κ = 3/2.For very large values of κ, on the other hand, k r converges slowly to 1.47 from below.
The coupling coefficients Q 12 and Q 21 , which depend on both wavenumbers k 1 and k 2 , may also be either positive or negative.While P j and Q jj have values of the order of unity, the values of Q 12 and Q 21 may also become very large positive or negative, implying very strong coupling for the interacting electrostatic wavepackets in the plasma fluid.All this prescribes a multifaceted dynamical profile with strong interconnection among the various physical parameters.A compatibility condition was derived through standard modulational instability analysis, in the form of a fourth degree polynomial in the frequency of the perturbation.By numerically calculating the roots of the compatibility condition and identifying the one with the largest imaginary part (growth rate), modulationally stable and unstable areas on the k 2 − K plane have been identified.The instability growth rate Γ, mapped on that plane, reveals that modulational instability occurs in areas of the k 2 − K plane (say, for fixed k 1 ), actually for all four stability regimes of the decoupled CNLS equation (i.e., those resulting by setting Q 12 = Q 21 = 0), for all values of κ used here.
Modulational instability in fact may occur even if both waves of the decoupled system are stable (stable-stable regime), since their strong nonlinear coupling may potentially destabilize the system in regions with large perturbation wavenumbers K.The instability that results from the nonlinear coupling between the two wavepackets is more intense in the case of equal amplitudes of the two waves (as those used in the present work); instability areas also occupy larger areas ("windows") on the k 2 − K plane in the latter case.When one or both the waves of the decoupled system are unstable, their nonlinear coupling leads again to an unstable pair of wavepackets.A comparison of the growth rates of the two waves of the decoupled system with that of the waves in the coupled one reveals that the latter is always larger than the former; furthermore, the unstable areas on the k 2 − K plane are larger than those of the waves in the decoupled system.Also, as it becomes obvious from Figs. 5-10, the most dramatic variation of the Γ patterns on the k 2 − K plane occurs for values of the spectral index κ withing the strongly nonthermal region, i.e., in the interval 1.5 < κ < 2.5.Moreover, as the numerical calculation of the growth rate Γ indicate, the CNLS equations are modulationally unstable in almost all of the k 2 − K plane shown in the figures below κ ∼ 1.9.For these values of κ, i.e., for κ < 1.9, the values of Γ increase fast with decreasing κ and diverge in the limit of κ → 1.5.
Extensive numerical search support the conclusion that the nonlinear wave modes of the CNLS equations are apparently always partly modulationally unstable, i.e., that there is an interval in the perturbation wavenumber K for which modulational instability occurs.Although it cannot be considered as a definite proof our investigations shows that it is highly possible.As mentioned above, the instability is probably due to the large values of the coupling coefficients Q 12 and Q 21 (as compared with the values of the coefficients P j and Q jj ).Note that it is actually the product of the two, i.e., Q 12 Q 21 which enters into Ω 4 c that appears in the right-hand-side of the compatibility condition Eq. ( 76).However, it is possible to have modulationally stable nonlinear wave modes by "eliminating" Ω 4 c from the right-hand-side of Eq. (76).This is possible if we set one of the two coeffcients Q 12 and Q 21 equal to zero (in the considered model it is not possible to find k 1 and k 2 as both these coefficients become zero).In fact, Q 12 is for example, zero along the curves shown in Fig. 3(b) or as the contours shown in 4 on the k 1 − k 2 plane.Along these curves Q 12 = 0 so that the CNLS equations are still at least partially coupled, but their perturbation fremquencies can be calculated from Eq. (80) which are the same as for the two decoupled CNLS equations.Thus, for modulational stability in this special case we can assert that P 1 Q 11 < 0 and P 2 Q 2 < 0 simultaneously.This condition can be clearly satisfied within the considered model for k 1 , k 2 < k r (κ).A bare look at Figs. 3(b) and 4 is enough to see that there are available k 1 and k 2 that satisfy the latter inequality.
The occurence of very broad modulational instability areas in the parameter space of k 1 , k 2 , and K, may on the other hand favor the emergence of vector solitons which may emerge in combinations of bright and dark components.The possible existence of such vector solitons, the determination of their parameters such as their width and amplitudes, as well as their stability, and how these are affected by the spectral index κ, will be the focus of future work.

Acknowledgements
Authors IK and NL gratefully acknowledge financial support from Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates, via the project CIRA-2021-064 (8474000412).IK also acknowledges financial support from KU via the project FSU-2021-012 (8474000352) as well as from KU Space and Planetary Science Center, Abu Dhabi, United Arab Emirates, via grant No. KU-SPSC-8474000336.
This work was completed during a long research visit by author IK to the National and Kapodistrian University of Athens, Greece.During the same period, IK also held an Adjunct Researcher status at the Hellenic Space Center, Greece.The hospitality of both hosts, represented by Professor D.J. Frantzeskakis and Professor I. Daglis, respectively, is warmly acknowledged.where j = 1, 2, and k j , ℓ j defined as ℓ 1 = 1, ℓ 2 = −1, k 1 = 2, and k 2 = 1.
-leads us to an equation in the form:

Fig. 3 (
Fig. 3 (a) The root kr of Q 11 as a function of the spectral index κ.For large values of κ, that root converges asymptotically to the value k 1 = kr = 1.47, i.e., its value in the Maxwell-Boltzmann case.The inset shows a magnified part of the main figure for physically relevant values of κ.The orange-dashed line is located at kr = 1.47.(b) The zeroes of the coupling coefficient Q 12 on the k 1 − k 2 plane are shown (in ascending order, top to bottom) for κ = 2 (black), κ = 3 (red), κ = 4(green), κ = 5 (blue), κ = 6 (orange), κ = 100 (brown).The coupling coefficient Q 12 is negative in the area of the k 1 − k 2 plane above the corresponding curve (and positive underneath).Note that the coupling coefficient Q 21 can be obtained from Q 12 by interchanging k 1 and k 2 in the corresponding expressions.The corresponding curves for the zeros of Q 21 can be obtained by a reflection with respect to the diagonal k 1 = k 2 (and thus Q 21 would be positive in the area above the curves, and negative underneath).
3(b), in which Q 12 is negative in the areas of the k 1 − k 2 plane above the boundary.It is obvious from the above Figure that smaller values of the 2nd wavenumber (k 2 )

Fig. 4
Fig. 4 The nonlinear coupling coefficient Q 12 mapped on the plane of the wavenumbers of the carrier waves k 1 and k 2 , for (a) κ = 2, (b) κ = 3, (c) κ = 4, (d) κ = 5, (e) κ = 6, (f) κ = 100.In (f), the κ distribution practically coincides with the Maxwell-Boltzmann distribution.The black contour indicates the boundary between positive and negative values of Q 12 .Note that the plots above are shown from k 2 = 0.1 to 2 when comparing the boundaries (the contours) with the curves in Fig. 3(b).At low k 2 values, Q 12 becomes very large and has been omitted, for clarity.The corresponding coupling coefficients Q 21 are obtained upon a reflection around the line k 1 = k 2 (i.e., upon a permutation between k 1 and k 2 .).

Fig. 5
Fig. 5 Maps of the growth rate Γ on the plane of the wavenumber of the second carrier wave k 2 and the wavenumber of the perturbation K for k 1 = 0.4, Ψ 1,0 = Ψ 2,0 = 0.15, and (a) κ = 2; (b) κ = 3; (c) κ = 4; (d) κ = 5; (e) κ = 6; (f) κ = 100.For the selected value of k 1 , the two single NLS equations resulting from decoupling the CNLS equations (e.g., by setting Q 12 = Q 21 = 0 in Eqs.(63) and (64) are either both in the stable regime (for k 2 < kr where P 1 Q 11 < 0 and P 2 Q 22 < 0), or the first is in the stable regime and the second in the unstable (for k 2 > kr where P 1 Q 11 < 0 and P 2 Q 22 > 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig.3(a).

Fig. 6
Fig. 6 Maps of the growth rate Γ on the k 2 − K plane for (a) κ = 1.75;(b) κ = 2.00; (c) κ = 2.25; (d) κ = 2.50.The other parameters as in Fig.5.The vertical black-solid line indicates the boundary between two different regimes of the decoupled CNLS equations, i.e., the stable-stable regime for which P 1 Q 11 < 0 and P 2 Q 22 < 0 (k 2 < kr, left from the vertical line) and the stable-unstable regime for which P 1 Q 11 < 0 and P 2 Q 22 > 0 (k 2 > kr, right from the vertical line).Note that kr depends weakly on the spectral index κ, but the differences cannot be observed in this scale.

Fig. 7
Fig. 7 Maps of the growth rate Γ on the plane of the wavenumber of the second carrier wave k 2 and the wavenumber of the perturbation K for k 1 = 0.7, Ψ 1,0 = Ψ 2,0 = 0.15, and (a) κ = 2; (b) κ = 3; (c) κ = 4; (d) κ = 5; (e) κ = 6; (f) κ = 100.For the selected value of k 1 , the two single NLS equations resulting from decoupling the CNLS equations (e.g., by setting Q 12 = Q 21 = 0 in Eqs.(63) and (64) are either both in the stable regime (for k 2 < kr where P 1 Q 11 < 0 and P 2 Q 22 < 0), or the first is in the stable regime and the second in the unstable (for k 2 > kr where P 1 Q 11 < 0 and P 2 Q 22 > 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig.3(a).

Fig. 8
Fig. 8 Maps of the growth rate Γ on the k 2 − K plane for (a) κ = 1.75;(b) κ = 2.00; (c) κ = 2.25; (d) κ = 2.50.The other parameters as in Fig.7.The vertical black-solid line indicates the boundary between two different regimes of the decoupled CNLS equations, i.e., the stable-stable regime for which P 1 Q 11 < 0 and P 2 Q 22 < 0 (k 2 < kr, left from the vertical line) and the stable-unstable regime for which P 1 Q 11 < 0 and P 2 Q 22 > 0 (k 2 > kr, right from the vertical line).The stability situation in these maps is the same as that in figs.5 and 6.Note that kr depends weakly on the spectral index κ, but the differences cannot be observed in this scale.

Fig. 9
Fig. 9 Maps of the growth rate Γ on the plane of the wavenumber of the second carrier wave k 2 and the wavenumber of the perturbation K for k 1 = 1.6, Ψ 1,0 = Ψ 2,0 = 0.15, and (a) κ = 2; (b) κ = 3; (c) κ = 4; (d) κ = 5; (e) κ = 6; (f) κ = 100.For the selected value of k 1 , the two single NLS equations resulting from decoupling the CNLS equations (e.g., by setting Q 12 = Q 21 = 0 in Eqs.(63) and (64) are either both in the unstable regime (for k 2 > kr where P 1 Q 11 > 0 and P 2 Q 22 > 0), or the first is in the unstable regime and the second in the stable (for k 2 < kr where P 1 Q 11 > 0 and P 2 Q 22 < 0).Note that kr depends on the particular value of the spectral index κ in accordance with the inset of Fig.3(a).

Fig. 10
Fig. 10 Maps of the growth rate Γ on the k 2 − K plane for (a) κ = 1.75;(b) κ = 2.00; (c) κ = 2.25; (d) κ = 2.50.The other parameters as in Fig.9.The vertical black-solid line indicates the boundary between two different regimes of the decoupled CNLS equations, i.e., the unstable-stable regime for which P 1 Q 11 > 0 and P 2 Q 22 < 0 (k 2 < kr, left from the vertical line) and the unstable-unstable regime for which P 1 Q 11 > 0 and P 2 Q 22 > 0 (k 2 > kr, right from the vertical line).Note that kr depends weakly on the spectral index κ, but the differences cannot be observed in this scale.