Evolution of Subsonic Shock Waves Associated with Reconnection Jets in Earth’s Magnetotail

Motivated by the signatures of nonlinear electrostatic waves observed by the Magnetospheric Multiscale spacecraft mission in reconnection jet regions of Earth's magnetotail, we have explored the dynamical features of ion-acoustic shock waves in the magnetotail. In this investigation, we have examined the dynamics and characteristics of ion-acoustic subsonic shock waves in non-Maxwellian space plasma comprising of two counterstreaming ion beams with suprathermal electrons, assumed to follow a kappa (κ) distribution. A reductive perturbation technique has been adopted to establish an evolution equation for small amplitude electrostatic shock structures. Importantly, subsonic waves only exist when the beam velocity exceeds a certain threshold, beyond which supersonic and subsonic waves may coexist. The combined effects of the beam velocity and the non-Maxwellian electron statistics have been analyzed to examine the characteristics of subsonic shock waves. Both symmetric and asymmetric (in relative beam density) models have been considered, leading to distinct possibilities in the evolution of subsonic shock waves. The findings of the investigation will help unfold the relatively unexplored dynamical characteristics of subsonic shock waves that may form and propagate in the magnetosphere.

The investigation of shock waves in collisionless plasmas has garnered significant attention over recent decades, both from experimental (Nakamura et al. 1999;Heinrich et al. 2009) and theoretical (Washimi & Taniuti 1966;Malfliet & Hereman 1996;Kourakis et al. 2012aKourakis et al. , 2012b;;Sultana et al. 2012) perspectives.The increasing level of detail provided by both laboratory measurements and space-based observations underscores the necessity for continuous refinement of the underlying theoretical framework.Many theoretical investigations in the literature have relied on either the Burgers equation or the hybrid Korteweg-de Vries-Burgers (KdVB) equation, which elegantly captures the intricate interplay between nonlinearity, dispersion, and dissipation effects in the generation and evolution of shock structures.Nonetheless, the diverse range of plasma environments where shock waves may arise has posed challenges in developing a unified theoretical description for these inherently nonlinear phenomena.
The presence of both fast and slow ion-acoustic (IA) modes has been confirmed in multi-ion (warm ion) plasma models utilizing linear wave theory, subsequently leading to intriguing nonlinear investigations of plasmas permeated by two opposing ion-beam flows (Lakhina et al. 2020;Verheest & Hellberg 2021).These studies, employing a multi-plasma-fluid framework, have revealed that beam-permeated plasmas may sustain not just the "traditional" supersonic electrostatic solitary wave structures, but also a subsonic IA nonlinear mode.This subsonic mode propagates at a velocity lower than the speed of sound but above a specific threshold (for the Mach number value), say M min (Lakhina et al. 2020;Verheest & Hellberg 2021).A noteworthy development is the observation that the amplitude of SWs demonstrates an intriguing growth trend as the Mach number decreases toward the critical value of M min , below which the existence of SWs becomes unfeasible.A recent study by Lakhina et al. (2021) proceeded further by investigating the generation mechanism of these waves in the context of reconnection jets.This investigation considered the influence of beam thermal effects, an aspect not covered in a prior study by Liu et al. (2019).That study assumed the presence of ion beams coexisting with Maxwellian electrons within the background plasma.Lakhina et al. (2020Lakhina et al. ( , 2021) ) employed a plasma-fluid model that featured two counterstreaming ion beams and Maxwellian electrons.They were pioneers in demonstrating the possibility of SWs existing at Mach numbers below a critical threshold, provided that specific criteria were satisfied.In a subsequent study, Verheest & Hellberg (2021) re-examined this same model, omitting thermal ion effects and shedding light on the pivotal role played by the velocity of the ion beams, rather than thermal factors, in the formation of subsonic SWs.Their meticulous investigation conclusively verified that these slowmode SWs indeed propagate at subsonic speeds in the laboratory frame, both in symmetric and asymmetric beamplasma models.
Various satellite missions have provided compelling evidence for the prevalence of energetic particles in a variety of space and astrophysical settings.In these scenarios, the electron velocity distribution exhibits a distinctive long-tailed pattern at high energies, giving rise to a significant suprathermal component.This distribution diverges from the conventional "textbook" model of a thermal Maxwell-Boltzmann distribution, as discussed by Livadiotis (2017Livadiotis ( , 2018)).Suprathermal particles have been documented in various locations, including Earth's magnetosphere (Feldman et al. 1975) and the auroral region (Mendis & Rosenberg 1994;Lazar et al. 2008), as well as in Mercury's magnetosphere, i.e., detected by MESSEN-GER mission data (Ho et al. 2011).The kappa distribution was initially introduced by Vasyliunas (1968) as a heuristic model to interpret data from OGO 1 and OGO 3 spacecraft within Earth's magnetosphere.Subsequently, kappa-type (non-Maxwellian) distribution functions have been applied to characterize suprathermal particle populations in the solar wind (Armstrong et al. 1983) and in planetary magnetospheres, e.g., Earth's, Saturn's, and Jupiter's (Leubner 1982;Singh et al. 2021Singh et al. , 2022Singh et al. , 2023;;Varghese et al. 2023).
Magnetic reconnection is a mechanism through which magnetic energy is converted into plasma kinetic energy, accompanied by changes in the magnetic field configuration (Yamada et al. 2010;Fu et al. 2017).In Earth's magnetotail, this process occurs when two oppositely directed plasmas become interconnected with magnetic field lines, leading to the ejection of these reconnected plasmas at high velocities from the reconnection site.This phenomenon is commonly referred to as a reconnection jet (Cao et al. 2013).It is now wellestablished that reconnection jets play a pivotal role in energizing plasmas in space and astrophysical contexts, including phenomena like solar flares, pulsar winds, and active galactic nuclei (Kirk & Skjaeraasen 2003;Chen et al. 2019;Liu et al. 2019;Lakhina et al. 2021).Within Earth's magnetotail, reconnection jets can give rise to various types of plasma waves and instabilities, in addition to ESWs.Liu et al. (2019) presented the first observational evidence, using data from the Magnetospheric Multiscale (MMS) spacecraft, that elucidated the development of ESWs associated with cold ion beams within reconnection jets in Earth's magnetotail.
The magnetotail region is extensively permeated by counterstreaming ion beams, exhibiting relative velocities of up to 2000 km s −1 , in conjunction with hot electrons.These ion beams, as highlighted in the study by Liu et al. (2019), can serve as a source of free energy, contributing to the generation of ESWs.Over the past decade, slow IA electron SWs have been documented within the plasma sheet boundary layer (Norgren et al. 2015) and at reconnection sites (Graham et al. 2015).More recently, Lakhina et al. (2021) investigated the formation of both fast and slow ESWs within reconnection jets, taking into account the presence of two ion beams, thermal effects, and Maxwellian electrons, as mentioned above.
The research presented in this paper was initiated in response to the aforementioned factors and marked an inaugural (unprecedented, to the best of our knowledge) investigation of subsonic shock waves within the reconnection jet region(s) of Earth's magnetotail (plasma).Such subsonic shocks were previously detected by various satellite observations (Zhou et al. 2018;Gingell et al. 2019), but the underlying theoretical setting has so far remained unexplored in the existing scientific literature.Motivated by these considerations, our study at hand focuses on the presence of two counterstreaming (cold) ion populations surrounded by suprathermal electrons.To identify and characterize the plasma regime(s) where subsonic shock waves may exist, our approach in this article relies on a nonlinear reductive perturbation technique, specifically the KdVB equation, as our main analytical tool.
To the best of our knowledge, this investigation represents a novel endeavor as it concentrates on a realistic scenario within space plasmas to predict subsonic shock-like structure occurrence.Specifically, we explore the collision of two counterstreaming plasmas in the presence of nonthermal (kappa-distributed) electrons, a prevalent condition widely observed in Earth's magnetosphere.This phenomenon may also have relevance in planetary magnetospheres (e.g., the Martian one) influenced by interactions with the solar wind, which serves as a continuous source of streaming ions and energized electrons.Previous research primarily centered on supersonic SWs, with subsonic wave solutions typically ruled out in a quiescent plasma without a beam (Dubinov 2009).

