On the Existence of Subsonic Solitary Waves Associated with Reconnection Jets in Earth’s Magnetotail

The Magnetospheric Multiscale Spacecraft (MMS) has detected the signature of electrostatic solitary waves (ESWs) occurring in the reconnection jet site of the Earth’s magnetotail (Liu et al.). These observations have motivated us to explore the mechanism underlying the formation of fast- and slow-mode ion-acoustic solitary waves in the magnetotail region. To this end, we have formulated a three-component magnetized plasma model consisting of nonthermal electrons and two cold ion beams streaming parallel and antiparallel to the magnetic field, respectively. In this work, we have examined the existence conditions for ion-acoustic subsonic waves in a suprathermal space plasma comprising two counterstreaming (drifting) ion beams interacting with a suprathermal electron background. An exact (nonperturbative) nonlinear technique has been adopted to examine the role of the beam velocity and the spectral index on the evolution of subsonic waves. Linear analysis reveals that subsonic waves are unstable when the beam velocity is lower than a threshold value; hence in this regime, only conventional supersonic (fast) waves are formed. On the other hand, when the beam velocity exceeds the threshold, either supersonic or subsonic waves may exist. The combined impact of the beam velocity and electron superthermality on the characteristics of subsonic solitary waves has been analyzed. Our results are shown to be in good agreement with observations of slow ESWs by the MMS spacecraft. Our findings will help to unfold the so-far unexplored dynamical characteristics of subsonic waves that may occur in the reconnection site of Earth’s magnetotail.


