[1]N. Lazarides

[1]Department of Mathematics, Khalifa University of Science and Technology,
Abu Dhabi, P.O. Box 127788, United Arab Emirates

2]Space and Planetary Science Center, Khalifa University of Science and
Technology, Abu Dhabi, P.O. Box 127788, United Arab Emirates

3]Hellenic Space Center, Leoforos Kifissias 178, Athens, GR-15231, Greece

4]Department of Physics, National and Kapodistrian University of Athens,
Zografou, Athens, GR-15784, Greece

# Coupled Nonlinear Schrödinger (CNLS) Equations for two interacting
electrostatic wavepackets in a non-Maxwellian fluid plasma model

<EMAIL_ADDRESS>Ioannis Kourakis<EMAIL_ADDRESS>*
[  [  [

###### Abstract

The nonlinear dynamics of two co-propagating electrostatic wavepackets in a
one-dimensional non-magnetized plasma fluid model is considered, from first
principles. The coupled waves are characterized by different (carrier)
wavenumbers and amplitudes. A plasma consisting of non-thermalized
($\kappa-$distributed) electrons evolving against a cold (stationary) ion
background is considered. The original model is reduced, by means of a
multiple-scale perturbation method, to a pair of coupled nonlinear Schrödinger
(CNLS) equations for the dynamics of the wavepacket envelopes. For arbitrary
wavenumbers, the resulting CNLS equations exhibit no known symmetry and thus
intrinsically differ from the Manakov system, in general.

Exact analytical expressions have been derived for the dispersion, self-
modulation (nonlinearity) and cross-modulation (coupling) coefficients
involved in the CNLS equations, as functions of the wavenumbers ($k_{1}$,
$k_{2}$) and of the spectral index $\kappa$ characterizing the electron
profile. An analytical investigation has thus been carried out of the
modulational instability (MI) properties of this pair of wavepackets, focusing
on the role of the intrinsic (variable) parameters. Modulational instability
is shown to occur in most parts of the parameter space. The instability
window(s) and the corresponding growth rate are calculated numerically 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$ index (values) suggest that deviation from thermal
(Maxwellian) equilibrium, for low $\kappa$ values, leads to enhances 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
whose coefficients present no obvious symmetries for arbitrary $k_{1}$ and
$k_{2}$. Although we have focused on electrostatic wavepacket propagation in
nonthermal (non-Maxwellian) plasma, the results of this study 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.

###### keywords:

Plasma fluid model, Multiscale perturbative reduction, Coupled nonlinear
Schrödinger equations, Modulational instability, Kappa distribution

##  1 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 Maxwellian-like “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 Korteweg - de 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].

The consideration of two co-propagating wavepackets in a particular plasma
fluid and their mutual interaction, results in a pair of coupled NLS (CNLS)
equations, whose coefficients generally depend on the carrier wavenumber(s) of
the respective carrier waves. Electron plasma (Langmuir) waves and ion
acoustic waves [29, 30, 31], amidst other studies that focused mostly on
electromagnetic modes[32, 33, 34, 35, 36, 37, 38]. Beyond plasma science,
general forms of CNLS equations have been investigated in recent years with
respect to the existence of vector solitons, breathers, and rogue wave
solutions, in contexts including higher-order coupled NLS systems [39, 40],
coupled mixed-derivative NLS equations [41, 42, 43], variable coefficients
CNLS equations [44], non-autonomous CNLS equations [45, 46, 47], systems
involving four-wave mixing terms [48], coupled cubic-quintic NLS equations
[49], space-shifted CNLS equations [50] and non-autonomous partially non-local
CNLS equations [51], among others.

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 left-handed 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 $\kappa$). 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 $\kappa$ 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 $\kappa$ 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 $\kappa$ (besides the externally
determined $k_{1}$ and $k_{2}$). This constitutes advantage to other models
since it isolates the effect of $\kappa$ on the investigated MI properties.
Moreover, a thorough MI analysis is performed as a function of $\kappa$ 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 $\kappa$ (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

###  2.1 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]

$\displaystyle\frac{\partial n}{\partial t}+\frac{\partial(nu)}{\partial
x}=0,$ (1) $\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial
u}{\partial x}=-\frac{\partial\phi}{\partial x},$ (2)
$\displaystyle\frac{\partial^{2}\phi}{\partial x^{2}}=n_{e}-n,$ (3)

respectively, with $\phi$ 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 $\omega_{pi}^{-1}$ and to the Debye screening
length $\lambda_{D}$, respectively, where

$\displaystyle\lambda_{D}=\left(\frac{\varepsilon_{0}k_{B}T_{e}}{z_{i}n_{i0}e^{2}}\right)^{1/2},$
(4)
$\displaystyle\omega_{pi}=\left[\frac{n_{0,i}(z_{i}e)^{2}}{\varepsilon_{0}m_{i}}\right]^{1/2}\,.$
(5)

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

$c_{s}=\lambda_{D}\,\omega_{pi}=\left(\frac{z_{i}k_{B}T_{e}}{m_{e}}\right)^{1/2},$
(6)

i.e. essentially the (ion-acoustic) sound speed. Finally, the electrostatic
potential $\phi$ is normalized to $k_{B}T_{e}/e$. Adopting a standard
notation, the symbols $\varepsilon_{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]

$n_{e}=\left(1-\frac{\phi}{\kappa-\frac{3}{2}}\right)^{-\kappa+\frac{1}{2}}.$
(7)

Note that this function is characterized by a real parameter $\kappa$ (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 $\kappa$ in each
particular case. Possible values of $\kappa$ in Eq. (7) lie in the interval
$(3/2,\infty)$, while for $\kappa=+\infty$ the kappa distribution reverts to
the Maxwell-Boltzmann one. Lower values of $\kappa$ 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 $\kappa$ into a near- ($3/2<\kappa<2.5$) and a far-
($2.5<\kappa<\infty$) equilibrium region. In particular, kappa distributions
with $2<\kappa<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 $\kappa$ evenly distributed in the interval $2<\kappa<6$ as
well as a large value of $\kappa=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 $\phi$, i.e. for
$|\phi|\ll 1$ (in normalized units), the kappa distribution Eq. (6) can be
expanded in a Mc Laurin series in powers of $\phi$ as

$n_{e}\approx 1+c_{1}\phi+c_{2}\phi^{2}+c_{3}\phi^{3}+\cdots\,,$ (8)

from which we shall keep terms up to order $\phi^{3}$. In that case, the first
three expansion coefficients are given by

$\displaystyle c_{1}$ $\displaystyle=$
$\displaystyle\frac{\kappa-\frac{1}{2}}{\kappa-\frac{3}{2}}\,,$ (9)
$\displaystyle c_{2}$ $\displaystyle=$
$\displaystyle\frac{\left(\kappa-\frac{1}{2}\right)\left(\kappa+\frac{1}{2}\right)}{2!\left(\kappa-\frac{3}{2}\right)^{2}}\,,$
(10) $\displaystyle c_{3}$ $\displaystyle=$
$\displaystyle\frac{\left(\kappa-\frac{1}{2}\right)\left(\kappa+\frac{1}{2}\right)\left(\kappa+\frac{3}{2}\right)}{3!\left(\kappa-\frac{3}{2}\right)^{3}}\,.$
(11)

###  2.2 Application of the reductive perturbation method

We introduce fast and slow variables so that the derivatives are approximated
as

$\displaystyle\frac{\partial}{\partial
t}\rightarrow\sum_{k=0}^{\infty}\varepsilon^{k}\frac{\partial}{\partial
t_{k}},~{}~{}~{}\frac{\partial}{\partial
x}\rightarrow\sum_{k=0}^{\infty}\varepsilon^{k}\frac{\partial}{\partial
x_{k}}$ (12)

where $\varepsilon$ is a small quantity, and $x_{k}=\varepsilon^{k}x$,
$t_{k}=\varepsilon^{k}t$ (for $k=0,1,2,3,...$). The dependent variables $n$,
$u$, and $\phi$ can be expanded around their equilibrium value as

$\displaystyle n=1+\sum_{k=1}^{\infty}\varepsilon^{k}n_{k},$ (13)
$\displaystyle u=\sum_{k=1}^{\infty}\varepsilon^{k}u_{k},$ (14)
$\displaystyle\phi=\sum_{k=1}^{\infty}\varepsilon^{k}\phi_{k}.$ (15)

The expansions of $\\{n,u,\phi\\}$ 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 $\varepsilon$ up to order
$\varepsilon^{3}$.

The fluid equations obtained at first order ($\propto\varepsilon^{1}$) read

$\displaystyle\frac{\partial n_{1}}{\partial
t_{0}}+\frac{\partial(u_{1})}{\partial x_{0}}=0,$ (16)
$\displaystyle\frac{\partial u_{1}}{\partial
t_{0}}+\frac{\partial\phi_{1}}{\partial x_{0}}=0,$ (17)
$\displaystyle\frac{\partial^{2}\phi_{1}}{\partial
x_{0}^{2}}+n_{1}-c_{1}\phi_{1}=0,$ (18)

where $x_{0}=x$ and $t_{0}=t$. In these equations, the following Ansatz will
be introduced

$\displaystyle
S_{1}=S_{1,1}^{(1)}\,e^{+i\theta_{1}}+S_{1,2}^{(1)}\,e^{+i\theta_{2}}+c.c.,$
(19)

where $S_{1}=\\{n_{1},u_{1},\phi_{1}\\}$ and the phases are
$\theta_{j}=k_{j}x-\omega_{j}t$ with $k_{j}$ and $\omega_{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}^{(1)}$ –
– (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$)

$\displaystyle-\omega_{j}n_{1,j}^{(1)}+k_{j}u_{1,j}^{(1)}=0,$
$\displaystyle-\omega_{j}u_{1,j}^{(1)}+k_{j}\phi_{1,j}^{(1)}=0,$ (20)
$\displaystyle+n_{1,j}^{(1)}-(k_{j}^{2}+c_{1})\phi_{1,j}^{(1)}=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\times 3$ systems with unknowns
$\\{n_{1,j}^{(1)},u_{1,j}^{(1)},\phi_{1,j}^{(1)}\\}$ for $j=1$ and $j=2$,
respectively, are obtained as ($j=1,2$)

$\displaystyle\begin{bmatrix}-i\omega_{j}&+ik_{j}&0\\\
0&-i\omega_{j}&+ik_{j}\\\
1&0&-(k_{j}^{2}+c_{1})\end{bmatrix}\begin{bmatrix}n_{1,j}^{(1)}\\\
u_{1,j}^{(1)}\\\ \phi_{1,j}^{(1)}\end{bmatrix}=\begin{bmatrix}0\\\ 0\\\
0\end{bmatrix},$ (21)

For these systems to have a non-trivial solution, it is required that their
determinants should vanish, i.e.,

$\displaystyle D_{j}=\begin{vmatrix}-\omega_{j}&+k_{j}&0\\\
0&-\omega_{1}&+k_{j}\\\ 1&0&-(k_{j}^{2}+c_{1})\end{vmatrix}=0,$ (22)

which provides the linear frequency dispersion relation ($\omega$) and the
corresponding group ($v_{g}$) and phase ($v_{ph}$) velocity for the $j-$th
wave ($j=1$ or 2) as

$\displaystyle\omega_{j}$ $\displaystyle=$
$\displaystyle\frac{k_{j}}{[k_{j}^{2}+c_{1}(\kappa)]^{1/2}},$ (23)
$\displaystyle v_{g,j}$ $\displaystyle=$
$\displaystyle\frac{c_{1}(\kappa)}{[k_{j}^{2}+c_{1}(\kappa)]^{3/2}},$ (24)
$\displaystyle v_{ph,j}$ $\displaystyle=$
$\displaystyle\frac{\omega_{j}}{k_{j}}=\frac{1}{[k_{j}^{2}+c_{1}(\kappa)]^{1/2}},$
(25)

where the dependence of the expansion coefficient $c_{1}$ on the spectral
index $\kappa$ is explicitly indicated. The systems of Eqs. (21) can be solved
in terms of the variables ($j=1,2$)

$\phi_{1,j}^{(1)}=\Psi_{j}$ (26)

to give

$n_{1,j}^{(1)}=\left(\frac{k_{j}}{\omega_{j}}\right)^{2}\Psi_{j},~{}~{}~{}~{}u_{1,j}^{(1)}=\frac{k_{j}}{\omega_{j}}\Psi_{j}.$
(27)

Figure 1:  The frequency dispersion $\omega_{1}$ (a), the corresponding group
velocity $v_{g,1}$ (b) as calculated from Eq. (24), and the phase velocity
$v_{ph,1}=\omega_{1}/k_{1}$ (c) as calculated from Eq. (25), as a function of
the wavenumber $k_{1}$ for $\kappa=2$ (black), $\kappa=3$ (red), $\kappa=4$
(green), $\kappa=5$ (blue), $\kappa=6$ (orange), $\kappa=100$ (brown).

The frequency dispersion $\omega_{1}$ and the group velocity $v_{g,1}$ which
have been calculated 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 $\kappa$. The
frequency dispersion curves seem to have the same qualitative features for any
value of $\kappa$. However, the frequency (or fixed $k_{1}$) increases with
increasing $\kappa$ for any $k_{1}$ in the (physically relevant) interval
shown, while it saturates for large spectral index $\kappa$ to its Maxwell-
Boltzmann value. Note that the curve obtained for $\kappa=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}\lesssim 0.9$) the group velocity
$v_{g,1}$ increases with increasing $\kappa$ just as the frequency does in
Fig. 1(a). However, for high $k_{1}$ (i.e., $k_{1}\gtrsim 0.9$) that tendency
is reversed, and the group velocity $v_{g,1}$ decreases with increasing
$\kappa$. Note that for very low $k_{1}$ (i.e., for $k_{1}\rightarrow 0$), the
group velocity increases from $v_{g,1}\simeq 0.59$ for $\kappa=2$ to
$v_{g,1}=1$ for $\kappa=100$ (practically the same as in the case of Maxwell-
Boltzmann distributed electrons, as expected for large values of $\kappa$).
This constitutes a relative increase of up to $\simeq 70\%$ in group velocity
(value) between the extreme values of $\kappa$ adopted here. We draw the
conclusion that lower values of $\kappa$, 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. $\gtrsim$ the Debye length, more or less). The
inverse trend is observed for short (sub-screening-radius) wavelengths (i.e.
for wavenumber exceeding unity, roughly), since lower $\kappa$ 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 $\kappa$ 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 $\kappa$, 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_{k_{1}\rightarrow 0}v_{g,1}=\lim_{k_{1}\rightarrow
0}v_{ph,1}=1/\sqrt{c_{1}(\kappa)}$, 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\rightarrow 2$.