A Double Beam Plasma Model
We have considered a plasma consisting of non-Maxwellian electrons (mass m e , charge −e) and double 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 as we are dealing with proton beams, i.e., hydrogen nuclei.) The fluid model equations describe the plasma state in terms of the two proton fluid densities (n j ) and the fluid speeds (u j )-where j = 1 or 2, respectively, for each of the two proton fluids-consist of the continuity equation(s): and Poisson's equation: where f ˜is the electrostatic potential and ò 0 is the permittivity of vacuum.
The electron density entering Poisson's equation can be obtained by integrating the kappa velocity distribution (Hellberg et al. 2009) as ,0 3 2 At equilibrium, the charge neutrality condition is where n j for ( j = 1 or 2) denote(s) the unperturbed number density of the ions, respectively.For the sake of analytical convenience, we shall adopt the following normalization: n n n is a characteristic length inspired by the expression for the Debye length in electron-ion (e-i) plasma), and adopts the expression for the ion plasma frequency in e-i plasma).Furthermore, the ion viscosity is rescaled as ), while the drifting speed of the two ion fluids (i.e., the beams) is also normalized as u u C j j ,0 ,0 0 = ˜/ .Applying the above normalization scheme, Equations (1)-( 4) take the form where Note that charge neutrality at equilibrium imposes the condition: reflect the relative density ratio (with respect to the total ion density n 0 ) of the ion beams.