Introduction
Electrostatic solitary waves (ESWs) are a common occurrence in Space plasma observations.In their most common realization, these are identified in measurements as bipolar electric field waveforms, associated with pulse-shaped electrostatic potential excitations and localized electron/ion density disturbances that (co-)propagate along the ambient magnetic field lines, mainly.From a modeling perspective, the main analytical tools for the study of ESWs were elaborated by Sagdeev (1966) and Washimi & Taniuti (1966).Over the last years, various theoretical and experimental investigations have focused on solitary waves occurring in different plasma environments (Baboolal et al. 1990;Berthomier et al. 1998;Hellberg & Verheest 2008;Lakhina et al. 2011;Mahmood & Ur-Rehman 2013;Ur-Rehman et al. 2014;Varghese & Ghosh 2020;Varghese et al. 2022).ESWs have been commonly observed in the Earth's magnetosphere by spacecraft missions, e.g., Viking (Boström et al. 1988), Geotail (Matsumoto et al. 1994), the Five-hundred-meter Aperture Spherical Telescope (Ergun et al. 1998), Cluster (Pickett et al. 2004).In addition to near-Earth space plasma environments, planetary missions such as Cassini and MAVEN have also observed similar electrostatic structures (ESWs) in the magnetospheres of Saturn (Pickett et al. 2015) and Mars (Andersson et al. 2015;Kakad et al. 2022), respectively.These observations have inspired plasma physicists to examine the formation of ESWs in various space and planetary environments by theoretical models, followed by simulation studies of solitary waves in multicomponent plasmas (Kakad et al. 2014(Kakad et al. , 2016b;;Lotekar et al. 2016;Singh et al. 2020Singh et al. , 2021Singh et al. , 2022)).
The existence of fast and slow ion-acoustic modes has been established in multi-warm ion plasma models based on linear wave theory, followed by fascinating nonlinear analyses of plasmas permeated by two counterstreaming ion beams (Lakhina et al. 2020;Verheest & Hellberg 2021).Recent plasma studies, adopting a plasma-fluid formalism, have revealed that beam-permeated plasmas may not only support the conventional supersonic (superacoustic) electrostatic structures (solitary waves, SWs) but also a subsonic (subacoustic) ion-acoustic nonlinear mode propagating at a speed below the sound speed but above a certain Mach number threshold (M min ; Lakhina et al. 2020;Verheest & Hellberg 2021).Interestingly, Papadopoulos et al. (1971) studied the heating of counterpropagating ion beams propagating across an ambient field, based on linear theory in combination with computer simulations.It was shown that the ion-ion beam instability is saturated above a certain threshold (i.e., when the beam velocity exceeds the ion sound speed).In a review article, Bret (2009) has laid out a unified description of various types of instability (filamentation, two-stream, Buneman, Bell, etc.) by considering a full three-dimensional dielectric tensor for a cold relativistic electron beam in a cold plasma, thus accounting for a guiding magnetic field: that description is valid for any orientations of the wavevector, which incorporates a variety of unstable modes.
An interesting new feature is that the SW amplitude may increase as the Mach number reduces toward M min , below Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
which no solitary wave may exist.Recently, Lakhina et al. (2021) also studied the generation mechanism of such waves in reconnection jets by taking into account beam thermal effects (although these were not discussed by Liu et al. 2019), assuming ion beams with Maxwellian electrons in the background.
A plasma-fluid model consisting of two hot counterstreaming ion beams and Maxwellian electrons was adopted by Lakhina et al. (2021Lakhina et al. ( , 2020)), who were the first to show that SWs may exist below critical Mach number, provided that certain conditions are met.Verheest & Hellberg (2021) subsequently revisited the same model by neglecting thermal ion effects and pinpointed the crucial role played by the beam speed(s)-rather than thermal effects-in subsonic solitary wave formation.In a meticulous way, they confirmed that such slow-mode SWs are indeed subsonic, as observed in the laboratory frame, in both symmetric and asymmetric beam plasma models.
Different satellite missions have established the predominance of energetic particles in various space and astrophysical environments, wherein the electron velocity distribution features a long-tailed trend in large arguments, accounting for a significant suprathermal component, thus diverging from the widely used "textbook" scenario of a thermal (Maxwell-Boltzmann) distribution (Livadiotis 2017(Livadiotis , 2018)).Suprathermal particles have been observed inter alia in the Earth's magnetospheric (Feldman et al. 1975) and auroral region(s) (Lazar et al. 2008;Mendis & Rosenberg 1994), and also in the magnetosphere of Mercury, as revealed by MESSENGER data (Ho et al. 2011).A kappa-type (non-Maxwellian) distribution was introduced for the first time as a heuristic formula by Vasyliunas (Vasyliunas 1968), to model OGO 1 and OGO 3 spacecraft data recorded in the Earth's magnetosphere.Suprathermal distributions were subsequently adopted to model particle distributions in the solar wind (Armstrong et al. 1983) and in planetary magnetospheres (e.g., Earth's, Saturn's, and Jupiter's; Leubner 1982).
The magnetic reconnection mechanism transforms magnetic energy into plasma kinetic energy, to be further accompanied by alteration in the magnetic field configuration (Fu et al. 2017;Yamada et al. 2010).In the Earth's magnetotail, the situation arises when two oppositely directed confined plasmas are coupled with magnetic field lines, and then these reconnected plasmas are ejected at high velocities from the reconnection site.This whole scenario is referred to as reconnection jets (Cao et al. 2013).It is now established that reconnection jets have a vital role in plasma energization in space and in astrophysical plasmas, e.g., in solar flares, pulsar winds, and active galactic nuclei (Chen et al. 2019;Kirk & Skjaeraasen 2003;Lakhina et al. 2021;Liu et al. 2019).Reconnection jets in the Earth's magnetotail can exhibit different kinds of plasma waves and instabilities, in addition to electrostatic solitary waves.Recently, Liu et al. (2019) reported the first observational evidence, by the Magnetospheric Multiscale Spacecraft (MMS) spacecraft, that unfolded the evolution of ESWs associated with cold ion beams within reconnection jets of Earth's magnetotail.The magnetotail is abundantly permeated by counterstreaming ion beams with relative velocity up to 2000 km s −1 , along with hot electrons, and these can act as a free energy source for the formation of ESWs, as discussed by Liu et al. (2019).In the last decade, slow ionacoustic electron solitary waves were reported in the plasma sheet boundary layer (Norgren et al. 2015) and at the reconnection site (Graham et al. 2015).Subsequently, Kakad et al. (2016a) adopted fluid simulations to examine the formation of slow ion-acoustic ESWs in the plasma sheet boundary layer.Recently, Lakhina et al. (2021) investigated the formation of fast and slow ESWs in reconnection jets by considering the two ion beams with thermal effects and Maxwellian electrons.The observations of Liu et al. (2019) have confirmed the presence of cold ion beams but give no information about the particles' distribution function.
High-speed jets are usually produced by the magnetic reconnection process and are one of the chief sources of suprathermal electrons in space & astrophysical plasmas (Vaivads et al. 2021).Many observations and computational studies reveal that suprathermal electrons are the primary driver in the reconnection process that allows electron acceleration in the reconnection as well as in the outburst regions (Birn et al. 2012;Vaivads et al. 2021).
The work presented in the paper at hand was motivated by the above considerations and focused for the first time on subsonic ESWs occurring in reconnection jet regions in Earth's magnetotail plasma, where subsonic ESWs were observed by satellites.Admittedly, very little can be found in the existing literature along this challenging line of research.Liu et al. (2019) did not report strong currents nor hot ions (with a temperature of ∼10 keV) during the occurrence of ESWs.Accordingly, in our study, we have considered two counterstreaming (cold) ion beams in the presence of suprathermal electrons.
In this article, we shall rely on a (Sagdeev-type) nonlinear pseudopotential technique, in order to explore the existence regime of subsonic waves in a plasma with suprathermal electrons.
This investigation is, to the best of our knowledge, the first of its kind, in that it has focused on a realistic Space plasma situation where two counterstreaming plasmas collide, against a background of nonthermal (kappa-distributed) electrons; this is a ubiquitous scenario in the Earth's magnetosphere and may arguably also occur in other planetary magnetospheres, affected by the solar wind, as a constant source of streaming ions in addition to energized electrons.Earlier studies have focused exclusively on supersonic solitary waves; as a matter of fact, subsonic solutions are ruled out in a quiescent plasma-i.e., in the absence of a beam-despite certain studies claiming the opposite in the past, due to a misinterpretation of the sound speed concept, as discussed, e.g., by Dubinov (2009).However, slow (subacoustic) solitary waves were reportedly detected in satellite observations, and a theoretical framework for these was so far lacking in the literature.Our work was based on the surprising, and arguably pioneering, recent prediction furnished by Lakhina et al. (2021Lakhina et al. ( , 2020)), and later refined by Verheest & Hellberg (2021), that subsonic electrostatic soliton may, in fact, exist in a beam-permeated plasma.
Our results presented in this study should assist in the interpretation of observational data from the MMS mission in reconnection jet areas in Earth's magnetotail.

A Double Beam Plasma Model
Let us consider a plasma consisting of non-Maxwellian electrons (mass m e , charge −e) and two (H + ) ion-beam fluids (mass m 1 = m 2 = m i , charge q 1 = q 2 = + e. (We have considered a charge state Z i = + 1, assuming we are dealing with proton beams, i.e., hydrogen nuclei; this does not affect the generality of our results below, which may apply to any other type of ions, upon a trivial change of scale.)Note that the two beams consist of positive ions of the same kind; however, these are treated as different fluids, in that they differ in equilibrium fluid speed (beam velocity).
We consider electrostatic waves propagating in a direction parallel to the ambient magnetic field (which allows us to neglect the Larmor force in the fluid equations of motion, for simplicity).The fluid model equations, describe the plasma state in terms of the two ion fluid densities (n j ) and the fluid speeds (u j )-where j = 1 or 2 respectively, for each of the two ion fluids-consist of the continuity equation(s): the momentum equation(s): and Poisson's equation:


where f is the electrostatic potential and ò 0 is the permittivity of the vacuum.(Note that the tilde, here used to denote physicali.e., dimensional-quantities, will be dropped later, after rescaling has been applied to derive a dimensionless system of evolution equations).
The electron density entering Poisson's equation can be obtained by integrating the kappa velocity distribution (Hellberg et al. 2009) as At equilibrium, the charge neutrality condition is where n j for ( j = 1 or 2) denotes the unperturbed number density of the ions, respectively.For analytical simplicity in manipulating the above algebraic equations, we shall now apply the following normalization: , (where


); finally, the drifting (streaming) speed of the two ion fluids (i.e., the beams) is also normalized as = U U C j j I A .Note that the scales adopted for the fluid speed, time, and space represent the sound speed, the plasma period, and the Debye length for ionacoustic wave propagation in a "textbook" e-i plasma configuration (i.e., for a single ion component).Let us emphasize, for clarity, that these scales are not to be mistaken for the analogous quantities in our plasma configuration, i.e., C IA is not the (true) sound speed in our case, as will be discussed below, and so on.
Upon applying the above normalization, Equations (1)-( 4) become The charge neutrality requirement (at equilibrium) imposes the constraint: as the equilibrium values of the density variables may differ, in principle, i.e., n 1,0 ≠ n 2,0 .

Linear Analysis
Before investigating the nonlinear analysis of ion-acoustic (IA) waves, we will examine the linear stability of the system briefly.By linearizing Equations ( 6)-( 9) and assuming harmonic oscillations of angular frequency ω and wavenumber k, we obtain the linear dispersion relation This is a quartic equation (i.e., polynomial of 4th degree) with four roots, among which two are associated with fast ionacoustic mode and two are slow modes.As a limiting case, if we assume for a minute that there are no streaming velocities, i.e., U 1 = U 2 = 0, the above relation becomes a quadratic equation in the form: (recall that δ 1 + δ 2 = 1).Remember that the quantity a 1 is actually related to the (kappa-dependent) charge screening length λ Debye and the sound speed C s in non-Maxwellian e-i plasmas, as l ~C a 1 s Debye 1 , as discussed by Kourakis et al. (2012) and independently by Livadiotis & McComas (2014).(a 1 tends to unity in the Maxwellian limit κ e → ∞.) Equation (10) possesses two fast IA modes (and two slow modes), in fact propagating in opposite directions, i.e., one propagates toward the right and the other one toward the left.
Inspired by the analysis in Lakhina et al. (2020) and Verheest & Hellberg (2021), we proceed by assuming a symmetric bi-ion plasma, i.e., δ 1 = δ 2 = δ = 1/2, with U 1 = − U, and U 2 = U.Note that this choice ensures a vanishing total charge current (as expressed in the laboratory frame), in account of the electrostatic approximation (precluding magnetization, via Ampere's law).Substituting into the dispersion relation, we are led to and the corresponding solutions are The subscripts F and S denote the fast and slow IA modes corresponding to the ± sign, respectively.

Nonlinear Analysis
The possibility for the existence of solitons in the plasma of our focus can be studied by employing the well-known (so-called Sagdeev-type) pseudopotential technique (Sagdeev 1966; also see in Verheest & Hellberg 2009 for an exhaustive review of the methodology in the context of Space plasma).Equations ( 7)-( 10) can be expressed into a stationary frame of reference that is moving at speed V (representing the solitary wave speed), i.e., ξ = (x − Mt) where M = V/C IA is the normalized value of the pulse speed.(We have to point out, for rigor, that the term "Mach number" has been adopted here, in agreement with earlier works, even though this is clearly not the true Mach number, as C IA is not the true sound speed.) Now, integrating the continuity and momentum Equations ( 6) and (7) to obtain the ion densities, and taking into account vanishing boundary conditions (for localized solutions) viz.as |ξ| → ± ∞ , we obtain: The ion density functions become infinite for limiting potentials This (twofold) infinite compression limit will enable us to determine the existence domain for soliton pulses, as this is a limit not to be exceeded, in order for all state variables to be real: see that the right-hand side of Equation ( 14) and/or the right-hand side of Equation (15) become imaginary, if f exceeds f l,1 and/or f l,2 .
Substituting Equations ( 14) and (15) in Equation (8) along with Equation (9) and then integrating with respect to f, we find the energy balance equation: where the pseudopotential function S is given by For localized modes (solitary waves) to exist, the Sagdeev pseudopotential S(f, M) must satisfy the following conditions: 0 (where the root f 0 represents the soliton amplitude), (iii) S(f, M) < 0 for 0 < |f| < |f 0 |.In particular, the second derivative yields This will give the acoustic speeds M s for the given plasma system.Note that the latter equation could have been derived from Equation (10) by dividing both sides by k 2 , setting ω/k = v ph , and then taking the limit k = a 1 .It is thus obvious that the lower bound for the soliton speed-satisfying (20), that iscoincides with ( ) w  k lim k 0 , i.e., essentially with the true sound speed in the given plasma configuration.

Symmetric Model
Let us now consider beams of equal strength and absolute drift speed, by taking which can be rewritten as Two roots for M 2 can be computed from the above biquadratic equation, in the form physically representing two values for the acoustic (sound) speed.(Recall that a 1 = (κ e − 1/2)/(κ e − 3/2).)Among the (four) roots thus obtained for M s , we shall retain the (two) positive roots (the remaining two negative solutions stand for wave propagation to the left and will be neglected).Indeed, we need not discuss the case M < 0, as it just follows the same behavior as M > 0 in the symmetric case.Note that the subscript "s" has been used in the latter two expressions, to remind us that these relations essentially provide us with the sound speed of the fast and slow modes, respectively.In the symmetric case, the limiting values Equations ( 16)-( 17) can also be rewritten as For the reality of the state variables, one needs to ensure that both inequalities S(f l,1 ) 0 and S(f l,2 ) 0 hold.As (M + U) 2 > (M − U) 2 , and thus f l,2 < f l,1 , only f l,2 plays a role in limiting the range of values for the soliton speed M, viz.S(f l,2 ) 0 is sufficient a condition to be imposed for reality to be guaranteed.The Sagdeev pseudopotential in the symmetric case becomes which for f = f l,2 can be simplified as An analogous relation appears in (Verheest & Hellberg 2021).This tells us that, for a given value of U, the faster soliton regime, lying between the sound speed M s,1 and the upper limit M l,1 , satisfies U < M s,1 < M < M l,1 .Obviously, for M > M l,1 no solitons can be found.However, there is also a group of slower solitons between M s,2 and M l,2 , that satisfies U > M s,2 > M > M l,2 > 0; these slow solitons cannot propagate with speed M < M l,2 .Such solutions will therefore travel at subsonic speed, a possibility that did not exist in the absence of the beam.Note that there will be no solitary wave between M s,1 and M s,2 .This gap indicates a stopband in the Mach number domain.

Asymmetric Model
As a matter of fact, subsonic solitons owe their very existence to the beam-and not to the symmetric composition of the bi-ion mixture as above, as discussed by Verheest & Hellberg (2021).Inspired by the latter study, we shall consider a specific composition with δ 1 = 2/3, δ 2 = 1/3, U 1 = − U, and U 2 = 2U.This ensures that the equilibrium is charge-and current-neutral, i.e., δ 1 + δ 2 = 1 and δ 1 U 1 + δ 2 U 2 = 0. (As a matter of fact, any combination of values satisfying ∀α ä (0, 1) would work just as well, here, to proceed with.We have chosen α = 2/3, for comparison with the latter reference.) The reality condition now becomes which is, again, a biquadratic equation in M. Depending on U, there are either two or four real roots.Contrary to the fast beam mode, which is always present-i.e., even with U = 0, the slow mode only exists after a certain threshold value of the beam speed U, below which it is unstable.This is also true in the symmetric model.This fact is illustrated in Figures 2, 4, and 10, to be discussed below.
The limiting potential values that will be used in formulating a reality now read: In the asymmetric case, both of these conditions are important to evaluate the upper Mach limit from the requirements S (f l1, M) = 0 and S(f l2, M) = 0. Therefore, for the lower branch, in which the solutions are mostly negative, the limit is given by Equation (30), whereas positive roots are limited by (31).

Parametric Analysis
Numerical analysis will now be used to illustrate the above algebraic results, regarding the existence conditions of slow mode subsonic solitons in non-Maxwellian plasma.

Linear Analysis
We have considered the symmetric model in dispersion relation Equation ( 13) and obtained four roots, consisting of two fast and two slow linear modes.The two panels in Figure 1 depict the variation in the angular frequency ω F (for the fast IA mode) in the k-U plane, corresponding to a) κ e = 20 (quasi-Maxwellian) and b) κ e = 2 (i.e., a strongly nonthermal case).The fast mode monotonically increases with the drifting velocity U.Both the frequency and the phase speed of this (fast IA) mode are lower for a non-Maxwellian (kappa-) distribution (low κ e )-see Figure 1(b))-than in the thermal case (shown in Figure 1(a)).Remember that only the positive solution ω > 0 is taken into account, as ω < 0 will simply be a mirror image, i.e., a symmetric curve on the negative side.
The imaginary part of the angular frequency (i.e., the linear instability growth rate) Im(ω S ) for the slow (unstable) IA mode is shown in Figure 2 (a,b), as it varies on the k−U plane.The same (two) cases are considered, as in the respective panels in Figure 1.The growth rate is notably lower for low κ e (non-Maxwellian case), and the instability window is also narrower on the streaming velocity (U) axis.(Again, only the positive root for ω S was considered.) Figure 3 depicts the real part of the angular frequency Re(ω S ) in the k−U plane, for the slow (stable) IA mode (the same two κ e values have been considered, as in the previous figures).This (slow linear) mode exists (only) above a certain threshold (value) of the beam speed U, and in fact monotonically increases with U. Perhaps counterintuitively, the phase speed of the slow IA mode is higher for lower κ e (nonthermal) than in the Maxwellian case (for equal k and U, that is), in the region where it does exist.respectively, depict the (fast mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of κ e .We have taken U = 1.5, M = 2.3 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).