The equations obtained at the second order in $\varepsilon$ are

$\displaystyle\frac{\partial n_{2}}{\partial
t_{0}}+\frac{\partial(u_{2})}{\partial x_{0}}={\cal F}_{1},$
$\displaystyle\frac{\partial u_{2}}{\partial
t_{0}}+\frac{\partial\phi_{2}}{\partial x_{0}}={\cal F}_{2},$ (28)
$\displaystyle\frac{\partial^{2}\phi_{2}}{\partial
x_{0}^{2}}-c_{1}\phi_{2}+n_{2}={\cal F}_{3},$

where ${\cal F}_{1}$, ${\cal F}_{2}$, and ${\cal F}_{3}$ are known functions
of $\Psi_{j}$ which are given in Appendix B. To solve Eqs. (28), we use the
following ansatz

$\displaystyle
S_{2}=S_{2}^{(0)}+\left[S_{2,1}^{(1)}e^{i\theta_{1}}+S_{2,1}^{(2)}e^{2i\theta_{1}}+S_{2,2}^{(1)}e^{i\theta_{2}}+S_{2,2}^{(2)}e^{2i\theta_{2}}\right.$
$\displaystyle\left.+S_{2,+}^{(1)}e^{i(\theta_{1}+\theta_{2})}+S_{2,-}^{(1)}e^{2i(\theta_{1}-\theta_{2})}+c.c.\right],$
(29)

where $S=n,u$ or $\phi$; by “$c.c$” we have denoted the complex conjugate of
the expression within the curly brackets. Upon substitution of Eqs. (2.2) into
Eqs. (28), we obtain the systems of equations

$\displaystyle-\omega_{j}n_{2,j}^{(1)}+k_{j}u_{2,j}^{(1)}=\mu_{1,j},$
$\displaystyle-\omega_{j}u_{2,j}^{(1)}+k_{j}\phi_{2,j}^{(1)}=\mu_{2,j},$ (30)
$\displaystyle-(k_{j}^{2}+c_{1})\phi_{2,j}^{(1)}+n_{2,j}^{(1)}=\mu_{3,j},$

resulting from all terms proportional to $e^{i\theta_{1}}$ and
$e^{i\theta_{2}}$, respectively (to be discussed later; note that $\mu_{1,j}$,
$\mu_{2,j}$, and $\mu_{3,j}$, involving temporal and spatial derivatives of
$\Psi_{j}$, for $j=1,2$, are given in Appendix B),

$\displaystyle-2\omega_{j}n_{2,j}^{(2)}+2k_{j}u_{2,j}^{(2)}=-2k_{j}n_{1,j}^{(1)}u_{1,j}^{(1)},$
$\displaystyle-2\omega_{j}u_{2,j}^{(2)}+2k_{j}\phi_{2,j}^{(2)}=-k_{j}u_{1,j}^{(1)}{{}^{2}},$
(31)
$\displaystyle-(4k_{j}^{2}+c_{1})\phi_{2,j}^{(2)}+n_{2,j}^{(2)}=c_{2}\Psi_{j}^{2},$

for $j=1$ and $2$, the latter resulting from the terms proportional to
$e^{2i\theta_{1}}$ and $e^{2i\theta_{2}}$, respectively. The second harmonic
amplitudes are given by solutions in the form

$\displaystyle n_{2,j}^{(2)}=C_{n,2,j}^{(2)}\Psi_{j}^{2},$ $\displaystyle
u_{2,j}^{(2)}=C_{u,2,j}^{(2)}\Psi_{j}^{2},$ (32)
$\displaystyle\phi_{2,j}^{(2)}=C_{\phi,2,j}^{(2)}\Psi_{j}^{2},$

where the coefficients $C_{n,2,j}^{(2)}$, $C_{u,2,j}^{(2)}$, and
$C_{\phi,2,j}^{(2)}$ are given in Appendix A, and finally

$\displaystyle-(\omega_{1}+\omega_{2})n_{2,+}^{(1)}+(k_{1}+k_{2})u_{2,+}^{(1)}=f_{1}^{(+)}\Psi_{1}\Psi_{2},$
(33)
$\displaystyle-(\omega_{1}+\omega_{2})u_{2,+}^{(1)}+(k_{1}+k_{2})\phi_{2,+}^{(1)}=f_{2}^{(+)}\Psi_{1}\Psi_{2},$
(34)
$\displaystyle-\left[(k_{1}+k_{2})^{2}+c_{1}\right]\phi_{2,+}^{(2)}+n_{2,+}^{(1)}=f_{3}^{(+)}\Psi_{1}\Psi_{2},$
(35)

and

$\displaystyle-(\omega_{1}-\omega_{2})n_{2,-}^{(1)}+(k_{1}-k_{2})u_{2,-}^{(1)}=f_{1}^{(-)}\Psi_{1}\Psi_{2}^{\star},$
(36)
$\displaystyle-(\omega_{1}-\omega_{2})u_{2,-}^{(1)}+(k_{1}-k_{2})\phi_{2,-}^{(1)}=f_{2}^{(-)}\Psi_{1}\Psi_{2}^{\star},$
(37)
$\displaystyle-\left[(k_{1}-k_{2})^{2}+c_{1}\right]\phi_{2,-}^{(2)}+n_{2,-}^{(1)}==f_{3}^{(-)}\Psi_{1}\Psi_{2}^{\star},$
(38)

where $f_{1}^{(\pm)}$, $f_{2}^{(\pm)}$, and $f_{3}^{(\pm)}$ 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(\theta_{1}+\theta_{2})}$ and
$e^{i(\theta_{2}-\theta_{2})}$, respectively. Their solution reads:

$\displaystyle
n_{2,+}^{(1)}=C_{n,2,+}^{(1)}\Psi_{1}\Psi_{2},~{}~{}~{}~{}n_{2,-}^{(1)}=C_{n,2,-}^{(1)}\Psi_{1}\Psi_{2}^{\star},$
(39) $\displaystyle
u_{2,+}^{(1)}=C_{u,2,+}^{(1)}\Psi_{1}\Psi_{2},~{}~{}~{}~{}u_{2,-}^{(1)}=C_{u,2,-}^{(1)}\Psi_{1}\Psi_{2}^{\star},$
(40)
$\displaystyle\phi_{2,+}^{(1)}=C_{\phi,2,+}^{(1)}\Psi_{1}\Psi_{2},~{}~{}~{}~{}\phi_{2,-}^{(1)}=C_{\phi,2,-}^{(1)}\Psi_{1}\Psi_{2}^{\star},$
(41)

where again the coefficients $C_{n,2,\pm}^{(1)}$, $C_{u,2,\pm}^{(1)}$, and
$C_{\phi,2,\pm}^{(1)}$ are given in Appendix A.

In the same (second) order ($\sim\varepsilon^{2}$), we also obtain an equation
involving the “constant” terms (zeroth order harmonics) in the form

$-c_{1}\phi_{2}^{(0)}+n_{2}^{(0)}=+2c_{2}\left(|\Psi_{1}|^{2}+|\Psi_{2}|^{2}\right),$
(42)

which will be used in the next subsection.

Before moving on to the equations which are third order in $\varepsilon$, we
return to the discussion of the solution of Eqs. (30). The latter, which
actually consists of two $3\times 3$ uncoupled systems for the variables
$\\{n_{2,1}^{(1)},u_{2,1}^{(1)},\phi_{2,1}^{(1)}\\}$ and
$\\{n_{2,2}^{(1)},u_{2,2}^{(1)},\phi_{2,2}^{(1)}\\}$, can be written in matrix
form as