Linear Analysis
Before examining the nonlinear analysis of IA waves, we will study the linear stability 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 fourth-order polynomial equation, which gives two fast and two slow modes.In the limiting case, without streaming velocities, i.e., u 1,0 = u 2,0 = 0, the above relation becomes where c 1 was defined above.Note that the quantity c 1 is related to both the (kappa-dependent) charge screening length λ Debye and sound speed C s in suprathermal e-i plasmas, both of which scale as c , C 1 , as discussed by Kourakis et al. (2012a) and independently (for the Debye length) by Bryant (1996) and Livadiotis & McComas (2014).(As expected, c 1 tends to unity in the Maxwellian limit as κ → ∞ ).
Equation (10) has four roots, thus prescribing two fast IA modes (in addition to two slow modes), which in fact propagate 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 may now proceed by assuming a symmetric bi-ion plasma, i.e., δ 1 = δ 2 = δ = 1/2, with u 1,0 = − u 0 , and u 2,0 = u 0 .Substituting into the above dispersion relation, we are led to and the corresponding solutions are given by The subscripts F and S denote the fast and slow IA modes corresponding to the ± sign, respectively.A short algebraic calculation reveals that there is a threshold value in the beam velocity above which subsonic IA waves exist, while the interval below is an unstable region.This was discussed in recent work by Singh et al. (2023; see Figures 1-3 therein).

Nonlinear Analysis
The standard reductive perturbative technique, as originally introduced by Washimi & Taniuti (1966) can be employed at this stage.The method relies on introducing (independent) stretching coordinates as where ò = 1 is a small (real) parameter reflecting the strength of nonlinearity.Here, λ is the phase speed of the solution (wave, shock form), whose value is left to be determined later, eventually imposed by algebraic considerations.The dependent variables n j , u j , and f can be expanded as a power of ò: We have considered a weakly damped medium, assuming an infinitesimal (value of the) kinematic viscosity, as where η j0 is a finite quantity.The latter assumption is justified a posteriori by the fact that dispersive, nonlinear, and dissipative terms will arise on equal footing, in the perturbative analysis, at a certain order.
Substituting Equations ( 9)-( 15) into Equations ( 6)-( 8) leads to the following expressions, upon solving the first-order equations: n u u u 1 and 1 .17 An algebraic compatibility condition leads to the following expression for the wave phase speed λ: Note that the right-hand side (RHS) is equal to c 1 (essentially the signature of the kappa distribution in the linear results), which was defined previously.It is interesting to point out that, had one considered the dispersion relation (Equation ( 10)) in the infinite-wavelength (i.e., vanishing wavenumber) limit k → 0, one would have found precisely Equation (18), via a different approach.One concludes, therefore, that parameter λ represents the phase speed of low-k linear waves, i.e., essentially the true IA ("sound") speed in the given plasma configuration, as physically and mathematically expected.
Proceeding to the next higher order of ò, one finds the following equations: One may now eliminate the second-order quantities from Equations ( 19)-( 21) by simply using Equations (17)-( 18).
After some straightforward algebraic manipulation, one obtains a closed partial-differential equation (PDE) in the form of a hybrid KdVB equation: The coefficients appearing in the KdVB Equation ( 22) consist of the nonlinearity coefficient and the dissipation coefficient Recall that c 2 is a kappa-dependent parameter related to quadratic nonlinearity in f, which was defined earlier; see Equation (9) above.
It is interesting to point out that both the dispersion coefficient B and the damping coefficient C, given in Equations ( 24) and (25) above, respectively, may acquire negative values, due to an interplay between the values of the beam speeds u 1,0 and u 2,0 (and also η 1,0 and η 2,0 , as regards C).This new possibility, nonexistent in beam-free plasmas-since, upon setting u 1,0 = u 2,0 = 0, both B and C become positive definite-calls for a certain caution in the physical interpretation of analytical solutions for physically meaningful quantities, such as, e.g., f 1 , governed here by Equation (22) above.
In the investigation that follows, we have considered a pair of proton (H + ) beams, i.e., consisting of identical microscopic particles (ions).Therefore, we may set η 1,0 = η 2,0 = η 0 to proceed.The dissipation coefficient is thus given by C = η 0 /2 (>0).Note that C is positive, therefore, throughout the remainder of the paper.However, B may either be positive (as in the beam-free case) or negative, depending on the value (s) of the beam speed(s) entering the expression (Equation ( 24)) above.