Symmetric Model
The four panels in Figure 4 illustrate the existence domain of solitary waves associated with the fast and slow IA modes, versus the beam speed (U), for four (4) different-descendingvalues of κ e , assuming a symmetric bi-ion mixture.The fast mode is depicted in blue color, while the slow mode is shown in red.The solid curves correspond to the sonic speed limit M s and the dashed curves correspond to the infinite compression limit M l , while the black dotted curve (for M = U) separates the fast and slow regime(s).The existence domain for the faster IA solitary waves (given in the blue curves) is constrained between M s,1 (lower limit) and M l,1 (upper limit), i.e., U < M s,1 < M < M l,1 .No soliton will exist for M > M l,1 .On the other hand, the existence domain of the slow-mode IA solitary waves (expressed by the red curves) will be confined between M s,2 and M l,2 , viz.U > M s,2 > M > M l,2 > 0. No such soliton will exist for M < M l,2 .Interestingly, the existence region for both fast and slow IA solitary waves actually shrinks, for lower κ e (nonthermal case) in comparison with the quasi-Maxwellian (quasi-thermal) case shown in panel 4(a).(See that the latter panel practically matches the earlier result obtained by Verheest & Hellberg 2021; it also agrees well with Lakhina et al. 2020, if thermal effects are ignored.) Figure 5 shows the variation in the Sagdeev pseudopotential profile for different values of the spectral index (κ e ), still for a beam-free plasma (U = 0).Note that an e-i quiescent plasma is essentially obtained in this case, so this model essentially recovered the standard original case studied by Sagdeev and coworkers (Sagdeev 1966;Verheest & Hellberg 2009).Note that, within the plasma model adopted in our study, only positive polarity solitary structures are predicted.Larger values of the root-in the potential f axis-occur for smaller values of κ e , due to the presence of highly energetic electrons.This trend is reflected in the corresponding potential pulse and (bipolar) or the E-field profiles, respectively shown in panels (b) and (c) of the same Figure.
Figure 6 shows the variation in the Sagdeev pseudopotential profile associated with fast-mode IA solitary waves, for different values of the Mach number, while the beam speed remains fixed (i.e., U = 1.5 here).In this case, both the fast and slow modes exist.For the fast mode, the amplitude of fast IA solitary waves is enhanced as the Mach number increases.A similar trend is also portrayed in the corresponding potential and electric field profiles.
Figure 7 depicts the variation in the Sagdeev pseudopotential profile for different values of the spectral index (κ e ) with fixed streaming speed (U = 1.5, here).Higher amplitudes of IA solitary waves occur for lower κ e , i.e., for strongly non-Maxwellian electrons, due to the presence of energetic suprathermal electron component.The same trend is witnessed in the corresponding potential and electric field profiles.respectively, depict the (slow mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of the Mach number.We have taken U = 1.5, κ e = 3 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).
Figure 8 illustrates the variation in the Sagdeev pseudopotential profile associated with slow IA solitary waves, for different values of the Mach number (for a fixed value of the beam speed, i.e., U = 1.5 in this plot).Rather counterintuitively, for the slow mode, the potential pulse amplitude (i.e., the root of the S − function) of IA solitary waves reduces as the Mach number increases.This surprising trend is also reflected in the corresponding potential and electric field profiles, obtained numerically: see panels (b) and (c) in the same Figure.
Figure 9 portrays the variation in the Sagdeev pseudopotential profile associated with slow-mode IA solitary waves, for different values of the spectral index (κ e ), keeping the beam speed fixed at U = 1.5.Larger values of the potential pulse amplitude of slow-mode IA solitary waves are obtained for lower κ e , as in the case of the fast mode discussed above.This fact is also confirmed by obtaining the corresponding potential and electric field profiles, shown in panels (b) and (c) of the same Figure.