$\displaystyle\begin{bmatrix}-\omega_{j}&+k_{j}&0\\\ 0&-\omega_{j}&+k_{j}\\\
1&0&-(k_{j}^{2}+c_{1})\end{bmatrix}\begin{bmatrix}n_{2,j}^{(1)}\\\
u_{2,j}^{(1)}\\\ \phi_{2,j}^{(1)}\end{bmatrix}=\begin{bmatrix}\mu_{1,j}\\\
\mu_{2,j}\\\ \mu_{3,j}\end{bmatrix},$ (43)

where $\mu_{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:

$\displaystyle D_{j}^{\prime}=\begin{vmatrix}\mu_{1,j}&+k_{j}&0\\\
\mu_{2,j}&-\omega_{1}&+k_{j}\\\
\mu_{3,j}&0&-(k_{j}^{2}+c_{1})\end{vmatrix}=0.$ (44)

These determinants result from the replacement of the first column of the
determinant $D_{j}$ with the vector
${\bf\mu}_{j}=[\mu_{1,j}\,\mu_{2,j}\,\mu_{3,j}]^{T}$. After calculations we
obtain that in order for $D_{j}^{\prime}$ to vanish, the following
compatibility conditions

$\frac{\partial\Psi_{j}}{\partial
t_{1}}=-v_{g,j}\frac{\partial\Psi_{j}}{\partial x_{1}}$ (45)

(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].

###  2.3 Derivation of a system of coupled NLS equations

The equations in order $\varepsilon^{3}$ are

$\displaystyle\frac{\partial n_{3}}{\partial
t_{0}}+\frac{\partial(u_{3})}{\partial x_{0}}={\cal G}_{1},$
$\displaystyle\frac{\partial u_{3}}{\partial
t_{0}}+\frac{\partial\phi_{3}}{\partial x_{0}}={\cal G}_{2},$ (46)
$\displaystyle\frac{\partial^{2}\phi_{3}}{\partial
x_{0}^{2}}-c_{1}\phi_{3}+n_{3}={\cal G}_{3},$

where ${\cal G}_{1}$, ${\cal G}_{2}$, and ${\cal G}_{3}$ are expressions
involving lower ($\leq 2$) order harmonics of $n$, $u$, and $\phi$, which are
given in Appendix B. We shall seek a multi-harmonic solution by substituting
into the system of Eqs. (46) the following Ansatz

$\displaystyle
S_{3}=S_{3}^{(0)}+\left[]S_{3,1}^{(1)}e^{i\theta_{1}}+S_{3,1}^{(2)}e^{2i\theta_{1}}+S_{3,1}^{(3)}e^{3i\theta_{1}}+S_{3,2}^{(1)}e^{i\theta_{2}}+S_{3,2}^{(2)}e^{2i\theta_{2}}+S_{3,2}^{(3)}e^{3i\theta_{2}}\right.$
$\displaystyle\left.+S_{3,+}^{(1)}e^{i(\theta_{1}+\theta_{2})}+S_{3,-}^{(1)}e^{2i(\theta_{1}-\theta_{2})}+S_{3,+12}^{(1)}e^{i(\theta_{1}+2\theta_{2})}+S_{3,-12}^{(1)}e^{i(\theta_{1}-2\theta_{2})}\right.$
$\displaystyle\left.+S_{3,+21}^{(1)}e^{i(2\theta_{1}+\theta_{2})}+S_{3,-21}^{(1)}e^{i(2\theta_{1}-\theta_{2})}+c.c.\right],$
(47)

where $S=n,u$, and $\phi$ and $c.c.$ denotes the complex conjugate of the
expression within the curly brackets.

Isolating the $\ell=0$ contributions, we obtain two equations for the zeroth
harmonic amplitudes which can be combined with Eq. (42) to provide a $3\times
3$ system to be satisfied by the (3) unknowns $n_{2}^{(0)}$, $u_{2}^{(0)}$,
and $\phi_{2}^{(0)}$, i.e.

$\displaystyle\frac{\partial n_{2}^{(0)}}{\partial t_{1}}+\frac{\partial
u_{2}^{(0)}}{\partial x_{1}}$ $\displaystyle=$
$\displaystyle-2\frac{\partial}{\partial
x_{1}}\left[\left(\frac{k_{1}}{\omega_{1}}\right)^{3}|\Psi_{1}|^{2}+\left(\frac{k_{2}}{\omega_{2}}\right)^{3}|\Psi_{2}|^{2}\right],$
(48) $\displaystyle\frac{\partial u_{2}^{(0)}}{\partial
t_{1}}+\frac{\partial\phi_{2}^{(0)}}{\partial x_{1}}$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial
x_{1}}\left[\left(\frac{k_{1}}{\omega_{1}}\right)^{2}|\Psi_{1}|^{2}+\left(\frac{k_{2}}{\omega_{2}}\right)^{2}|\Psi_{2}|^{2}\right]$
(49) $\displaystyle n_{2}^{(0)}$ $\displaystyle=$ $\displaystyle
c_{1}\phi_{2}^{(0)}+2c_{2}\left(|\Psi_{1}|^{2}+|\Psi_{2}|^{2}\right)\,.$ (50)

In writing down the latter equations, we have also used the relations (27). To
solve the $3\times 3$ system thus obtained for the zeroth harmonics, we shall
introduce the Ansatz

$\displaystyle\phi_{2}^{(0)}=C_{\phi,2,1}^{(0)}|\Psi_{1}|^{2}+C_{\phi,2,2}^{(0)}|\Psi_{2}|^{2},$
$\displaystyle
u_{2}^{(0)}=C_{u,2,1}^{(0)}|\Psi_{1}|^{2}+C_{u,2,2}^{(0)}|\Psi_{2}|^{2},$ (51)
$\displaystyle
n_{2}^{(0)}=C_{n,2,1}^{(0)}|\Psi_{1}|^{2}+C_{n,2,2}^{(0)}|\Psi_{2}|^{2},$

into Eqs. (48)-(50) and use Eq. (15) to determine the coefficients
$C_{\phi,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 $\varepsilon$) for the
(3rd-rder) first harmonic mode $\ell=1$ read

$\displaystyle-i\omega_{j}n_{3,j}^{(1)}+ik_{j}u_{3,j}^{(1)}={\cal R}_{1,j},$
$\displaystyle-i\omega_{j}u_{3,j}^{(1)}+ik_{j}\phi_{3,j}^{(1)}={\cal
R}_{2,j},$ (52)
$\displaystyle-\left(k_{j}^{2}+c_{1}\right)\phi_{3,j}^{(1)}+n_{3,j}^{(1)}={\cal
R}_{3,j},$

where the functions ${\cal R}_{1,j}$, ${\cal R}_{2,j}$, and ${\cal R}_{3,j}$
are given in Appendix B.

A lengthy algebraic manipulation of Eqs. (52) – aimed at eliminating the
variables $n_{3,j}^{(1)}$, $u_{3,j}^{(1)}$, and $\phi_{3,j}^{(1)}$ – leads us
to an equation in the form:

$-i\omega_{j}{\cal R}_{1,j}-ik_{j}{\cal R}_{2,j}+\omega_{j}^{2}{\cal
R}_{2,j}=0,$ (53)

where the expressions for ${\cal R}_{m,j}$ ($m=1,2,3$ and $j=1,2$) are given
by Eqs. (B)-(B) in Appendix B.

Finally, by substituting the expressions for ${\cal R}_{m,j}$ into Eq. (53)
and using the obtained solutions for the harmonic amplitudes $n_{1,j}^{(1)}$,
$u_{1,j}^{(1)}$, $\phi_{1,j}^{(1)}=\Psi_{j}$, $n_{2,j}^{(2)}$,
$u_{2,j}^{(2)}$, $\phi_{2,j}^{(2)}$, $n_{2}^{(0)}$, $u_{2}^{(0)}$, and
$\phi_{2}^{(0)}$, we obtain after a tedious and lengthy calculation with
equations

$\displaystyle 4\omega_{1}k_{1}\frac{\partial^{2}\Psi_{1}}{\partial
t_{1}\partial
x_{1}}-\frac{k_{1}^{2}}{\omega_{1}^{2}}\frac{\partial^{2}\Psi_{1}}{\partial
t_{1}^{2}}+(1-\omega_{1}^{2})\frac{\partial^{2}\Psi_{1}}{\partial
x_{1}^{2}}+2i\frac{k_{1}^{2}}{\omega_{1}}\frac{\partial^{2}\Psi_{1}}{\partial
t_{2}^{2}}+2ik_{1}(1-\omega_{1}^{2})\frac{\partial\Psi_{1}}{\partial x_{2}}$
$\displaystyle+\left(\tilde{Q}_{11}|\Psi_{1}|^{2}+\tilde{Q}_{12}|\Psi_{2}|^{2}\right)\Psi_{1}=0,$
(54) $\displaystyle 4\omega_{2}k_{2}\frac{\partial^{2}\Psi_{2}}{\partial
t_{1}\partial
x_{1}}-\frac{k_{2}^{2}}{\omega_{2}^{2}}\frac{\partial^{2}\Psi_{2}}{\partial
t_{1}^{2}}+(1-\omega_{2}^{2})\frac{\partial^{2}\Psi_{2}}{\partial
x_{1}^{2}}+2i\frac{k_{2}^{2}}{\omega_{2}}\frac{\partial^{2}\Psi_{2}}{\partial
t_{2}^{2}}+2ik_{2}(1-\omega_{2}^{2})\frac{\partial\Psi_{2}}{\partial x_{2}}$
$\displaystyle+\left(\tilde{Q}_{21}|\Psi_{1}|^{2}+\tilde{Q}_{22}|\Psi_{2}|^{2}\right)\Psi_{2}=0,$
(55)

where

$\tilde{Q}_{jj}=-2\frac{k_{j}^{3}}{\omega_{j}}\left(C_{u,2,j}^{(2)}+C_{u,2,j}^{(0)}\right)-k_{j}^{2}\left(C_{n,2,j}^{(2)}+C_{n,2,j}^{(0)}\right)+2c_{2}\omega_{j}^{2}\left(C_{\phi,2,j}^{(2)}+C_{\phi,2,j}^{(0)}\right)+3c_{3}\omega_{j}^{2},$
(56)

and