Symmetric Model
Let us consider a pair of counterpropagating beams of equal density and (absolute) drift velocity, by taking δ 1 = δ 2 = 1/2 and u 2,0 = − u 1,0 = u 0 , along the lines introduced by Lakhina et al. (2020) and Verheest & Hellberg (2021).Equation (18) thus reduces to which can be rewritten as Two roots for λ 2 can be computed from the above biquadratic equation, in the form of From the above four roots thus obtained for λ F,S , we shall retain the (two) positive roots for numerical analysis.Indeed, we need not discuss the case λ < 0, as it just follows the same behavior as λ > 0 in the symmetric case.It is straightforward to show that the latter two expressions could also be recovered from the dispersion relation (Equation ( 13)) in the small wavenumber limit, i.e., by assuming w and setting λ = ω/k for the phase speed.
In this (symmetric) case, the nonlinearity, dispersion and dissipation coefficients in the KdVB equation can be respectively rewritten as Iterating an earlier comment, note that C > 0, while B may either be positive (as in the beam-free case) or negative, depending on the value(s) of u 0 entering the expression (Equation (31)) above.
which can be rewritten as This quartic polynomial equation in λ has four roots.In this (asymmetric) case, the coefficients in the KdVB equation can be rewritten as We notice that, as in the symmetric case above, C > 0 here too, while B may either be positive (as in the beam-free case) or negative, depending on the value(s) of u 0 entering the expression (Equation (36)) above.

Shock-wave Solution of the KdVB Equation
To obtain a stationary-profile solution of the KdVB equation, one may consider the variable transformation to the moving coordinate c x t = - (while the variation, i.e., the derivative with respect to time τ will be suppressed).Here, λ is the speed of the electrostatic shock wave in the moving reference frame (i.e., moving at the sound speed), normalized by C 0 ; in other words, the shock speed will be equal to v shock l = +   .Consider Φ = f 1 (ξ, τ), obeying Equation (22), which now becomes Φ(χ) in the moving reference frame, whose dynamics is the shock amplitude and is related to the shock width, i.e., its spatial extension (the larger the value of Δ, the more stretched the shock form will be).Furthermore, the shock speed  is prescribed, as 6C 25B .

= 
We should point out that the shock excitation obtained for the electrostatic potential will be associated with a monopolar (pulse-like) form expected to be observed for the electric field.
Recalling that E = − ∇Φ (in the electrostatic approximation), one obtains the following expression (in 1D geometry) for the E-field: represents the maximum/minimum (for A > 0 or for A < 0, respectively) value of the monopolar E-field form.The former case (E 0 > 0) will be the standard case for IA waves ("positive polarity") in monopolar E-field space observations, while the latter case (E 0 < 0) will be referred to as a "negative polarity" excitation.This prediction matches, qualitatively speaking, certain space observations of propagating monopolar electric field waveforms (Tsurutani et al. 1998;De Keyser et al. 2010).
It is interesting to point out that, since all of the coefficients (A, B, and C) may be either positive or negative (see the discussion above), their sign may affect the sign(s) of Φ m , Δ, and  .As a consequence, different possibilities exist, as regards the shock behavior.Recall that B > 0 in the beam-free case; therefore, some of the possibilities open in our case did not exist in a beam-free model.
Let us assume that Δ > 0. Note that the limit of RHS of Equation (39) for χ → m ∞ is Φ m and 0, respectively, in this case.Therefore, the latter expression represents a positivevalued function, monotonically decreasing from Φ m > 0 to zero, for AB > 0.
On the other hand, the inverse monotonic behavior will be witnessed for AB < 0: a negative-valued function is obtained that increases from Φ m < 0 to zero as χ varies from − ∞ to ∞.
Inversely, assuming that Δ < 0, the limit of RHS of Equation (39) for χ → m ∞ is 0 and Φ m , respectively, in this case.Therefore, the latter expression represents a positivevalued function, monotonically increasing from zero to Φ m > 0, for AB > 0. On the other hand, for AB < 0, a negative-valued function is obtained that decreases from zero to Φ m < 0 as χ varies from − ∞ to ∞.
Interestingly, the sign of B determines the sign of the propagation speed  also: a positive B value (as in the beamfree case u 0 = 0) implies a rightward propagating excitation: an anti-kink soliton form is thus obtained if A > 0, while a kink soliton is obtained in A < 0. Inversely, a negative B value implies a leftward propagating excitation: however, a stepshaped function (Φ m < 0) propagating backwards is tantamount to an anti-kink propagating to the left (this is obvious, upon a space-reversal).Likewise, a descending-step-like function (Φ m > 0) propagating backward is tantamount to a kink excitation propagating to the left.
The different possibilities that exist in terms of our plasma parameters will be exposed in detail in Section 5.