Asymmetric Model
We shall now consider an asymmetric bi-ion beam plasma mixture.Figure 10(a-d) illustrates the existence domain(s) for IA solitary waves associated with either of the (fast and slow) modes, depicted versus the beam speed U, for different values of the spectral index κ e .Four different roots exist in this case, as shown earlier.The positive and negative domains are not the same in this case, hence the existence regions found (numerically) for negative and positive values of M are not mirror-symmetric anymore: see Figure 10, where values of M on either side of the horizontal axis were considered.
The fast mode is shown in blue color while the slow mode is in red, in Figure 10.The solid curves correspond to M s , while the dashed curves correspond to M s .As discussed in the symmetric case earlier, the existence regimes of both fast and slow IA solitary waves shrink.Note that, in the Maxwellian limit (of very large κ e ), our results agree well with those in Verheest & Hellberg (2021), as expected.
Figure 11 illustrates the Sagdeev pseudopotential profile associated with slow IA solitary waves, for different values of the Mach number, in the case of an asymmetric beam pair.In the case of the slow mode, we see that the amplitude of solitary waves increases as the Mach number increases.This fact is reflected in the corresponding potential (pulse) and electric field profiles, shown in the latter two panels, in the same Figure.
Figure 12 portrays the variation in the Sagdeev pseudopotential profile describing IA solitary waves associated with the slow mode, for different values of the spectral index κ e .Larger amplitude solitary waves will occur for lower κ e , as also witnessed and discussed in the symmetric case earlier.(Also see panels (b) and (c), which confirm this trend.) In this section, we explored the existence domains of subsonic IA solitary waves and their characteristics in a more general way by adopting symmetric as well as asymmetric ionbeam models.Now, in the next section, we will compare the  11.Slow-mode solitons in the asymmetric model, for different M. (a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c), respectively, depict the (slow mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of the Mach number.We have taken U = 1.5, κ e = 3, δ 1 = 2/3, and δ 2 = 1/3 in this plot.The same curve style and color has been adopted among the three panels, as a guide to the eye).characteristics of Slow IA waves with the observation of ESWs in reconnection jets by Liu et al. (2019).