$\displaystyle\tilde{Q}_{12}=-2\frac{k_{1}^{3}}{\omega_{1}}C_{u,2,2}^{(0)}-k_{1}\frac{k_{2}}{\omega_{2}}\left(\omega_{1}\frac{k_{2}}{\omega_{2}}+k_{1}\right)\left(C_{u,2,+}^{(1)}+C_{u,2,-}^{(1)}\right)$
$\displaystyle-
k_{1}^{2}C_{n,2,2}^{(0)}-\omega_{1}k_{1}\frac{k_{2}}{\omega_{2}}\left(C_{n,2,+}^{(1)}+C_{n,2,-}^{(1)}\right)$
$\displaystyle+2c_{2}\omega_{1}^{2}\left(C_{\phi,2,2}^{(0)}+C_{\phi,2,+}^{(1)}+C_{\phi,2,-}^{(1)}\right)+6c_{3}\omega_{1}^{2},$
(57)
$\displaystyle\tilde{Q}_{21}=-2\frac{k_{2}^{3}}{\omega_{2}}C_{u,2,1}^{(0)}-k_{2}\frac{k_{1}}{\omega_{1}}\left(\omega_{2}\frac{k_{1}}{\omega_{1}}+k_{2}\right)$
$\displaystyle\left(C_{u,2,+}^{(1)}+C_{u,2,-}^{(1)}\right)-k_{2}^{2}C_{n,2,1}^{(0)}-\omega_{2}k_{2}\frac{k_{1}}{\omega_{1}}\left(C_{n,2,+}^{(1)}+C_{n,2,-}^{(1)}\right)$
$\displaystyle+2c_{2}\omega_{2}^{2}\left(C_{\phi,2,1}^{(0)}+C_{\phi,2,+}^{(1)}+C_{\phi,2,-}^{(1)}\right)+6c_{3}\omega_{2}^{2}.$
(58)

Note that the terms proportional to the variables $n_{2,j}^{(1)}$,
$u_{2,j}^{(1)}$, and $\phi_{2,j}^{(1)}$ were eliminated with the help of Eqs.
(30).

From Eqs. (2.3) and (2.3), it is straightforward to obtain a pair of coupled
NLS equations in the familiar form

$\displaystyle i\left(\frac{\partial\Psi_{1}}{\partial
t_{2}}+v_{g,1}\frac{\partial\Psi_{1}}{\partial
x_{2}}\right)+P_{1}\frac{\partial^{2}\Psi_{1}}{\partial
x_{1}^{2}}+\left(Q_{11}|\Psi_{1}|^{2}+Q_{12}|\Psi_{2}|^{2}\right)\Psi_{1}=0,$
(59) $\displaystyle i\left(\frac{\partial\Psi_{2}}{\partial
t_{2}}+v_{g,2}\frac{\partial\Psi_{2}}{\partial
x_{2}}\right)+P_{2}\frac{\partial^{2}\Psi_{2}}{\partial
x_{1}^{2}}+\left(Q_{21}|\Psi_{1}|^{2}+Q_{22}|\Psi_{2}|^{2}\right)\Psi_{2}=0\,,$
(60)

where (for $j=1,2$)

$P_{j}=-\frac{3}{2}c_{1}\frac{k_{j}}{(k_{j}^{2}+c_{1})^{5/2}},$ (61)

are the (2) dispersion coefficients (analogous to the group-velocity-
dispersion/GVD terms in nonlinear optics; note that
$P_{j}=\frac{1}{2}\frac{\partial\omega_{j}}{\partial k_{j}}$, and

$Q_{jj}=\frac{\omega_{j}}{2k_{j}^{2}}\tilde{Q}_{jj},~{}~{}~{}~{}Q_{12}=\frac{\omega_{1}}{2k_{1}^{2}}\tilde{Q}_{12},~{}~{}~{}~{}Q_{21}=\frac{\omega_{2}}{2k_{2}^{2}}\tilde{Q}_{21}$
(62)

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.)

Figure 2:  The frequency dispersion coefficient $P_{1}$, and the nonlinearity
coefficient $Q_{11}$, calculated from Eqs. (61) and the first of Eqs. (62), as
a function of the wavenumber $k_{1}$ for $\kappa=2$ (black), $\kappa=3$ (red),
$\kappa=4$ (green), $\kappa=5$ (blue), $\kappa=6$ (orange), $\kappa=100$
(brown).

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 $\kappa$. In the Maxwell-Boltzmann case (i.e.,
corresponding practically to the value of $\kappa=100$), the minimum value of
$P_{1}$ is located at $k_{1}\simeq 0.5$, while it shifts to high values of
$k_{1}$ with increasing $\kappa$.

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 $\kappa$ 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}\equiv k_{r}=1.47$, in agreement with
earlier results [8, 9, 61].

#### Critical carrier wavenumbers for non-Maxwellian plasmas.

Generally, for arbitrary $\kappa$, 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};\kappa)=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 $\kappa$.

Figure 3:  (a) The root $k_{r}$ of $Q_{11}$ as a function of the spectral
index $\kappa$. For large values of $\kappa$, that root converges
asymptotically to the value $k_{1}=k_{r}=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 $\kappa$. The orange-dashed line is located
at $k_{r}=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
$\kappa=2$ (black), $\kappa=3$ (red), $\kappa=4$ (green), $\kappa=5$ (blue),
$\kappa=6$ (orange), $\kappa=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).

The dependence of the root $k_{r}$ on the spectral index $\kappa$ is shown in
Fig. 3(a) for a large interval of $\kappa$ values. As it can be observed in
Fig. 3(a), the root $k_{r}$ diverges as $\kappa\rightarrow 3/2$. The value of
the root $k_{r}$ decreases rapidly as $\kappa$ increases from the value $3/2$,
then reaches a minimum around $\kappa\simeq 2.5$ close to the border of far-
to near-equilibrium, and then, with further increasing $\kappa$, it increases
slowly to converge to $k_{r}=1.47$ from below for very large values of
$\kappa$ (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 $\kappa$ in Fig. 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}$) that might be associated with negative $Q_{12}$ for
Maxwellian electrons (i.e. large $\kappa$) now switch to a positive $Q_{12}$
in the case of suprathermal electrons, i.e. for smaller $\kappa$. 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 $\kappa$. (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 $\kappa=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 $\kappa$, i.e., for
$\kappa=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].

Figure 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) $\kappa=2$,
(b) $\kappa=3$, (c) $\kappa=4$, (d) $\kappa=5$, (e) $\kappa=6$, (f)
$\kappa=100$. In (f), the $\kappa$ 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}$.).

###  2.4 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 $\Psi_{1}$ and $\Psi_{2}$. Let us
introduce the new independent variables $\xi$ and $\tau$ through $\xi=x-vt$
and $\tau=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

$\displaystyle
i\left(\frac{\partial\Psi_{1}}{\partial\tau}+\delta\frac{\partial\Psi_{1}}{\partial\xi}\right)+P_{1}\frac{\partial^{2}\Psi_{1}}{\partial\xi^{2}}+\left\\{Q_{11}|\Psi_{1}|^{2}+Q_{12}|\Psi_{2}|^{2}\right\\}\Psi_{1}=0,$
(63) $\displaystyle
i\left(\frac{\partial\Psi_{2}}{\partial\tau}-\delta\frac{\partial\Psi_{2}}{\partial\xi}\right)+P_{2}\frac{\partial^{2}\Psi_{2}}{\partial\xi^{2}}+\left\\{Q_{21}|\Psi_{1}|^{2}+Q_{22}|\Psi_{2}|^{2}\right\\}\Psi_{2}=0,$
(64)

where $\delta=(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 $\xi_{1}$, while the higher-order differentiation in the
beginning of the left-hand-side are with respect to $\xi_{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
$\Psi_{j}=\bar{\Psi}_{j}\exp\left[i\left(\frac{\delta^{2}}{4P_{j}}\tau-\frac{\delta}{2P_{j}}\xi\right)\right]$,
where $j=1,2$, we obtain

$\displaystyle
i\frac{\partial\bar{\Psi}_{1}}{\partial\tau}+P_{1}\frac{\partial^{2}\bar{\Psi}_{1}}{\partial\xi^{2}}+\left\\{Q_{11}|\bar{\Psi}_{1}|^{2}+Q_{12}|\bar{\Psi}_{2}|^{2}\right\\}\bar{\Psi}_{1}=0,$
(65) $\displaystyle
i\frac{\partial\bar{\Psi}_{2}}{\partial\tau}+P_{2}\frac{\partial^{2}\bar{\Psi}_{2}}{\partial\xi^{2}}+\left\\{Q_{21}|\bar{\Psi}_{1}|^{2}+Q_{22}|\bar{\Psi}_{2}|^{2}\right\\}\bar{\Psi}_{2}=0,$
(66)

where $\bar{\Psi}_{j}$ are complex functions of the new variables $\xi$ and
$\tau$.

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.

##  3 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

$\Psi_{j}=\Psi_{j,0}\,e^{i\tilde{\omega}_{j}\tau},$ (67)

where $\Psi_{j,0}$ ($j=1,2$) is a constant amplitude (here assumed to be real)
and $\tilde{\omega}_{j}$ is a frequency, to be determined. By inserting Eq.
(67) into Eqs. (63) and (64), we obtain

$\displaystyle\tilde{\omega}_{1}=Q_{11}\Psi_{1,0}^{2}+Q_{12}\Psi_{2,0}^{2},$
(68)
$\displaystyle\tilde{\omega}_{2}=Q_{21}\Psi_{1,0}^{2}+Q_{22}\Psi_{2,0}^{2}.$
(69)

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
$\varepsilon_{j}$ to their amplitudes ($|\varepsilon_{j}|\ll|\Psi_{j,0}|$), as

$\Psi_{j}=(\Psi_{j,0}+\varepsilon_{j})\,e^{i\tilde{\omega}_{j}\tau}.$ (70)

By inserting the perturbed $\Psi_{j}$ into Eqs. (63) and (64) we obtain

$\displaystyle
i\left(\frac{\partial\varepsilon_{1}}{\partial\tau}+\delta\frac{\partial\varepsilon_{1}}{\partial\xi}\right)+P_{1}\frac{\partial^{2}\varepsilon_{1}}{\partial\xi^{2}}+Q_{11}\Psi_{1,0}^{2}(\varepsilon_{1}+\varepsilon_{1}^{\star})+Q_{12}{\Psi}_{1,0}\Psi_{2,0}(\varepsilon_{2}+\varepsilon_{2}^{\star})=0,$
(71) $\displaystyle
i\left(\frac{\partial\varepsilon_{2}}{\partial\tau}-\delta\frac{\partial\varepsilon_{2}}{\partial\xi}\right)+P_{2}\frac{\partial^{2}\varepsilon_{2}}{\partial\xi^{2}}+Q_{22}\Psi_{2,0}^{2}(\varepsilon_{2}+\varepsilon_{2}^{\star})+Q_{21}\Psi_{2,0}\Psi_{1,0}(\varepsilon_{1}+\varepsilon_{1}^{\star})=0.$
(72)

Then, by setting $\varepsilon_{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

$\displaystyle\left\\{\left(\frac{\partial}{\partial\tau}+\delta\frac{\partial}{\partial\xi}\right)^{2}+P_{1}\left(P_{1}\frac{\partial^{2}}{\partial\xi^{2}}+2Q_{11}\Psi_{1,0}^{2}\right)\frac{\partial^{2}}{\partial\xi^{2}}\right\\}g_{1}$
$\displaystyle+2P_{1}Q_{12}\Psi_{2,0}\Psi_{1,0}\frac{\partial^{2}}{\partial\xi^{2}}g_{2}=0\,,$
(73)
$\displaystyle\left\\{\left(\frac{\partial}{\partial\tau}-\delta\frac{\partial}{\partial\xi}\right)^{2}+P_{2}\left(P_{2}\frac{\partial^{2}}{\partial\xi^{2}}+2Q_{22}\Psi_{2,0}^{2}\right)\frac{\partial^{2}}{\partial\xi^{2}}\right\\}g_{2}$
$\displaystyle+2P_{2}Q_{21}\Psi_{1,0}\Psi_{2,0}\frac{\partial^{2}}{\partial\xi^{2}}g_{1}=0\,.$
(74)

Eventually, by setting

$g_{j}=g_{j,0}\,e^{i(K\xi-\Omega\tau)}+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

$\left[(\Omega-\delta K)^{2}-\Omega_{1}^{2}\right]\left[(\Omega+\delta
K)^{2}-\Omega_{2}^{2}\right]=\Omega_{c}^{4},$ (76)

where ($j=1,2$)

$\displaystyle\Omega_{j}^{2}$ $\displaystyle=$ $\displaystyle
P_{j}K^{2}\left(P_{j}K^{2}-2Q_{jj}\Psi_{j,0}^{2}\right),$ (77)
$\displaystyle\Omega_{c}^{4}$ $\displaystyle=$ $\displaystyle
4P_{1}P_{2}Q_{12}Q_{21}\Psi_{1,0}^{2}\Psi_{2,0}^{2}K^{4}.$ (78)

Figure 5:  Maps of the growth rate $\Gamma$ 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$, $\Psi_{1,0}=\Psi_{2,0}=0.15$, and (a) $\kappa=2$; (b) $\kappa=3$;
(c) $\kappa=4$; (d) $\kappa=5$; (e) $\kappa=6$; (f) $\kappa=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}<k_{r}$ 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}>k_{r}$ where $P_{1}Q_{11}<0$ and
$P_{2}Q_{22}>0$). Note that $k_{r}$ depends on the particular value of the
spectral index $\kappa$ in accordance with the inset of Fig. 3(a).

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 $\Gamma$. 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, $\Gamma$, i.e.,