Parametric Analysis
Based on the above algebraic results, we may now carry out a parametric investigation to study the characteristics of IA shocks associated with fast and slow modes in non-Maxwellian plasma.

Symmetric Model
Fast mode. Figure 1 shows the dependence of (a) the phase velocity (λ F ) (i.e., the shock speed, to leading order), (b) the nonlinearity coefficient A, and (c) the dispersion coefficient B associated with the fast mode in a symmetric bi-ion-beam plasma, on the beam speed (u 0 ) and the electron spectral index (κ).Notice that all of the coefficients (A, B, and C) are positive in this case, which prescribes positive values of (all of) Φ m , Δ, and  (see the definitions above).Note that the phase speed increases with either a larger beam velocity (value) or with a higher spectral index.In other words, stronger deviation from Maxwellian statistics (i.e., lower kappa) leads to a slower acoustic ("sound")-as expected from earlier works (Kourakis et al. 2012a)-and, consequently, to the possibility for a slower shock speed.An analogous effect regards the beam velocity, which actually seems to energize the shocks: the faster the beam, the faster a shock may be.It turns out that the nonlinearity coefficient (A) is slightly higher (for finite beam speed) for smaller κ (i.e., for stronger deviation from the Maxwellian), and also higher for higher beam speed (see Figure 1 (b)).Conversely, the dispersion coefficient decreases (slightly) for larger u 0 , although it also decreases for smaller spectral index (value), as shown in Figure 1(c).Recall, in passing, that the dissipation coefficient in Equation ( 22) is directly proportional to the medium's viscosity (coefficient).
Figure 2 depicts an electrostatic shock profile associated with the fast mode, emphasizing its variation with different parameters.Note that only positive polarity (monopolar) electric field profiles associated with "shocks" (or anti-kink soliton forms) moving toward the right are predicted in this case, since both A and B are positive, for any u 0 and κ; see the definition of Φ m in Equation (39).It is clearly seen that the amplitude of IA shocks associated with the fast mode in the symmetric model decreases for higher values of either the spectral index and the beam velocity, as shown in panels (a) and (b), respectively.Lower values of kappa (i.e., stronger deviation from the Maxwellian behavior) will therefore result in stronger shocks.Note that larger amplitude shocks develop for higher values of viscosity coefficient, as seen in Figure 2(c): this is due to the amplitude being related to C 2 , which in turn is directly proportional to η. Panels (d), (e), and (f) in the same figure illustrate the variation of the monopolar electric field (E) pulse profile of shocks corresponding to panels (a), (b), and (c), respectively.
Slow mode. Figure 3 depicts the variation of ((a), (b)) the phase speed (λ S ), (c) the nonlinearity coefficient A, and (d) the dispersion coefficient B, for electrostatic shocks associated with the slow mode in a symmetric bi-ion-beam plasma.Notice that coefficients A and B are both negative in this case (while C > 0), which prescribes positive values of Φ m , but negative Δ and negative  (see the definitions above).Note that slow IA shocks only occur for values of the beam velocity exceeding a certain threshold value, as discussed in Singh et al. (2023).The phase speed increases for higher beam speed (value), and also for larger electron spectral index (kappa).Note that the nonlinear coefficient is negative for the slow mode.The absolute value of the nonlinear coefficient first decreases up to κ ≈ 4 and then increases (see Figure 3(c)).We also notice that the nonlinearity coefficient increases in absolute value as the beam speed increases.The dispersion coefficient is also negative and increases (in absolute value) for larger u 0 and the spectral index as shown in Figure 3(d).
Figure 4 illustrates the variation of an electrostatic shock profile associated with the slow IA mode.Only negative polarity monopolar electric field profiles and anti-kink soliton forms (for f) moving toward the left are observed, since B is negative in this case; see the definition of the propagation speed  following Equation (39).
One clearly sees that the amplitude of IA shocks associated with the slow mode in the symmetric model decreases for higher values of spectral index and beam velocity, as shown in panels