Comparison with Observations in Reconnection Jets in the Earth's Magnetotail
Relying on the above considerations, we have considered in our model cold counterstreaming ion beams with suprathermal electrons population in the background.The following parameters corresponding to reconnection jets are used as a numerical application, to corroborate our model: 1. the density of the first (H + ) ion beam: n 1 = 2.6 × 10 4 m −3 with drift speed u 1 = − 900 km s −1 ; 2. the density of the second (H + ) ion beam: n 2 = 0.9 × 10 4 m −3 with drift speed u 2 = 950 km s −1 and density of electrons n e = 3.5 × 10 4 m −3 with temperature T e = 2.86 keV.The sound speed of e-i plasma is C IA = 523 km s −1 and the Debye length λ De = 2.12 km.Hence, the normalized parameters of the reconnection jet are δ 1 = 0.74, δ 2 = 0.26, U 1 = − 1.72, and U 2 = 1.82.
The features of ESWs propagating antiparallel to the ambient magnetic field B that were observed by the MMS spacecraft in the reconnection jet region, as reported by Liu et al. (2019), were bipolar electric fields, E ∼ (5-30) mV m −1 , with positive potential f ∼ (50-200) V, velocity antiparallel to B = −650 km s −1 and width W ∼ 20 km.In Table 1, we have numerically computed the features of ESWs associated with reconnection jet parameters (Liu et al. 2019) for the four different modes that may exist-as discussed above-i.e., fast and slow mode(s) for parallel or antiparallel propagation with respect to the magnetic field.Table 1 illustrates that both the characteristics of fast and slow modes of IA solitary waves corresponding to parallel (antiparallel) waves with respect to B are in good agreement with observations by Liu et al. (2019).It is important to highlight the wave characteristics of subsonic waves for both cases: highly suprathermal case (κ e = 3) and quasi-Maxwellian case (κ e = 20) are in an acceptable range as observed by Liu et al. (2019).(Note, for rigor, that in the Maxwellian case, Lakhina et al. (2021) studied this model with thermal effects, yet even then our results are in good agreement).
The present investigation proposes an analytical mechanism for the generation of subsonic electrostatic solitary waves in reconnection jet sites in the Earth's magnetotail.In a nutshell, it shows that the drifting cold ion beams, along with suprathermal electrons, may play a major role in the evolution of subsonic waves in reconnection jets.
Despite the (limited) mismatches in the speed of the fast mode, both the predicted pulse amplitude and bipolar electric fields are acceptable in the range of magnitude, though notably in the direction parallel to B (notice that Liu et al. 2019 only captured signals of ESWs with antiparallel orientation).Other reasons for this (limited) discrepancy in the agreement may be (a) the waveforms received by different satellites were not Figure 12.Slow-mode solitons in an asymmetric beam plasma system, for different spectral index κ e .(a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c), respectively, depict the (slow mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of κ e .We have taken U = 1.5, M = 0.8, and δ = 0.5 in these plots.The same curve style and color has been adopted among the three panels, as a guide to the eye).
synchronized, (b) the ion beams observed may not be the primary beam, and (c) the missing oxygen ion species (Lakhina et al. 2021).
We draw the conclusion that our ab initio investigation of subsonic ESWs is in good agreement with the observations in the reconnection jet region(s) in Earth's magnetotail that were reported by Liu et al. (2019).