$\Gamma=max[Im(\Omega_{r})],$ (79)

where $\Omega_{r}$ denotes one (each) of the four roots of the compatibility
condition.

Figure 6:  Maps of the growth rate $\Gamma$ on the $k_{2}-K$ plane for (a)
$\kappa=1.75$; (b) $\kappa=2.00$; (c) $\kappa=2.25$; (d) $\kappa=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}<k_{r}$, 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}>k_{r}$, right from the
vertical line). Note that $k_{r}$ depends weakly on the spectral index
$\kappa$, but the differences cannot be observed in this scale.

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 (3) and (3) 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 $\Omega$ as a function of the perturbation wavenumber
$K$ are given by

$\Omega=+\delta K\pm\Omega_{1},~{}~{}~{}~{}\Omega=-\delta K\pm\Omega_{2},$
(80)

for the first and the second NLS equation that results from decoupling Eqs.
(3) and (3), 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}(\kappa)$ while $Q_{jj}<0$ for $k_{j}>k_{r}(\kappa)$, 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}(\kappa)$, 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\rightarrow 2$, i.e. it will be stable for
$k_{2}<k_{r}(\kappa)$ 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 $PQ<0$ and instability occurs otherwise.

##  4 Modulational instability: numerical parametric analysis

We have calculated the growth rate $\Gamma$ numerically for the CNLS equations
by finding the roots of the polynomial in $\Omega$ 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 $\kappa$, i.e., for $\kappa=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 $\Gamma$ is mapped on the $k_{2}-K$
plane for $k_{1}=0.4<k_{r}(\kappa)$ and the six values of the spectral index
$\kappa$ 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}(\kappa)$
($P_{2}Q_{22}<0$) while it is modulationally unstable for
$k_{2}>k_{r}(\kappa)$ ($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 $\Gamma$ 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.

Figure 7:  Maps of the growth rate $\Gamma$ 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$, $\Psi_{1,0}=\Psi_{2,0}=0.15$, and (a) $\kappa=2$; (b) $\kappa=3$;
(c) $\kappa=4$; (d) $\kappa=5$; (e) $\kappa=6$; (f) $\kappa=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}<k_{r}$ 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}>k_{r}$ where $P_{1}Q_{11}<0$ and
$P_{2}Q_{22}>0$). Note that $k_{r}$ depends on the particular value of the
spectral index $\kappa$ in accordance with the inset of Fig. 3(a).  Figure 8:
Maps of the growth rate $\Gamma$ on the $k_{2}-K$ plane for (a) $\kappa=1.75$;
(b) $\kappa=2.00$; (c) $\kappa=2.25$; (d) $\kappa=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}<k_{r}$, 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}>k_{r}$, right from the vertical line). The
stability situation in these maps is the same as that in figs. 5 and 6. Note
that $k_{r}$ depends weakly on the spectral index $\kappa$, but the
differences cannot be observed in this scale.

Generally speaking, the decrease of the spectral index $\kappa$ increases the
areas on the $k_{2}-K$ plane in which MI occurs. Moreover, the maximum values
of the growth rate $\Gamma$ seem to increase considerably with decreasing
$\kappa$ (note the different scales in the colorbars). Note also the
difference in the mudulational instability patterns for values of $\kappa$ in
the far- and near-equilibrium region, i.e., for $\kappa=2$ and
$\kappa=3,4,5,6,100$, respectively, which are separated by a rough empirical
boundary at around $\kappa\simeq 2.5$. Indeed, for $\kappa=2$, the growth rate
$\Gamma$ 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 $\kappa$, i.e., for $\kappa=3$,
the CNLS system is modulationally stable in the most part of the area in which
the decoupled CNLS equations are in the stable-stable regime. With further
increasing the value of $\kappa$, 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 $\Gamma$ 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 $\Gamma$ seems to occur for values of
the spectral index $\kappa$ within the strongly nonthermal region, i.e., for
$1.5<\kappa<2.5$. To the contrary, the variation of the $\Gamma$ pattern for
$\kappa=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 $\Gamma$ in these areas are concerned. In Fig. 6, maps of $\Gamma$ on the
$k_{2}-K$ plane are shown for more values of $\kappa$ within the strongly
nonthermal region, i.e., for $\kappa=1.75$, $2.25$, and $2.50$. In the figure,
the map of the already shown map for $\kappa=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 $\Gamma$ patterns with the
variation of $\kappa$. For $\kappa=1.75$ (Fig. 6(a)), in particular,
modulational instability accompanied by very high values of $\Gamma$ 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 $\kappa$ up to $1.90$, where
where some modulationally stable islands begin to appear on the $k_{2}-K$
plane. When $\kappa$ 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 $\Gamma$
occur. With further increasing $\kappa$, 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 $\kappa=3$).

Figure 9:  Maps of the growth rate $\Gamma$ 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$, $\Psi_{1,0}=\Psi_{2,0}=0.15$, and (a) $\kappa=2$; (b) $\kappa=3$;
(c) $\kappa=4$; (d) $\kappa=5$; (e) $\kappa=6$; (f) $\kappa=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}>k_{r}$ 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}<k_{r}$ where $P_{1}Q_{11}>0$ and
$P_{2}Q_{22}<0$). Note that $k_{r}$ depends on the particular value of the
spectral index $\kappa$ in accordance with the inset of Fig. 3(a).

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}(\kappa)$, while they are in the stable-unstable regime for
$k_{2}>k_{r}(\kappa)$ ($P_{1}Q_{11}<0$ and $P_{2}Q_{22}>0$). In comparison
with Fig. 5, we observe that for $\kappa=3,4,5,6,100$, areas of strong
instability with relatively high values of the growth rate $\Gamma$ 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 $\kappa=3$ to $100$. A similar situation can be described for
$\kappa=2$, where high values of the growth rate $\Gamma$ occur at low $k_{2}$
as in the case of higher $\kappa$; however, the values of $\Gamma$ 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 $\Gamma$ 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 $\kappa\geq 3$.

As above, we also present maps of $\Gamma$ for a few more values of $\kappa$
in the strongly nonthermal region, i.e., in the interval $1.5<\kappa<2.5$ in
Fig. 8. Roughly, the remarks made for the patterns of $\Gamma$ in Fig. 6 also
hold here. For $\kappa=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 $\kappa$, the modulationally stable areas on the $k_{2}-K$ plane
grow to the expense of the modulationally unstable ones. The pattern of
$\Gamma$ 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
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}(\kappa)$.
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 $\kappa$ continuous for $\kappa>2.5$. Note that
in this case the strongly modulationally unstable area for low $k_{2}$
persist, although becoming much smaller with increasing $\kappa$, even for
very large values of $\kappa$, i.e., $\kappa=100$ (Fig. 7(f)).

Figure 10:  Maps of the growth rate $\Gamma$ on the $k_{2}-K$ plane for (a)
$\kappa=1.75$; (b) $\kappa=2.00$; (c) $\kappa=2.25$; (d) $\kappa=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}<k_{r}$, 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}>k_{r}$, right from the
vertical line). Note that $k_{r}$ depends weakly on the spectral index
$\kappa$, but the differences cannot be observed in this scale.

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}(\kappa)$. 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}(\kappa)$, 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}(\kappa)$, 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 $\kappa\geq 3$. For all values of
$\kappa$, we again observe strong modulational instability with high values of
$\Gamma$ for low $k_{2}$. For $\kappa=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 unstable-unstable regime.
For $\kappa\geq 3$ we again observe some shrinkage of the instability areas
with increasing $\kappa$. We can also observe for these values of $\kappa$
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 $\Gamma$
on the $k_{2}-K$ plane for a few values of $\kappa$ in the strongly nonthermal
region, i.e. for small values of $\kappa$. 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 $\kappa$ values within the strongly nonthermal region.