Asymmetric Model
The variation of the phase speed (for both fast and slow modes) versus the spectral index of electron (κ) for different values of beam velocity is shown in Figure 5. Recall that only the fast mode exists below a certain threshold in the beam speed (value), as shown in Figure 5; see panels (a) and (b) therein.Above this threshold, both fast and slow modes may occur in the κ domain.Conversely to the previous case (slow mode), the fast mode phase speed increases for lower κ (see Figure 6 shows the variation of the nonlinearity (A) and dispersion (B) coefficients versus the spectral index (kappa), for a fixed value of the beam speed.Note that both coefficients are negative for the fast mode and positive for the slow mode.
Fast mode.Notice that the coefficients A and B are both negative in this case (while C > 0), which prescribes positive values of Φ m , but negative Δ and negative  (see the definitions above).Figure 7 illustrates the variation of an IA shock profile associated with the fast mode.Only negative polarity electric field profiles associated with shocks (or antikink solitons) moving toward the left are observed, since B < 0 in this case; note the definition of the propagation speed  following Equation (39).We see that the amplitude of IA shocks associated with the fast mode in the symmetric model  (2019) in reconnection jet regions in Earth's magnetotail, we have considered two cold counterstreaming ion beams with non-Maxwelian electrons.We have adopted the plasma parameters reported in the latter reference.The concentration of the first proton beam I is n 1 = 2.6 × 10 4 m −3 , while its streaming velocity is u 1 = − 900 km s −1 .The concentration of the second proton is n 2 = 0.9 × 10 4 m −3 , with streaming velocity u 2 = 950 km s −1 .The electron concentration is n e = 3.5 × 10 4 m −3 , with temperature T e = 2.86 keV.The IA sound speed for these values (in e-i plasma) is C S = 523 km s −1 and the Debye length is λ De = 2.12 km.Thus, the dimensionless parameters in our model take the values δ 1 = 0.74, δ 2 = 0.26, U 1 = − 1.72, and U 2 = 1.82.Note that we have considered viscosity η 0 = 0.5 (as an arbitrary, indicative value).
The signatures of electrostatic waves 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), with potential f ∼ (50-200) V, velocity antiparallel to B = −650 km s −1 and width Δ ∼ 20 km.Table 1 presents the attributes of both fast and slow modes of IA shock waves corresponding to parallel (antiparallel) waves relative to the magnetic field.Our findings demonstrate good agreement with the findings reported by Liu et al. (2019).
In this work, we have introduced an analytical framework modeling the occurrence of shock waves, with implications in the reconnection jet regions in Earth's magnetotail.In summary, our observations suggest that streaming beams, suprathermal electrons, and fluid viscosity combined could bear a significant impact on the existence and also on the propagation characteristics of subsonic waves in the reconnection jet region of the magnetotail.Notwithstanding the (limited) disparity in the (value of the) velocity of the fast mode-in comparison with the observations reported by Liu et al. (2019) -the predicted amplitude and phase velocity falls within an acceptable magnitude range.This is true, particularly, in the direction parallel to the magnetic field.It is worth noting that Liu et al. (2019) exclusively detected signals of electrostatic nonlinear waves with an antiparallel orientation.
The results of this fundamental investigation, from first principles, is in good agreement with the observations of nonlinear excitations in the reconnection jet region(s) in Earth's magnetotail that were reported by Liu et al. (2019) and may be further confirmed by some future space missions.