Conclusions
We have studied, from first principles, the occurrence of both fast and slow IA subsonic waves in a non-Maxwellian space plasma comprising two counterstreaming (drifting, at equal absolute speed U) ion beams with kappa-distributed electrons.Our analysis relied on the exact (nonperturbative) Sagdeevtype pseudopotential technique, to derive an energy balance equation, based on which we have delineated the existence domains for subsonic solitary waves.
The impact of the beam velocity U and of the spectral index (κ e ) on the existence of subsonic waves has been investigated in detail.Linear analysis reveals that subsonic waves are unstable when the beam velocity is less than a certain value (that actually differs, for different values of κ e ); in this regime, only conventional supersonic (fast) waves are formed.On the other hand, when the beam speed exceeds the threshold, either supersonic or subsonic waves may (co-) exist.
The combined parametric dependence of the beam and of the electron statistical distribution has been analyzed, with respect to subsonic solitary waves.Both symmetric and asymmetric (beam) models were considered, and the results were qualitatively not different.
Our theoretical investigation has successfully grasped the basic features of the recent observations of subsonic ESWs in the reconnection jet sites reported by Liu et al. (2019) in MMS measurements.We have shown that an adequate plasma model incorporating a suprathermal population of electrons can explain the observed generation of slow/subsonic ESWs in the magnetotail, both qualitatively and quantitatively.Such subsonic waves may thermalize the ion beams and give a new channel for ion heating inside the jet region, as suggested by Liu et al. (2019).
This investigation is, to the best of our knowledge, the first of its kind, in that it has focused on a realistic Space plasma situation where two counterstreaming plasmas collide, against a background of nonthermal (kappa-distributed) electrons; this is a ubiquitous scenario in the Earth's magnetosphere and may arguably also occur in other planetary magnetospheres, affected by the solar wind, as a constant source of streaming ions in addition to energized electrons.As discussed above, earlier modeling efforts have focused exclusively on supersonic (superacoustic) solitary waves.However, slow (subacoustic) solitary waves are apparently detected in satellite observations, and a theoretical framework for these was so far lacking in the literature.Our work was based on the main ideas introduced by Lakhina et al. (2021Lakhina et al. ( , 2020)), and later refined by Verheest & Hellberg (2021), who argued that the "slow" waves detected by Liu et al. (2019) were, in fact, subsonic electrostatic solitons, which may, in fact, exist in a beam-permeated plasma.The subsonic wave interpretation may arguably depend on the adopted reference frame: while the subsonic modes were seen in the laboratory frame, the mode may be seen as supersonic in the corresponding ion frame.
It may be added, for rigor, that we have relied on the working hypothesis the propagation was along the ambient magnetic field.Furthermore, magnetization was "neglected" in our model, as overall current neutrality was implicitly assumed.Whether this was satisfied by the actual measurements-within the actual error bars on the densities and currents/beams-was beyond our reach to verify; however, we note that Liu et al. (2019) stated explicitly that no currents were measured.
Our study clearly suggests that ESWs played a vital role in jet dynamics (Liu et al. 2019).Our findings will help to unfold the (mostly unexplored) dynamical characteristics of subsonic waves observed in the reconnection jets of Earth's magnetosphere (Lakhina et al. 2021;Liu et al. 2019).