Figure 11:  The growth rate $\Gamma$ of the CNLS system (black curves) and the
growth rates $\Gamma_{1}$ and $\Gamma_{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 $\kappa=3$,
$\Psi_{1,0}=\Psi_{2,0}=0.15$, and (a) $k_{1}=1$ and $k_{2}=0.6$, for which
$P_{1}Q_{11}<0$ and $P_{2}Q_{22}<0$ (stable-stable); (b) $k_{1}=1.6$ and
$k_{2}=1.35$, for which $P_{1}Q_{11}>0$ and $P_{2}Q_{22}<0$ (unstable-stable);
(c) $k_{1}=1.6$ and $k_{2}=1.45$, for which $P_{1}Q_{11}>0$ and
$P_{2}Q_{22}>0$ (unstable-unstable). Note the difference in the horizontal
scales. The growth rates $\Gamma_{1}$ and $\Gamma_{2}$ are always zero in (a),
while $\Gamma_{2}$ is always zero in (b).

#### Representative examples of growth rate variation.

It may be interesting to present some illustrative examples for the growth
rate $\Gamma$ 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 $\Gamma(K)$ curves are shown in three cases that differ in respect to
the stability regime of the decoupled CNLS equations, for $\kappa=3$. Along
with these curves, we also show the corresponding growth rates $\Gamma_{1}$
and $\Gamma_{2}$ obtained for the first and the second decoupled CNLS
equation, respectively. In Fig. 11(a), in which $k_{1}=1<k_{r}(\kappa)$ and
$k_{2}=0.6<k_{r}(\kappa)$, i.e., in the stable-stable regime of the decoupled
CNLS equations, we obtain a $\Gamma(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 $\Gamma_{1}$ and $\Gamma_{2}$ are both zero. In
Fig. 11(b), in which $k_{1}=1.6>k_{r}(\kappa)$ and $k_{2}=1.35>k_{r}(\kappa)$,
i.e., in the unstable-stable regime of the decoupled CNLS equations, we obtain
a $\Gamma(K)$ curve with two humps. Moreover, we obtain a non-zero
$\Gamma_{1}(\kappa)$ curve which has the commonly observed shape. A comparison
between the two indicates that besides the difference in shape, the values of
$\Gamma_{1}(\kappa)$ are always significantly smaller than those of
$\Gamma(\kappa)$. 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}(\kappa)$ and
$k_{2}=1.45>k_{r}(\kappa)=1.42$, i.e., in the unstable-unstable regime of the
decoupled CNLS equations, we obtain three curves for the growth rates
$\Gamma(K)$, $\Gamma_{1}$, and $\Gamma_{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.

##  5 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
$\kappa$ which characterizes the kappa distribution of the electrons, i.e.,
$k_{r}=k_{r}(\kappa)$. 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 $\kappa-$dependence
of $k_{r}$ is not monotonous but it exhibits a minimum around $\kappa\simeq
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}(\kappa)$
dependence is weak for most values of $\kappa$, while it becomes stronger for
$\kappa<2$ where it increases abruptly until its divergence at $\kappa=3/2$.
For very large values of $\kappa$, 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 $\Gamma$, 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 $\kappa$ 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 $\Gamma$ patterns
on the $k_{2}-K$ plane occurs for values of the spectral index $\kappa$
withing the strongly nonthermal region, i.e., in the interval
$1.5<\kappa<2.5$. Moreover, as the numerical calculation of the growth rate
$\Gamma$ indicate, the CNLS equations are modulationally unstable in almost
all of the $k_{2}-K$ plane shown in the figures below $\kappa\sim 1.9$. For
these values of $\kappa$, i.e., for $\kappa<1.9$, the values of $\Gamma$
increase fast with decreasing $\kappa$ and diverge in the limit of
$\kappa\rightarrow 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 $\Omega_{c}^{4}$ 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” $\Omega_{c}^{4}$ 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}$
such 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}(\kappa)$. 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 $\kappa$, will be
the focus of future work.

##  6 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.

##  7 Declarations

###  7.1 Funding

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) (PI Ioannis Kourakis). IK also acknowledges
financial support from KU via the project FSU-2021-012 (8474000352) (PI
Ioannis Kourakis) as well as from KU Space and Planetary Science Center, via
grant No. KU-SPSC-8474000336 (PI Mohamed Ramy Mohamed Elmaarry).

###  7.2 Competing Interests

The authors have no relevant financial or non-financial interests to disclose.

###  7.3 Author Contributions

I. Kourakis contributed to the study conception, design, and methodology,
while N. Lazarides contributed to the methodology, software development and
numerics. Both authors contributed to the analysis of the results. The first
draft was written by N. Lazarides and the final manuscript by I. Kourakis.
Both authors read and approved the final manuscript.

###  7.4 Data Availability

The datasets generated during and/or analysed during the current study are
available from the corresponding author on reasonable request.

##  Appendix A Coefficients - Analytical Expressions

$\displaystyle
C_{n,2,j}^{(2)}=\frac{(k_{j}^{2}+c_{1})}{6k_{j}^{2}}\left[3(4k_{j}^{2}+c_{1})(k_{j}^{2}+c_{1})-2c_{2}\right],$
(81) $\displaystyle
C_{u,2,j}^{(2)}=\frac{\sqrt{k_{j}^{2}+c_{1}}}{6k_{j}^{2}}\left[3(k_{j}^{2}+c_{1})(2k_{j}^{2}+c_{1})-2c_{2}\right],$
(82) $\displaystyle
C_{\phi,2,j}^{(2)}=\frac{1}{6k_{j}^{2}}\left[3(k_{j}^{2}+c_{1})^{2}-2c_{2}\right],$
(83)

$\displaystyle C_{n,2,\pm}^{(1)}=\frac{1}{D_{\pm}}\left\\{\left[(k_{1}\pm
k_{2})^{2}+c_{1}\right]\left[(k_{1}\pm
k_{2})f_{2}^{(\pm)}+(\omega_{1}\pm\omega_{2})f_{1}^{(\pm)}\right]+(k_{1}\pm
k_{2})^{2}f_{3}^{(\pm)}\right\\},$ $\displaystyle
C_{u,2,\pm}^{(1)}=\frac{1}{D_{\pm}}\left\\{(k_{1}\pm
k_{2})\left[f_{1}^{(\pm)}+(\omega_{1}\pm\omega_{2})f_{3}^{(\pm)}\right]+(\omega_{1}\pm\omega_{2})\left[(k_{1}\pm
k_{2})^{2}+c_{1}\right]f_{2}^{(\pm)}\right\\},$ (84) $\displaystyle
C_{\phi,2,\pm}^{(1)}=\frac{1}{D_{\pm}}\left\\{(\omega_{1}\pm\omega_{2})\left[f_{1}^{(\pm)}+(\omega_{1}\pm\omega_{2})f_{3}^{(\pm)}\right]+(k_{1}\pm
k_{2})f_{2}^{(\pm)}\right\\},$

where

$D_{\pm}=(k_{1}\pm k_{2})^{2}-(\omega_{1}\pm\omega_{2})^{2}\left[(k_{1}\pm
k_{2})^{2}+c_{1}\right].$ (85)

The quantities $f_{m}^{(\pm)}$ ($m=1,2,3$) are defined in Eqs. (33)-(33).

$C_{n,2,j}^{(0)}=\left(c_{1}C_{\phi,2,j}^{(0)}+2c_{2}\right),\qquad
C_{u,2,j}^{(0)}=\frac{1}{v_{g,j}}\left[C_{\phi,2,j}^{(0)}+\frac{k_{j}^{2}}{\omega_{j}^{2}}\right],\qquad
C_{\phi,2,j}^{(0)}=\rho_{j}\frac{1}{1-c_{1}v_{g,j}^{2}},$ (86)

where

$\rho_{j}=-(k_{j}^{2}+c_{1})-2c_{1}+2c_{2}v_{g,j}^{2}.$ (87)

##  Appendix B Right-Hand-Sides of the Perturbative Equations.

$\displaystyle{\cal F}_{1}=-\frac{\partial(n_{1})}{\partial
t_{1}}-\frac{\partial(u_{1})}{\partial
x_{1}}-\frac{\partial(n_{1}u_{1})}{\partial x_{0}},\qquad{\cal
F}_{2}=-\frac{\partial u_{1}}{\partial t_{1}}-\frac{\partial\phi_{1}}{\partial
x_{1}}-u_{1}\frac{\partial u_{1}}{\partial x_{0}},$ $\displaystyle{\cal
F}_{3}=+c_{2}\phi_{1}^{2}-2\frac{\partial^{2}\phi_{1}}{\partial x_{0}\partial
x_{1}}.$ (88)

$\displaystyle\mu_{1,j}=-\left(\frac{k_{j}}{\omega_{j}}\right)^{2}\frac{\partial\Psi_{j}}{\partial
t_{1}}-\frac{k_{j}}{\omega_{j}}\frac{\partial\Psi_{j}}{\partial x_{1}},$ (89)
$\displaystyle\mu_{2,j}=-\frac{k_{j}}{\omega_{j}}\frac{\partial\Psi_{j}}{\partial
t_{1}}-\frac{\partial\Psi_{j}}{\partial x_{1}},$ (90)
$\displaystyle\mu_{3,j}=-2ik_{j}\frac{\partial\Psi_{j}}{\partial x_{1}}.$

$\displaystyle
f_{1}^{(+)}=-(k_{1}+k_{2})\frac{k_{1}k_{2}}{\omega_{1}\omega_{2}}\left(\frac{k_{1}}{\omega_{1}}+\frac{k_{2}}{\omega_{2}}\right),$
(91) $\displaystyle
f_{2}^{(+)}=-(k_{1}+k_{2})\frac{k_{1}k_{2}}{\omega_{1}\omega_{2}},$ (92)
$\displaystyle f_{3}^{(+)}=+2c_{2}.$ (93)

$\displaystyle
f_{1}^{(-)}=-(k_{1}-k_{2})\frac{k_{1}k_{2}}{\omega_{1}\omega_{2}}\left(\frac{k_{1}}{\omega_{1}}+\frac{k_{2}}{\omega_{2}}\right),$
(94) $\displaystyle
f_{2}^{(-)}=-(k_{1}-k_{2})\frac{k_{1}k_{2}}{\omega_{1}\omega_{2}},$ (95)
$\displaystyle f_{3}^{(-)}=+2c_{2}.$ (96)

$\displaystyle{\cal G}_{1}=-\frac{\partial(n_{2})}{\partial
t_{1}}-\frac{\partial(u_{2})}{\partial x_{1}}-\frac{\partial(n_{1})}{\partial
t_{2}}-\frac{\partial(u_{1})}{\partial
x_{2}}-\frac{\partial(n_{1}u_{2}+n_{2}u_{1})}{\partial
x_{0}}-\frac{\partial(n_{1}u_{1})}{\partial x_{1}},$ (97) $\displaystyle{\cal
G}_{2}=-\frac{\partial u_{2}}{\partial t_{1}}-\frac{\partial\phi_{2}}{\partial
x_{1}}-\frac{\partial u_{1}}{\partial t_{2}}-\frac{\partial\phi_{1}}{\partial
x_{2}}-\frac{\partial u_{1}u_{2}}{\partial x_{0}}-u_{1}\frac{\partial
u_{1}}{\partial x_{1}},$ (98) $\displaystyle{\cal
G}_{3}=-2\frac{\partial^{2}\phi_{2}}{\partial x_{0}\partial
x_{1}}-2\frac{\partial^{2}\phi_{1}}{\partial x_{0}\partial
x_{2}}-\frac{\partial^{2}\phi_{1}}{\partial
x_{1}^{2}}+2c_{2}\phi_{1}\phi_{2}+c_{3}\phi_{1}^{3}.$ (99)