Conclusions
We have investigated the existence domains for both fast and slow IA subsonic shock waves in suprathermal space plasma comprising two counterstreaming (drifting) ion beams with non-Maxwellian (suprathermal) electrons.We have employed a perturbative technique to derive an evolution equation in the form of a hybrid KdVB type PDE for the electrostatic potential.Our ambition was to shed some light on the role of the beam velocity and the spectral index on the evolution of subsonic shock waves and also on their potential existence in the magnetotail.
We have considered, separately, two cases, namely a symmetric beam pair, and also an asymmetric case (of two different beams).In both of these cases, for beam velocity above the threshold, both supersonic as well as subsonic waves may exist.Note that only positive polarity monopolar electric field profiles-associated with shocks of anti-kink soliton shape-moving toward the right (i.e., in the positive direction) are predicted, for the symmetric fast mode and-independently -for the asymmetric slow mode.Inversely, only negative polarity electric field (monopolar) profiles-associated with shocks of anti-kink soliton shape-moving toward the left (i.e., in the positive direction) are found for the symmetric case fast modes and also-independently-in the asymmetric case slow modes.The combined effects of the beam velocity (value) and electron superthermality on the characteristics of subsonic shock waves have been analyzed.Both symmetric and asymmetric models are considered to explore the evolution of subsonic shock waves in non-Maxwellian plasma.
The findings of this investigation should help elucidate themostly unexplored-characteristics of subsonic shock waves observed in the magnetosphere.Our investigation establishes the possible occurrence of subsonic shocks inside reconnection jets of Earth's magnetosphere (Liu et al. 2019;Lakhina et al. 2021;Singh et al. 2023), based on observations by Liu et al. (2019).Note.The concentration of the first proton beam I is n 1 = 2.6 × 10 4 m −3 , while its streaming velocity is u 1 = − 900 km s −1 .The concentration of the second proton is n 2 = 0.9 × 10 4 m −3 , with streaming velocity u 2 = 950 km s −1 .The electron concentration is n e = 3.5 × 10 4 m −3 , with temperature T e = 2.86 keV.The IA sound speed for these values (in e-i plasma) is C S = 523 km s −1 and the Debye length is λ De = 2.12 km.Note that we considered viscosity η 0 = 0.5.Here, f is the shock amplitude (in volts), V is the phase speed (in kilometers per second), and |Δ| is the width of shocks (in kilometers).

Data Availability
Any data that support the findings of this study are included within the article.No new data were created or analyzed in this study.

Figure 1 .
Figure 1.Symmetric case: fast mode.The variation of (a) the phase speed λ F , based on Equation (28), (b) the nonlinearity coefficient A, and (c) the dispersion coefficient B with respect to κ and u 0 .We have taken δ = 0.5.

Figure 5 .
Figure 5. Asymmetric case.The phase velocity of the fast and slow mode(s) vs. the spectral index (κ), for (a) u 0 = 0; (b) u 0 = 0.5; (c) u 0 = 1; and (d) u 0 = 1.5, as well as for a fixed value of δ 1 = 2δ 2 = 2/3.The fast mode is shown in blue color and the slow mode is in red.

Figure 6 .
Figure 6.Asymmetric case.The variation of the coefficients A and B vs. the spectral index κ for the fast and the slow mode, respectively, are shown.In this plot, δ 1 = 2δ 2 = 2/3, η 0 = 0.5, and u 0 = 1.5 are kept fixed.
4(a) and (b), respectively.Note that larger amplitude IA shocks will evolve for high values of the viscosity coefficient, as shown in panel 4(c).Panels (d), (e), and (f) illustrate the variation of the monopolar electric field (E) pulse profile associated with the shocks shown in panels (a), (b), and (c), respectively.

Figures 5
Figures 5(c) and (d)), i.e., for stronger deviation from the Maxwellian picture.Figure6shows the variation of the nonlinearity (A) and dispersion (B) coefficients versus the spectral index (kappa), for a fixed value of the beam speed.Note that both coefficients are negative for the fast mode and positive for the slow mode.Fast mode.Notice that the coefficients A and B are both negative in this case (while C > 0), which prescribes positive values of Φ m , but negative Δ and negative  (see the definitions above).Figure7illustrates the variation of an IA shock profile associated with the fast mode.Only negative polarity electric field profiles associated with shocks (or antikink solitons) moving toward the left are observed, since B < 0 in this case; note the definition of the propagation speed  following Equation (39).We see that the amplitude of IA shocks associated with the fast mode in the symmetric model

Table 1
Liu et al. (2019)cs of IA Shock Waves Have Been Computed by Using the Data Sets Listed inLiu et al. (2019)