Figure 1 .
Figure 1.ω F in the k-U plane corresponding to the fast IA mode (a) for κ e = 20 (quasi-Maxwellian) and (b)for κ e = 2 (strongly nonthermal case) and for fixed δ = 0.5.Note that only ω F > 0 is taken into account (as ω F < 0 will simply be symmetric, lying on the negative axis).

Figure 2 .
Figure 2. Growth rate Im(ω S ) in the k-U plane corresponding to slow IA (unstable) mode (a) for κ e = 20 (quasi-Maxwellian) and (b)for κ e = 2 (strongly nonthermal case) and for fixed δ = 0.5.Note that only ω F > 0 is taken into account (as ω F < 0 will simply be symmetric, lying on the negative axis).

Figure 3 .
Figure 3. Re(ω S ) in the k-U plane corresponding to the (real part of the) slow (stable) IA mode (a) for κ e = 20 (quasi-Maxwellian) and (b) for κ e = 2 (highly superthermal regime), for a fixed value δ = 0.5.Note that only Re(ω S > 0) is taken into account as Re(ω S < 0) is symmetrical in the negative axis.

Figure 4 .
Figure 4. Existence diagrams for a symmetric bi-ion plasma mixture.Existence domain of fast and slow modes of IA solitary waves vs. the streaming velocity (beam speed (U), for (a) κ e = 20; (b) κ e = 4; (c) κ e = 3; (d) κ e = 2, and for a fixed value of δ = 0.5.The fast mode is shown in blue color and the slow mode is in red.The solid curves correspond to the sonic boundary M s , while the dashed curves correspond to M l .The black dotted curve representing M = U separates the fast from the slow regime.

Figure 5 .
Figure5.Fast-mode solitons in the symmetric model, for different κ e in a beam-free plasma.(a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c), respectively, depict the (fast mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of κ e .We have taken U = 0, M = 1.1 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).

Figure 6 .
Figure 6.Fast-mode solitons in the symmetric model, for different M. (a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c), respectively, depict the (fast mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of the Mach number.We have taken U = 1.5, κ e = 3 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).

Figure 7 .
Figure 7. Fast-mode solitons in the symmetric model, for different κ e .(a) Pseudopotential curve profile S(f) vs. the electrostatic potential f.Panels (b) and (c),respectively, depict the (fast mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of κ e .We have taken U = 1.5, M = 2.3 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).

Figure 8 .
Figure 8. Slow-mode solitons in the symmetric model, for different M. (a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c),respectively, depict the (slow mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of the Mach number.We have taken U = 1.5, κ e = 3 and δ = 0.5.The same curve style and color has been adopted among the three panels, as a guide to the eye).

Figure 9 .
Figure 9. Slow-mode solitons in the symmetric model, for different κ e .(a) Pseudopotential curve profile S(f) vs. electrostatic potential f.Panels (b) and (c),respectively, depict the (fast mode) IA electrostatic potential pulse (f) and the bipolar electric field (E) profiles vs. the traveling space coordinate ξ, for different values of κ e .We have taken U = 1.5, M = 0.65, and δ = 0.5 in this plot.The same curve style and color has been adopted among the three panels, as a guide to the eye.

Figure 10 .
Figure 10.Existence diagrams for an asymmetric bi-ion plasma mixture.The existence domain of solitary waves associated with the fast and with the slow IA mode (in the asymmetric model) are shown vs. the streaming (beam) speed U for (a) κ e = 20; (b) κ e = 4; (c) κ e = 3; (d) κ e = 2, and for fixed values of δ 1 = 2δ 2 = 2/3, U 1 = − U and U 2 = 2U.The fast mode is shown in blue color and the slow mode is in red.The solid curves correspond to M s , while the dashed curves correspond to M l .