$\displaystyle{\cal R}_{1,j}=-\left(\frac{\partial n_{2,j}^{(1)}}{\partial
t_{1}}+\frac{\partial u_{2,j}^{(1)}}{\partial x_{1}}+\frac{\partial
n_{1,j}^{(1)}}{\partial t_{2}}+\frac{\partial u_{1,j}^{(1)}}{\partial
x_{2}}\right)$ $\displaystyle-
ik_{j}\left(n_{1,j}^{(-1)}u_{2,j}^{(2)}+n_{1,j}^{(1)}u_{2}^{(0)}+n_{1,k_{j}}^{(-1)}u_{2,p}^{(1)}+n_{1,k_{j}}^{(1)}u_{2,m}^{(\ell_{j})}\right)$
$\displaystyle-
ik_{j}\left(u_{1,j}^{(-1)}n_{2,j}^{(2)}+u_{1,j}^{(1)}n_{2}^{(0)}+u_{1,k_{j}}^{(-1)}n_{2,p}^{(1)}+u_{1,k_{j}}^{(1)}n_{2,m}^{(\ell_{j})}\right),$
(100) $\displaystyle{\cal R}_{2,j}=-\left(\frac{\partial
u_{2,j}^{(1)}}{\partial t_{1}}+\frac{\partial\phi_{2,j}^{(1)}}{\partial
x_{1}}+\frac{\partial u_{1,j}^{(1)}}{\partial
t_{2}}+\frac{\partial\Psi_{j}}{\partial x_{2}}\right)$ $\displaystyle-
ik_{j}\left(u_{1,j}^{(-1)}u_{2,j}^{(2)}+u_{1,j}^{(1)}u_{2}^{(0)}+u_{1,k_{j}}^{(-1)}u_{2,p}^{(1)}+u_{1,k_{j}}^{(1)}u_{2,m}^{(\ell_{j})}\right),$
(101) $\displaystyle{\cal
R}_{3,j}=-2ik_{j}\left(\frac{\partial\phi_{2,j}^{(1)}}{\partial
x_{1}}+\frac{\partial\Psi_{j}}{\partial x_{2}}\right)$
$\displaystyle-\frac{\partial^{2}\Psi_{j}}{\partial
x_{1}^{2}}+2c_{2}\left(\Psi_{j}^{\star}\phi_{2,j}^{(2)}+\Psi_{j}\phi_{2}^{(0)}+\Psi_{k_{j}}^{\star}\phi_{2,p}^{(1)}+\Psi_{k_{j}}\phi_{2,m}^{(\ell_{j})}\right)$
$\displaystyle+3c_{3}\Psi_{j}\left(|\Psi_{j}|^{2}+2|\Psi_{k_{j}}|^{2}\right),$
(102)

where $j=1,2$, and $k_{j}$, $\ell_{j}$ defined as $\ell_{1}=1$, $\ell_{2}=-1$,
$k_{1}=2$, and $k_{2}=1$.

## References

  * [1] Livadiotis, G., McComas, D. J.: Understanding Kappa Distributions: A Toolbox for Space Science and Astrophysics. Space Sci. Rev. 175, 183–214 (2013). https://doi.org/10.1007/s11214-013-9982-9
  * [2] Nicolaou, G., Livadiotis, G., Owen, C. J. Verscharen, D., Wicks, R. T.: Determining the kappa distributions of space plasmas from observations in a limited energy range. Astrophys. J. 864, 3 (2018). https://doi.org/10.3847/1538-4357/aad45d
  * [3] Nicolaou, G., Livadiotis, G., Wicks, R. T.: On the determination of kappa distribution functions from space plasma observations. Entropy 22, 212 (2020). https://doi.org/10.3390/e22020212
  * [4] Livadiotis, G.: Collision frequency and mean free path for plasmas described by kappa distributions. AIP Adv. 9, 105307 (2019). https://doi.org/10.1063/1.5125714
  * [5] Saberian, E., Livadiotis, G.: Plasma oscillations and spectral index in non-extensive statistics. Physica A 593, 126909 (2022). https://doi.org/10.1016/j.physa.2022.126909
  * [6] Kalita, J., Das, R., Hosseini, K., Baleanu, D., Salahshour, S.: Solitons in magnetized plasma with electron inertia under weakly relativistic effect. Nonlinear Dyn. 111, 3701–3711 (2023). https://doi.org/10.1007/s11071-022-08015-z
  * [7] Madhukalya, B., Das, R., Hosseini, K., Baleanu, D., Hincal, E.: Effect of ion and negative ion temperatures on KdV and mKdV solitons in a multicomponent plasma. Nonlinear Dyn. 111, 8659–8671 (2023). https://doi.org/10.1007/s11071-023-08262-8
  * [8] Shimizu K., Ichikawa, Y. H.: Automodulation of ion oscillation modes in plasma. J. Phys. Soc. Jpn. 33, 789–792 (1972). https://doi.org/10.1143/JPSJ.33.789
  * [9] Kakutani, T., Sugimoto, N.: Krylov-Bogoliubov-Mitropolsky method for nonlinear wave modulation. Phys. Fluids 17, 1617–1625 (1974). https://doi.org/10.1063/1.1694942
  * [10] Kourakis, I., Shukla, P. K.: Exact theory for localized envelope modulated electrostatic wavepackets in space and dusty plasmas. Nonlinear Process. in Geophys. 12, 407–423 (2005). https://doi.org/10.5194/npg-12-407-2005
  * [11] Chowdhury, N. A., Mannan, A., Hossen, M. R., Mamun, A. A.: Modulational instability and generation of envelope solitons in four-component space plasmas. Contrib. Plasma Phys. 58, 870–877 (2018). https://doi.org/10.1002/ctpp.201700069
  * [12] Singh, K., Saini, N. S.: Breather structures and Peregrine solitons in a polarized Space dusty plasma. Front. Phys. 8, 602229 (2020). https://doi.org/10.3389/fphy.2020.602229
  * [13] Sarkar, J., Chandra, S., Goswami, J., Ghosh, B.: Formation of solitary structures and envelope solitons in electron acoustic wave in inner magnetosphere plasma with suprathermal ions. Contrib. Plasma Phys. 60, e201900202 (2020). https://doi.org/10.1002/ctpp.201900202
  * [14] Borhanian, J., Kourakis, I., Sobhanian, S.: Electromagnetic envelope solitons in magnetized plasma. Phys. Lett. A 373, 3667–3677 (2009). https://doi.org/10.1016/j.physleta.2009.08.010
  * [15] Borhanian, J.: Extraordinary electromagnetic localized structures in plasmas: Modulational instability, envelope solitons, and rogue waves. Phys. Lett. A 379, 595–602 (2015). http://dx.doi.org/10.1016/j.physleta.2014.12.018
  * [16] Veldes, G. P., Borhanian, J., McKerr, M., Saxena, V., Frantzeskakis, D. J., Kourakis, I.: Electromagnetic rogue waves in beam-plasma interactions. J. Opt. 15, 064003 (2013). 10.1088/2040-8978/15/6/064003
  * [17] Borhanian, J.: Dissipative ion-acoustic solitary and shock waves in a plasma with superthermal electrons. Plasma Phys. Control. Fusion 55, 105012 (2013). http://iopscience.iop.org/0741-3335/55/10/105012
  * [18] Sultana, S., Kourakis, I.: Electron-scale dissipative electrostatic solitons in multi-species plasmas. Phys. Plasmas 22, 102302 (2015). http://dx.doi.org/10.1063/1.4932071
  * [19] Sultana, S., Schlickeiser, R., Elkamash, I. S., Kourakis, I.: Dissipative high-frequency envelope soliton modes in nonthermal plasmas. Phys. Rev. E 98, 033207 (2018). https://doi.org/10.1103/PhysRevE.98.033207
  * [20] Singh, K., Saini, N. S.: The evolution of rogue wave triplets and super rogue waves in superthermal polarized space dusty plasma. Phys. Plasmas 26, 113702 (2019). https://doi.org/10.1063/1.5119894
  * [21] Arham, M., Khan, S. A., Khan, M.: Weak dissipation of electrostatic solitary structures in warm collisional pair-ion plasmas with non-Maxwellian electron population. Chin. J. Phys. 69, 77–88 (2021). https://doi.org/10.1016/j.cjph.2020.10.028
  * [22] Sharmin, B. E., Shikha, R. K., Tamanna, N. K., Chowdhury, N. A., Mannan, A., Mamun, A. A.: Modulational instability of dust-ion-acoustic waves and associated first and second-order rogue waves in a super-thermal plasma. Results Phys. 26, 104373 (2021). https://doi.org/10.1016/j.rinp.2021.104373
  * [23] Yahia, M. E., Tolba, R. E., Moslem, W. M.: Super rogue wave catalysis in Titanas’ ionosphere. Adv. Space Res. 67, 1412–1424 (2021). https://doi.org/10.1016/j.asr.2020.11.027
  * [24] Xie, Y., Li, L., Zhu, S.: Dynamical behaviors of blowup solutions in trapped quantum gases: Concentration phenomenon. J. Math. Anal. Appl. 468, 169–181 (2018). https://doi.org/10.1016/j.jmaa.2018.08.011
  * [25] Li, L., Xie, Y., Zhu, S.: New exact solutions for a generalized KdV equation. Nonlinear Dyn. 92, 215–219 (2018). https://doi.org/10.1007/s11071-018-4050-3
  * [26] Xie, Y., Li, Y., Kang, Y.: New solitons and conditional stability to the high dispersive nonlinear Schrödinger equation with parabolic law nonlinearity. Nonlinear Dyn. 103, 1011–1021 (2021). https://doi.org/10.1007/s11071-020-06141-0
  * [27] Li, L., Yan, Y., Xie, Y.: Localized excitation and folded solitary wave for an extended (3+1)-dimensional B-type Kadomtsev-Petviashvili equation. Nonlinear Dyn. 109, 2013–2027 (2022). https://doi.org/10.1007/s11071-022-07559-4
  * [28] Li, L., Nie, Y., Zhu, M., Xie, Y.: Variable separation solution and multi-valued soliton of an extended (3+1)-dimensional B-type Kadomtsev–Petviashvili equation. Nonlinear Dyn. 111, 4723–4736 (2023). https://doi.org/10.1007/s11071-022-08092-0
  * [29] Spatschek, K. H.: Coupled localized electron-plasma waves and oscillatory ion-acoustic perturbations. Phys. Fluids 21, 1032–1035 (1978). https://doi.org/10.1063/1.862323
  * [30] Som, B. K., Gupta, M. R., Dasgupta, B.: Coupled nonlinear Schröinger equation for Langmuir and dispersive ion-acoustic waves. Phys. Lett. 72A, 111–114 (1979). https://doi.org/10.1016/0375-9601(79)90663-7
  * [31] Tabi, C. B., Panguetna, C. S., Motsumi, T. G., Kofané, T. C.: Modulational instability of coupled waves in electronegative plasmas. Phys. Scripta 95, 075211 (2020). https://doi.org/10.1088/1402-4896/ab8f40
  * [32] McKinstrie, C. J., Bingham, R.: The modulational instability of coupled waves. Physics of Fluids B: Plasma Physics 1, 230 (1989). https://doi.org/10.1063/1.859095
  * [33] McKinstrie, C. J., Luther, G. G.: The modulational instability of colinear waves. Phys. Scripta 30, 31–40 (1990). 10.1088/0031-8949/1990/T30/005
  * [34] Luther, G. G., McKinstrie, C. J,: Transverse modulational instability of collinear waves. J. Opt. Soc. Am. B 7, 1125–1141 (1990). https://doi.org/10.1364/JOSAB.7.001125
  * [35] Luther, G. G., McKinstrie, C. J.: Transverse modulational instability of counterpropagating light waves. J. Opt. Soc. Am. B 9, 1047–1060 (1992). https://doi.org/10.1364/JOSAB.9.001047
  * [36] Kourakis, I., Shukla, P. K., Morfill, G.: Modulational instability and localized excitations involving two nonlinearly coupled upper-hybrid waves in plasmas. New J. Phys. 7, 153 (2005). https://dx.doi.org/10.1088/1367-2630/7/1/153
  * [37] Singh, V.: Modulation instability of two laser beams in plasma. Laser Part. Beams 31, 753–758 (2013). https://doi.org/10.1017/S0263034613000748
  * [38] Borhanian, J., Golijan, H. A.: Copropagation of coupled laser pulses in magnetized plasmas: Modulational instability and coupled solitons. Phys. Plasmas 24, 033116 (2017). https://doi.org/10.1063/1.4978576
  * [39] Guo, H.-D., Xia, T.-C.: Multi-soliton solutions for a higher-order coupled nonlinear Schrödinger system in an optical fiber via Riemann-Hilbert approach. Nonlinear Dyn. 103, 1805–1816 (2021). https://doi.org/10.1007/s11071-020-06166-5
  * [40] Zhou, Q., Xu, M., Sun, Y., Zhong, Y., Mirzazadeh, M.: Generation and transformation of dark solitons, anti-dark solitons and dark double-hump solitons. Nonlinear Dyn. 110, 1747–1752 (2022). https://doi.org/10.1007/s11071-022-07673-3
  * [41] Xiang, X. S., Zuo, D. W.: Breather and rogue wave solutions of coupled derivative nonlinear Schrödinger equations. Nonlinear Dyn. 107, 1195–1204 (2022). https://doi.org/10.1007/s11071-021-07050-6
  * [42] Jin, J., Zhang, Y.,Ye, R., Wu, L.: The breather and semi-rational rogue wave solutions for the coupled mixed derivative nonlinear Schrödinger equations. Nonlinear Dyn. 111, 633–643 (2023). https://doi.org/10.1007/s11071-022-07834-4
  * [43] Wu, X. H., Gao, Y. T., Yu, X., Li, L.-Q, Ding, C. C.: Vector breathers, rogue and breather-rogue waves for a coupled mixed derivative nonlinear Schrödinger system in an optical fiber. Nonlinear Dyn. 111, 5641–5653 (2023). https://doi.org/10.1007/s11071-022-08058-2
  * [44] Zhang, Y., Yang, C., Yu, W., Mirzazadeh, M., Zhou, Q., Liu, W.: Interactions of vector anti-dark solitons for the coupled nonlinear Schrödinger equation in inhomogeneous fibers. Nonlinear Dyn. 94, 1351–1360 (2018). https://doi.org/10.1007/s11071-018-4428-2
  * [45] Yu, W., Liu, W., Triki, H., Zhou, Q., Biswas, A., Belić, M. R.: Control of dark and anti-dark solitons in the (2+1)-dimensional coupled nonlinear Schrödinger equations with perturbed dispersion and nonlinearity in a nonlinear optical system. Nonlinear Dyn. 97, 471–483 (2019). https://doi.org/10.1007/s11071-019-04992-w
  * [46] Patel, A., Kumar, V.: Modulation instability analysis of a nonautonomous (3+1)-dimensional coupled nonlinear Schrödinger equation. Nonlinear Dyn. 104, 4355–4365 (2021). https://doi.org/10.1007/s11071-021-06558-1
  * [47] Yang, D. Y., Tian, B., Wang, M., Zhao, X., Shan, W. R., Jiang, Y.: Lax pair, Darboux transformation, breathers and rogue waves of an N-coupled nonautonomous nonlinear Schrödinger system for an optical fiber or a plasma. Nonlinear Dyn. 107, 2657–2666 (2022). https://doi.org/10.1007/s11071-021-06886-2
  * [48] Wang, L., Luan, Z., Zhou, Q., Biswas, A., Alzahrani, A. K., Liu, W.: Bright soliton solutions of the (2+1)-dimensional generalized coupled nonlinear Schröinger equation with the four-wave mixing term. Nonlinear Dyn. 104, 2613–2620 (2021). https://doi.org/10.1007/s11071-021-06411-5
  * [49] Yan, X. W., Zhang, J.: Coupled cubic-quintic nonlinear Schrödinger equation: novel bright–dark rogue waves and dynamics. Nonlinear Dyn. 100, 3733–3743 (2020). https://doi.org/10.1007/s11071-020-05694-4
  * [50] Ren, P., Rao, J.: Bright-dark solitons in the space-shifted nonlocal coupled nonlinear Schrödinger equation. Nonlinear Dyn. 108, 2461–2470 (2022). https://doi.org/10.1007/s11071-022-07269-x
  * [51] Wu, H. Y., Jiang, L. H.: Diverse excitations of two-component rogue waves for a nonautonomous coupled partially nonlocal nonlinear Schröinger model under a parabolic potential. Nonlinear Dyn. 109, 1993–2002 (2022). https://doi.org/10.1007/s11071-022-07510-7
  * [52] Meng, G.-Q., Qin, J.-L., Yu, G.-L.: Breather and rogue wave solutions for a nonlinear Schrödinger-type system in plasmas. Nonlinear Dyn. 81, 739–751 (2015). https://doi.org/10.1007/s11071-015-2024-2
  * [53] Onorato, M., Osborne, A. R., Serio, M.: Modulational Instability in Crossing Sea States: A Possible Mechanism for the Formation of Freak Waves. Phys. Rev. Lett. 96, 014503 (2006). https://doi.org/10.1103/PhysRevLett.96.014503
  * [54] Haelterman, M., Sheppard, A.: Bifurcation phenomena and multiple soliton-bound states in isotropic Kerr media. Phys. Rev. E 49 (4), 3376–3381 (1994). https://doi.org/10.1103/PhysRevE.49.3376
  * [55] Seve, E., Millot, G., Wabnitz, S., Sylvestre, T., Maillotte, H.: Generation of vector dark-soliton trains by induced modulational instability in a highly birefringent fiber. J. Opt. Soc. Am. B 16 (10), 1642–1650 (1999). https://doi.org/10.1364/JOSAB.16.001642
  * [56] Frisquet, B., Kibler, B., Morin, M., Baronio, F., Conforti, M., Millot, G., Wabnitz, S.: Optical Dark Rogue Wave. Sci. Rep. 6, 20785 (2016). https://doi.org/10.1038/srep20785
  * [57] Lazarides, N., Tsironis, G. P.: Coupled nonlinear Schrödinger field equations for electromagnetic wave propagation in nonlinear left-handed materials. Phys. Rev. E 71, 036614 (2005). https://doi.org/10.1103/PhysRevE.71.036614
  * [58] Stalin, S., Ramakrishnan, R., Senthilvelan, M., Lakshmanan, M.: Nondegenerate Solitons in Manakov System. Phys. Rev. Lett. 122, 043901 (2019). https://doi.org/10.1103/PhysRevLett.122.043901
  * [59] Yilmaz, E. U., Khodad, F. S., Ozkan, Y. S., Abazari, R., Abouelregal, A. E., Shaayesteh, M. T., Rezazadeh, H., Ahmad, H.: Manakov model of coupled NLS equation and its optical soliton solutions. Journal of Ocean Engineering and Science (2022). https://doi.org/10.1016/j.joes.2022.03.005
  * [60] Kourakis, I., Shukla, P. K.: Modulational instability in asymmetric coupled wave functions. Eur. Phys. J. B 50, 321–325 (2006). https://doi.org/10.1140/epjb/e2006-00106-1
  * [61] Lazarides, N., Veldes, G. P., Javed, A., Kourakis, I.: Modulational electrostatic wave-wave interactions in plasma fluids modeled by non-integrable coupled nonlinear Schrödinger (CNLS) equations. Chaos, Solitons and Fractals 175, 113974 (2023). https://doi.org/10.1016/j.chaos.2023.113974
  * [62] Pierrard, V., Lazar, M.: Kappa distributions: Theory and applications in Space plasmas. Solar Phys. 267, 153–174 (2010). https://doi.org/10.1007/s11207-010-9640-2
  * [63] Livadiotis, G., Kappa distributions: Theory and applications in plasmas. Amsterdam, Netherlands; Elsevier, 2017. isbn: 9780128046388.
  * [64] Baluku, T. K., Hellberg, M. A., Kourakis, I., Saini, N. S.: Dust ion acoustic solitons in a plasma with kappa-distributed electrons. Phys. Plasmas 17, 053702 (2010). https://doi.org/10.1063/1.3400229
  * [65] Hellberg, A., Mace, R. L., Baluku, T. K., Kourakis, I., Saini, N. S.: Comment on mathematical and physical aspects of Kappa velocity distributions. Phys. Plasmas 16, 094701/1-5 (2009). http://dx.doi.org/10.1063/1.3213388
  * [66] Livadiotis, G., McComas, T. J.: Beyond kappa distributions: exploiting Tsallis statistical mechanics in space. J. Geophys. Res. 114, 11105 (2009). 10.1029/2009JA014352
  * [67] Kourakis, I., Sultana, S., Hellberg, M. A.: Plasma Phys. Control. Fusion 54, 124001 (2012). 10.1088/0741-3335/54/12/124001
  * [68] Infeld, E., Rowlands, G.: Nonlinear Waves, Solitons and Chaos. Cambridge University Press, 1990. isbn: 0521379377.
  * [69] He, Y., Slunyaev, A., Mori, N., Chabchoub, A.: Experimental evidence of nonlinear focusing in standing water waves. Phys. Rev. Lett. 129, 144502 (2022). https://doi.org/10.1103/PhysRevLett.129.144502
  * [70] Ablowitz, M. J., Horikis, T. P.: Interacting nonlinear wave envelopes and rogue wave formation in deep water. Phys. Fluids 27, 012107 (2015). https://doi.org/10.1063/1.4906770
  * [71] Veldes, G. P., Cuevas, J., Kevrekidis, P. G., Frantzeskakis, D. J.: Coupled backward- and forward-propagating solitons in a composite right- and left-handed transmission line. Phys. Rev. E 88, 013203 (2013). https://doi.org/10.1103/PhysRevE.88.013203
  * [72] Kivshar, Y. S., Turitsyn, S. K.: Vector dark solitons. Opt. Lett. 18, 337–339 (1993). https://doi.org/10.1364/OL.18.000337
  * [73] Tyutin, V. V.: Extended vector solitons with significantly different frequencies of the polarization components. JETP Letters 115, 634–637 (2022). https://doi.org/10.1134/S0021364022600690
  * [74] Huang, W.-T., Liu, F.-F., Lü, X., Wang, J. P., Xu, H.-T.: Optical soliton and modulation instability in the high birefringence fiber. Nonlinear Dyn. 108, 2429–2445 (2022). https://doi.org/10.1007/s11071-022-07212-0
  * [75] Sulem, C., Sulem, P.-L.: The nonlinear Schrödinger equation: Self-focusing and wave collapse. Applied Mathematical Sciences, 139, Springer, 1999. isbn: 0387986111.