Ion Acceleration in Multi-Fluid Plasma: Including Charge Separation induced Electric Field Effects in Supersonic Wave Layers

: The need to understand the process by which particles, including solar wind and coronal ions as well as pickup ions, are accelerated to high energies (ultimately to become anomalous cosmic rays) motivate a multi-ﬂuid shock wave model which includes kinetic effects (e.g. ion acceleration) in an electromagnetically self-consistent framework. Particle reﬂection at the cross-shock potential leads to ion acceleration in the motional electric ﬁeld and thus anisotropic heating and pressure in the shock layer, with important consequences for the multi-ﬂuid dynamics. This motivates development of a multi-ﬂuid model of solar wind electrons and ions treated as ﬂuid, coupled self-consistently with a small population of kinetically treated ions (e.g. pickup ions.) Consideration of both the time dependent and steady state regimes, indicate that such a multi-ﬂuid approach is necessary for resolving the, Debye scale, particle reﬂecting cross-shock potential and subsequent dynamics. To study charge separation effects in narrow, supersonic wave layers we consider a reduction of the system to the steady state for cold ions and hot electrons and ﬁnd two types of solitary waves inherent to the reduced two-ﬂuid system in this limiting case.


Introduction
Investigations of multiply reflected ion (MRI) acceleration at perpendicular shocks, where ions are treated kinetically (as test particles with detailed trajectories calculated from the Lorentz force equation) have found MRI to be an efficient mechanism capable energizing ions up to energies on the order of ∼ 100keV, as needed for injection into diffusive shock acceleration mechanisms [1][2][3][4][5] which further accelerate them to become the so-called anomalous component of cosmic rays. Much of this work treats the electromagnetic fields as fixed and ions as test particles. For an improvement, we have presented a two-fluid solar wind (SW) simulation with kinetically treated pickup ions (PUIs) coupled quasi self-consistently via inclusion of the PUI anisotropic pressure-calculated directly from their precisely tracked trajectories [6]. Here we consider a self-consistent inclusion of the back-reaction, due to ion acceleration, including their action on the electromagnetic fields as well as their anisotropic pressure. The possibility of self-consistently including kinetic particles using a 1-fluid MHD model, with source terms included to account for shock accelerated ions, has been rejected because, as considered below, only a multi-fluid plasma treatment can resolve important charge separation effects in supersonic wave layers.
When charge separation gives rise to particle reflecting electrostatic fields, resolution on electron Debye scales is necessary to include important plasma dynamics. For example, studies of adiabatic solitary waves [7] and shock surfing at a two-fluid plasma model [8], indicate that the cross-wave electrostatic field structures form on electron Debye scales. Such fine scale, electric fields provide important force sights where suprathermal ions gain energy by acceleration along the motional electric field after reflection from an electrostatic cross-shock potential (CSP), thus providing a dissipation mechanism for super-critical collisionless shocks, which can explanation some observations. For example, reflected and accelerated pickup ions can explain the cooler than expected solar wind observed by Voyager 2 downstream of the heliospheric termination shock [1,2,4,9,10].
Such observations and simulations indicate that electrostatic solitary waves (ESWs) are dynamically coupled to shocks in multi-fluid plasmas [11][12][13]. In this context, we pursue the development of a self-consistent fluid dynamics model, that can incorporate the effects of reflected particles, in particular for shock-surfing pickup ions. Indeed, our two-fluid hybrid model, including dynamically coupled ions, treated as test particles [6], required a multi-fluid treatment without the common assumption of quasi charge-neutrality (e.g. used consistently in [Tidman and Krall, 1971] [14]) since particle reflecting electrostatic structures arise from charge separation when ion and electron densities get perturbed from balance in narrow wave layers.
As an illustrative simplification of the complex multi-fluid system, we treat the case of cold ions and hot electrons in the steady state of the wave frame. Solving the resulting coupled set of ordinary differential equations (ODEs), via numeric parameterization of the set, we find both accelerating and decelerating solutions of electrostatic solitary waves (ESWs) associated with both electrostatic potential hills and wells, much like the two types of ESWs observed by Geotail near Earth's bow shock [11].

Motivation for the Development of a Multi-fluid Plus Kinetic Ions Plasma Model
Multiplying the continuity equations (see e.g. [15]) for electrons and ions by their fundamental electron charge e and summing yields the charge continuity relation ∂q ∂t + ∇ · j = 0 where q and j are the total charge density and current. The usual neglect of the displacement current in MHD (magneto-hydrodynamics) yields in Ampere's law implying that ∇ · j = 0 and thus that the charge density q is everywhere constant, meaning that a 1-fluid MHD approach, which assumes charge neutrality everywhere, cannot resolve the cross-shock potential (CSP) inducing charge separation in the shock layer. As outlined below, in order to resolve charge separation, we present a two-fluid plasma model which includes a third population, composed of candidate ions which we term as shock accelerated particles (SAPs) likely to be MRI accelerated, to be treated kinetically as quasi-test particles-meaning that for the duration of a time-step (on the plasma fluid scale) the electromagnetic fields are fixed and the trajectories of the SAPs, along with their associated stress and energy tensors, are calculated. The SAP tensors are then fed step-wise into the 2-fluid plasma model as source terms. A brief description of shock surfing acceleration, with implications for source terms follows. Multiply reflected ion (MRI) acceleration can occur when plasma is suddenly decelerated at a shock wave where the more massive ions on average, having more momentum, are able to penetrate further into the density pile-up immediately downstream of the shock than do the less massive electrons. The subsequent charge separation results in the formation of a cross-shock potential (CSP) which acts as an electrostatic barrier to incoming ions. Some ions incident at the shock with insufficient shock-normal component of velocity to overcome the CSP will be reflected. Reflected ions, which can be thought of as kicked out of the SW frame where electric fields are negligible, are 3 of 35 acted on by the motional electric field (directed tangentially to the surface of a perpendicular shock) and can thus be significantly energized by MRI acceleration. Illustration of the idealized initial ion distribution, in the shock frame, before and just after reflection. (a) The unreflected initial distribution is composed of a pickup ions (PUIs) which reside on a shell in velocity space of radius u 0 (the solar wind (SW) speed) and SW ions residing in a core located at (u 0 , 0, 0). The dashed portion of the PI shell to the left of v spec on the v x -axes (vertical dashed line) has an insufficient x-component of velocity to over come the cross shock potential (CSP, Φ) and will be reflected. (b) The ion distribution immediately following refection is composed of the specularly reflected fractional part of the PUI shell and the remainder which, along with the SW core, is convected downstream without reflection. Not illustrated by this figure is the stretching along motional electric field ( E = u × B) direction that the reflected partial shell will undergo during the energization process. It is evident that the reflected PUI partial shell has gained energy after a single reflection, with respect to the fluid frame where the SW core is at the origin.
As illustrated in Figure 1 (which illustrates a shell distribution where the entire portion of the shell satisfying v x < v spec = √ 2eΦ/m has undergone a single reflection) ion reflection can occur when i.e. when the ion's x-component of velocity is insufficient to overcome the CSP barrier. For a PUI distribution incident on a simple step-function shock structure with a vanishingly narrow ramp width (such as considered by [2]) the entire portion of the distribution satisfying the reflection condition will undergo one or more reflections, resulting in very strong energization of the downstream PUIs. On the other hand the main ion population, the SW ions which lie on a relatively cold Maxwellian core distribution (the heavy central dot depicted in Figure 1) are transmitted downstream of the shock with virtually no reflection. In a recent article [1] we have performed a test-particle simulation of the termination shock (TS) showing that MRI acceleration can also be important at shocks with fairly large ramp widths but that still posses sufficiently narrow fine structure to support reflection. Because of the important role the shock structure plays in ion acceleration, we seek a 2-fluid plasma model that will allow us to resolve fine details of the electric and magnetic fields in the local vicinity of the shock. There are several important features of MRI acceleration that the 2-fluid plus source terms model must include. For example, as depicted in Figure 3, ion reflection at the shock ramp is expected to perturb state variables such as pressure and density (represented by the arbitrary perturbation quantity g i (x, t) in Figure 3). Furthermore, our simulations [1] have shown that PUI pressure behind the termination shock can be highly anisotropic. . Depiction (for illustration purposes) of a perturbed quantity such as ion density, ion pressure, etc. represented by g i (x, t 0 ) at some arbitrary time t 0 shown together with the cross shock potential (CSP). The qualitative plots for these two functions are superimposed over a qualitative sketch of an ion trajectory. The illustration shows how some incoming ions are trapped by reflection from the CSP and accelerated in the motional electric field ( E = − u × B) before being convected downstream. This trapping leads to a small ion density increase which is accounted for by the nonzero value of g i in a small region centered around the shock front. Figure 4 illustrates ion distributions ahead of and just behind the TS model. Note that virtually none of the upstream Maxwellian SW distribution lies to the left of the velocity v spec , whereas a significant portion of the upstream PUI filled-shell does lie to the left of v spec (i.e. in the reflecting region of velocity space). We have shown, in test-particle calculations, that (cold) Maxwellian distributions incident on a narrow, perpendicular shock pass downstream virtually unreflected, as expected. This justifies the assumption employed below in the construction of our numerical model that the 2-fluid SW is Maxwellian and isotropic everywhere. The right hand panel of Figure 4 is the phase-space density plot of an MRI accelerated PUI distribution just behind the TS-model and it is clear that this distribution is anisotropic with important non-zero off-diagonal terms in the associated PUI pressure tensor. respectively, which shows MRI enhancement of the p a yy component as well as the formation of shear pressure terms p a xy = p a yx . We seek to include the formation of an anisotropic ion distribution in the shock layer via the inclusion of source terms in a two-fluid plasma model as described below.

Two-fluid Solar Wind Plus Kinetically Treated Shock Acceleration Candidate Kinetic Ions Numeric Model
We develop a quasi-self-consistent system composed of a core two-fluid plasma solar wind (SW) background, which is acted on by a small population of ions to be treated kinetically though self-consistantly coupled to the core SW. We are interested in the case where a small, energetic population (denoted as shock accelerated particle (SAP) candidates) undergoes multiply reflected ion (MRI) acceleration due to reflection from a cross-shock potential leading to acceleration in the motional electric field. A two-fluid description is derived to calculate the back reaction on the fields, in essence creating a feedback loop between kinetic shock accelerated particles (SAPs) with SW electron and ion fluids.
There are two distinct characteristic scale regimes involved in this development: long time and length scales (δt, L) associated with the description of the two-fluid core SW plasma, and short time and length scales (δt kin , l) associated with SAP dynamics (particularly within the narrow shock structure). The prescription here, for constructing a closed system that can be solved numerically, is to treat the electric and magnetic field E and B as constant for each time period δt during which we solve the motion of the kinetic ions (treated in a manner much like particle in cell (PIC) code); the pressure and heat flux tensors resulting from the detailed ion trajectories are then feed into the 2-fluid plasma equations as source terms; the 2-fluid equations are used to update the values of E and 7 of 35 B over the longer time scale and the process is iterated. Below we present a system ready for such a numerical treatment.
The derivation of the plasma fluid equations is followed from Boyd and Sanderson, Plasma Dynamics [15], Chapter 3 (from now on referred to as B&S), however instead of treating the system as two groups of particles, ions and electrons (labeled in B&S by '+' and '-' respectively), we treat the system as composed of three groups of particles, labeled {s, a, e}, which are solar wind (SW) ions, shock acceleration candidate ions and electrons respectively. The shock acceleration candidate ions (e.g. pickup ions (PUIs)) will compose between 0 and 20% of the total ion population and are kept separate since they will be treated kinetically. Definitions: Ion number density composed of the SW ions plus SAP candidates is and the electron density is n e . The charge density is q = en + − en e (2) and ion fluid velocity is u + = mn s u s + mn a u a mn s + mn a The current density is j = en s u s + en a u a − en e u e Here m is the ion mass (which for now is a proton mass), e is the electron/proton charge. The random (thermal) velocity of each species is Pressure tensors are defined as p ij = mn c i c j (6) and the heat tensors are q ijk = mn c i c j c k where, as per the B&S convention, species labels are suppressed to simplify the notation. Since, in (5), both SW ions and SAP ions have thermal velocities defined with respect to the common ion bulk velocity, the total ion pressure tensor can be written as the sum of two partial pressures We make the closure assumption that SW ions and electrons can be treated everywhere as thermalized, isotropic distributions such that and q s ijk = q e ijk = 0. 8 of 35 The shock accelerated particle (SAP) pressure is taken to be composed of an isotropic and anisotropic part, p a ij = δ ij p a + π a ij . (10) and has a nonzero heat tensor q a ijk . All fluid equations are constructed from the transfer equation where the closure assumption that collision terms are negligible has been made. Continuity equations result from (11) when ψ = 1.
∂n e ∂t + ∇ · (n e u e ) = 0 (species labels s, a suppressed) Electron momentum (i'th component) where for electrons v i v j = p ij mn + u i u j Taking m e → 0, (15) becomes Ion energy ∂ ∂t where for ions we have used where for electrons since q e ijk = 0. Taking m e → 0, (20) becomes ∂ ∂t With the Maxwell curl equations, equations (12), (13), (17), (18), (21) plus (22) and (23) form a set of 16 equations in the 16 unknowns {n + , n e , u + , u e , p s , p e , E, B} to be solved numerically by a fluid code, and the included sources terms p a ij and q a ijk are to be determined by a kinetic treatment of the SAP candidate ions.

Conservation form of the two-fluid plasma equations
Summing the momentum equations (13) and (17) yields where we have used and the Maxwell stress tensor is defined as: Summing the energy equations (18) and (21) yields ∂ ∂t where we have used An alternate form of the electron energy equation can be formed by taking (17) multiplied by u ei and subtracting the result from (21) yielding

Dimensionless forms of the equations
To prepare the equations for numerical discretization define the following dimensionless variables: where L is an arbitrary length scale and, as justified below, and the ∞ label indicates the constant far upstream (boundary) values for the state variables. Also note that the dimensionless ion pressure can be expressed in terms of ion partial pressures The dimensionless derivatives of dimensionless state variablesψ(x,t) = ψ/ψ ∞ take the form To simplify the notation for the two-fluid plasma system expressed in dimensionless form, it is useful to define the following dimensionless plasma parameters which are expressed in terms of the far upstream ion gyro-frequency, ion plasma frequency and ion gyro-radius defined as respectively. In writing the following dimensionless equations it is convenient to drop over-bars and take it as understood that all variables and derivatives dimensionless. Also note that the components of the dimensionless bulk velocities for the ions and electrons are given by respectively. In dimensionless form the continuity equations (12) become Ion momentum (13) becomes, (i'th component) Electron momentum (17) becomes, (i'th component) Ion energy (18) becomes ∂ ∂t Electron energy (21) becomes ∂ ∂t The the total momentum equation (24) becomes where the dimensionless Maxwell stress tensor is

The total energy equation (25) becomes
∂ ∂t The Maxwell equations (22) and (23) can also be expressed in dimensionless form as where The Poisson's equation (∇ · E = 4πq) can also be expressed in dimensionless form as

Special Case: Quasi-One Dimensional, Exactly Perpendicular Shock
As an initial case, limit variations to the shock normal direction (x-direction) so that Consider an exactly perpendicular shock configuration, with where the plasma flow at the far left (upstream) boundary is The assumption of charge neutrality can be imposed at the left boundary so that the dimensionless boundary conditions at x → −∞ are The fluid plus Maxwell's equations can now be expressed in the dimensionless form where the state vectors in (44) are The two remaining equations are from the electron momentum (33). Thus equation (44)  on the z-components of velocity, we assume w = 0 and W = 0. In fact examination of the steady state system (presented below) indicates that if W = w = 0 at the left boundary (x → −∞) then the z-comonent of bulk velocity will remain zero, i.e. we can view the z-direction as ignorable for the special case under consideration.
An alternate 'quasi-conservative' form of system (44) (plus (45) and (46), expressed as can be formed by swapping out the ion momentum and energy equation with the total momentum and energy equations in the conservation forms given by (36) and (37), where state vectors in (47) are written as The above system must be supplemented by the electron equations where (50), the quasi-conservative form of the electron energy, results on summing the x-component electron momentum (45) and the electron energy equation. Additionally the system (47) must be supplemented by the remaining non-conservative components of Maxwell's written as

Steady state equations for 1D case with perpendicular magnetic field
In the steady state ∂/∂t → 0 and (47) reduces to From which it is seen that several integrals of the motion are readily obtained. It turns out that for the special case under consideration there are nine algebra relations (i.e. 9 integrals of the motion) that can be obtained from the plasma system in the steady state (SS). The variables of the special case system under consideration are (n, N, u, U, v, V, p, P s , E x , E y , B z ) for a total of 11 unknowns. This means that, in principle, the SS system can be determined from a coupled pair of ODEs in two of the unknowns. Having solved the coupled ODEs the remaining unknowns can be found from the algebra relations and the system thus completely determined. The first four algebra relations are The first two SS relations, expressed in (54) above, follow from the continuity relations (i.e. first two components of (53)) and the boundary conditions (43). The remaining two constants of the plasma system, expressed by (54) follow from the Faraday relation given by the last component of (53) along with the y-component of electron momentum expressed by (49). The third, fourth and fifth components of (53) express conservation of total momentum and energy in the SS as where (43) and (54) have been applied and P ∞ = P ∞ s + P ∞ a + p ∞ is the total (dimensionless) pressure at the upstream boundary which is taken to be completely isotropic there. The equation of state for adiabatic electrons follows from (50) and (54). Using algebra relation (56) in Poisson's' equation (40) yields which, upon using (54), is also given by the y-component of the ion moment-i.e. the fourth component of equation (53). This demonstrates explicitly that the Poisson's equation is a corollary of the plasma fluid equations and the Maxwell curl equations. Several additional relations can be obtained by systematic substitution of the algebra relations (54), (55), (56) and (57) back into the SS forms of the non-conservative equations given by (44). As mentioned this process recovers the Poisson's equation (59) above along with several ODEs in the electron density from which we can recover an additional algebra relation By using (55) and (57) to find U = U(n, V), the coupled set of ODEs (65) can (in principle) be solved numerically and the steady state solution thus determined.

Steady state, 1D two-fluid plasma system with cold ions
Though reduced to the steady state, the above coupled set of ODEs, nonetheless, represents a considerable challenge for use in evaluating valid solutions. Thus, as an initial step, we consider a further simplification to the case of cold ions and hot electrons. Previously we have solved such a simplified two fluid system as an integral over a single component of the plasma velocity [16], which has the advantage of providing an exact (albeit numerical) solution of the wave structure. Here we solve for the wave structure by introducing a parameterization of the coupled ODEs. Though this coupled ODEs method might not be exact as the integral form, it has the advantage of being extensible for use in more complex cases, such as the coupled set of ODEs presented in the previous section. Furthermore, unlike the integral solution which only resolves decelerating solutions, the parameterization method used here finds both decelerating and accelerating solutions (termed here as lower and upper solutions) of the two fluids associated with potential wells and spiked potential hills of the electrostatic potential. This is an important improvement over the integral form, in the sense that the parameterized solutions can explain Geotail observations of ESWs associated with plasma flows beaming in both the upstream and downstream directions from the bow-shock [11].
The solutions of the two-fluid plasma system presented below show that the plasma reacts vary strongly to perturbations from an upstream neutral background state (NBS) where complete charge neutrality and zero electric field is assumed. The plasma response to perturbation from the NBS is the formation of solitary waves over which the electric field (induced by strong charge-imbalance) prevents the runaway acceleration of fluid velocities. Two types of solitons are identified. The classical case, where both proton and electron velocities dip below then return to their NBS values, is labeled as a lower-solution. The lower-solution is considered in many papers including Zank and McKenzie, 1988 [17] and McKenzie 2002 [18] as well as [16], and is characterized by an upward electrostatic potential hill. As we demonstrate below, setting the electron mass to zero is an excellent approximation for lower-solutions and has little effect on either the strength or structure of soliton waves in this case. The second solitary wave type presented his is labeled as an upper solution, where both proton and electron velocities accelerate above then return to their NBS values. Unlike lower waves, upper waves depend critically on a nonzero electron mass to prevent unbounded wave growth.
We consider the steady state, two-fluid plasma system, viewed in the wave frame such that ∂/∂t → 0, where ion pressure and magnetic field are taken to be zero. The 1D continuity equations can be expressed as n p u p = n e u e = n 0 u 0 The energy equation for adiabatic electrons is where the 0-subscript implies a constant, neutral background state (NBS) taken to exist somewhere upstream (as x → −∞). The system can be put in dimensionless form by introducing N = n p /n 0 , n = n e /n 0 , where all variables (u, U, p, E x , x) are now dimensionless. In dimensionless form the continuity equations (66) simplify to Equations (68), (67) and (69) using (70) can be put in the dimensionless form: where Taking the arbitrary length-scale factor to be the Debye length l = λ D , where λ D = kT e0 4πn 0 e 2 and setting where M e0 = u 0 /c e0 is the constant upstream, electron sound speed Mach number. The dimensionless, coupled set of ODEs (79), (78) and (80), contain the conservation of total momentum and total energy integrals, which can be expressed as which expresses conservation of momentum, and which expresses conservation of energy. Note that in the above conservation laws (81) and (82) we ½ assume that the undisturbed upstream plasma is charge neutral, so that the steady state system

½ ¿
To understand the behavior of the system near this fixed point we examine the first order expansion of the SSS about the NBS, which yields where the origin of the state variables has been shifted to the NBS (i.e. U → U − 1, u → u − 1, Ex → Ex − 0). The solution of linearized system (83) can be constructed from the Eigen-system for the 3x3 matrix and has the generic form where w = (u, U, E x ), the c i are constants determined by initial conditions and λ i and p i are the eigenvalues and associated eigenvectors for the 3x3 matrix. The first eigenvalue appearing in equation (84) is λ 1 = 0 which has the associated eigenvector p 1 = (1, 1, 0); thus the NBS as a constant state is built into the two-fluid system and occurs whenever initial conditions require that c 2 = c 3 = 0. In other words, if the NBS is chosen as initial data, the solution for the SSS is just the constant NBS. The remaining two eigenvalues for (83) can be expressed as from which we see λ ± will be pure real whenever and pure imaginary otherwise. Thus the linearized SSS is characterized by exponential growth or bounded oscillation depending on the value of the Mach number.
To facilitate the numerical solution of the cold ion SSS, it is convenient to introduce a parametrization variable t (which is not the time) and reduce the SSS to the form The SSS, as represented by equations (87) through (90), has a fixed point at u γ+1 = 1/M 2 e0 , E x = 0. However expansion of the SSS about this fixed point, in consideration of the integrals (81) and (82), indicates that a constant state is the only solution for a system which must include this fixed point. The other fixed point of the SSS occurs at u γ+1 = 1/M 2 e0 , U = 0, which, since here ions are perfectly cold, could be considered a limiting physical state of infinite ion density. To understand the interdependence between the electron and ion velocities u and U, it is useful to divide equation (88) by (87) which yields Introducing the parametrization variable s, equation (91) can be expressed as The interdependence between the ion and electron velocities is presented graphically in Figure  5. A key feature to note here, for the nonzero electron mass case, is that the level curves form closed-loops which confine the fluid velocities to a finite range. The closed-loop curves, of the left hand panel, are plots of U = U(u) as expressed by equation (82). The three curves plotted correspond to three choices of electron sound speed Mach number: (the unphysical) M e0 = 2 for the blue curve, M e0 = 0.15 for the orange curve and M e0 = 0.1 for the red curve. The solid portions of the curves correspond to regions where E 2 x > 0, as expressed by equation (81), which we take as a necessary condition for the solution to be physically real. The dashed portions of the curves correspond to (non-physical) regions where E 2 x < 0 and/or U < 0. Overlaid on the U(u) curves in the left panel of Figure 5, is the direction field of the solution space for the coupled ODE set (92) displayed as swirling red arrows. This direction field, corresponding to the choice M e0 = 0.1 (i.e. the red U(u) curve) circles around the stationary point (U = 0, u γ+1 = 1/M 2 e0 ) in a fashion characteristic of harmonic oscillator type ODE systems. To prove that this stationary point is a fixed node, meaning that no solution curves pass through this point, a linear expansion of the ODE set (92) about the stationary point can be effected, which yields pure imaginary eigenvalues of the linearized system; thus the system has a harmonic oscillator character and the fixed point is a node. The direction field corresponding to the linearized ODE system for (92) is presented in the right panel of Figure 5. Evidently, solution curves arbitrarily close to the fixed point (U = 0, u γ+1 = 1/M 2 e0 ) still 'circle' around it in oscillator fashion, demonstrating graphically that the fixed point is a node.
In the numerical solutions of the SSS presented, we limit ourselves to a consideration of solutions connected to the NBS. As noted above, if the NBS is chosen as the initial data, the only solution is the constant NBS. To find non-trivial solutions, we look for initial data nearby and connected (via equations (81) and (82)) to the NBS in the solution space. x > 0, then only curves which are concave up where they touch the NBS (e.g. the thick red curves) can be associated with physical solutions. This limits the range of Mach numbers that can support the physical connection of the plasma to the NBS to α 0 /(α 0 + 1) < M e0 < 1. Figure 6 presents plots of curves E 2 x = E 2 x (u) which simultaneously satisfy both equation (81) and (82) and helps to illustrate graphically how initial data can be chosen that is near and connected to the NBS in the domain of physically real solutions. The curves of the left panel correspond to M e0 = 0.1 (thick red curve), the limiting case M e0 = 1 (orange curve) and M e0 = 2 (dashed blue curve). The curves of the right panel correspond to M e0 = 0.03 (thick red curve), the limiting case M e0 = α 0 /(α 0 + 1) = 0.0233304 (orange curve) and M e0 = 0.01 (dashed blue curve). Assuming that physically real solutions require E 2 x > 0, then initial data should be chosen from curves E 2 x = E 2 x (u) that are concave up where they connect to the NBS. Evidently the condition for finding initial conditions associated with physically real solutions is which is identical to the condition (86) which specifies the region in which the SSS linearized about the NBS posses pure real eigenvalues. In other words if E 2 x > 0 when E x ∼ 0 then the eigenvalues λ ± , associated with equation (83), are real-valued, yielding a linearized solution space characterized by exponential growth. A further implication of the identical conditions (86) and (93) is that the oscillatory (Langmuir wave type) solutions, associated with imaginary λ ± , are unobtainable unless we relax the condition that the initial data must be chosen such that solutions pass through the NBS. We present numerical solutions of the SSS where initial data is chosen such that solution curves are physically connected to the NBS by ensuring that the condition (93) is satisfied. Two distinct classes of solutions (or solution types) become apparent: One type of solution is denoted as an upper solution, meaning that the plasma velocity at the wave center has increased from its initial value. The other solution type presented is denoted as a lower solution, where the plasma velocity decreases from its initial value at the center of the wave.
An upper-solution of the SSS equations (87) through (90), for the choice M e0 = 0.571631, is plotted in Figure 7. The solid blue curves are parametric plots of the form (x(t), u(t)). The dashed red curves are plots of the corresponding solution (84), w = w(x), of the linearized SSS (83). As Evidenced by inspection of Figure 7, the state variables (u, U, E x ), for the nonlinear case, are double valued in the region on the x-axis sandwiched between the sonic points where u γ+1 = 1/M 2 e0 and ∂u/∂x → ∞ (see equation (78) Inspection of Figure 8 indicates how a unique solution can be constructed from the double valued parametric solution. The left panel of the figure shows the solution curve for (x(t), u(t)) (blue curve) plotted on the same axis with (x(t), U(t)) (black curve-nonphysically scaled-up for visual clarity). The vertical, dashed, green lines are tangent to the u-curve at the sonic points where u γ+1 = 1/M 2 e0 and u ′ (x) → ∞. As required by equation (82) (and illustrated graphically in Figure 5), the sonic points are also locations where U = U(x) must reach a maximum. Note that at the bisector between the sonic points the slopes u ′ (x) = 0 and U ′ (x) = 0, which, by inspection of the SSS, requires that E x = 0 at the bisector also. The value of x = x c at the bisector (or wave center) can found from the plot of x = x(t) (blue curve) shown in the right panel of Figure 5). This plot shows that x is everywhere increasing with decreasing parametrization variable t, except in the region near the sonic points, where x(t) is double valued. To construct the unique solution, note that a horizontal line x = x c drawn between the two sonic points (dashed, black, horizontal line of the right panel) separates the double valued x(t) regions into two lobes of equal area. This suggests that (analogous to Van der Waals isothermal construction) a single-valued solution can be formed by removing the lobes and joining the x(t) curve along the horizontal, dashed, line segment. This equal area construction requires that the line segment be inserted precisely at the bisector value of x = x c , directly in the center of the symmetrical wave structure and effectively 'scissor-cuts' the top, double-valued lobes off of the u and U curves and joins the lower parts of the curves together with an infinitesimally narrow horizontal line along which the slopes u ′ (x c ) = U ′ (x c ) = 0 and E x (x c ) = 0. This solution is unique since there is only one point on the x-axis where the segment is inserted such that it is entirely compatible with the original SSS equations. Note that at the point where the horizontal segment is inserted the charge density of the plasma behaves like a delta function since dE x (x c )/dx → ∞. 2, shows a close-up view of u = u(x) near the wave center, revealing the double valued character of the solution. In these double valued, decelerating cases, a unique solution can be inserted using the method described above, in which the value of x = x(t) at the wave center is found from an equal area construction (see Figure 8).

Setting the electron mass to zero
To understand the significance of electron mass in the SSS, we can consider the limit m e → 0, which effectively removes the electron-acoustic, sonic point from the system. The zero electron mass form of the SSS can readily be obtained by keeping only terms with a 1/M 2 e0 factor in the momentum equations and by dropping all terms with an α 0 factor from the conservation equations. This yields the m e = 0 governing system of ODEs which contain the integrals where the parameter α 2 1 is here most conveniently expressed in terms of the 'collective ion-acoustic Mach number' citeMcKen02 The m e = 0 SSS (94) through (96) again is stationary at the NBS.

¾
The SSS for m e = 0, linearized about the NBS, is As with the previous m e = 0 case, the solution of equation (99) has the general form (84) and also possesses the eigenvalue λ 1 = 0 with the associated eigenvector p 1 = (1, 1, 0), with the implication that the NBS as a constant state is built into the system. The remaining two eigenvalues for (99) can be expressed as thus the eigenvalues λ ± will be pure real whenever and pure imaginary otherwise. With m e = 0 equations (91) and (92) become and Figure 18, for m e = 0, is analogous to Figure 5 for the nonzero electron mass case. Comparison of these analogous Figures shows that a key difference between the zero and nonzero electron mass cases is that when m e = 0 the sonic point is effectively shifted to infinity on the u-axis. Thus the proton verses electron velocity dependence, which is characterized by closed-loop (oscillator type) curves in phase-space for nonzero electron mass, become open (on the right), unbounded electron velocity curves when m e = 0. Figures 18 and 19 express, graphically, the condition M ep > 1 required for solutions to connectible physically to the NBS. Figure 20 is the only upper-solution presented for m e = 0. As noted above, setting the electron mass to zero has effectively shifted the sonic point to infinity, meaning that the upward growing electron velocity is unbounded and thus unphysical in this case.
Inspection of Figures 20 through 24, corresponding to solutions of the SSS for m e = 0, indicates that the decelerating soliton-type waves do not vary significantly between the m e = 0 and m e = 9.1094 * 10 −28 g cases.   (84)). This figure illustrates that setting m e = 0 has effectively shifted the sonic point to infinity on the x-axis with the result that plasma acceleration is unbounded when the velocity growth is positive.

Conclusion
The need for a self consistent system of multi-fluid equations which dynamically includes a small population of hot, energetic particles, such as pickup ions (treated kinetically as is done in PIC codes) has been addressed. We previously presented a multi-fluid treatment which included a population of kinetically treated pickup ions (PUIs) coupled to a two-fluid treatment with quasi self-consistent coupling through the anisotropic PUI pressure tensor [6]. Here we considered an improvement of that model with a small population of kinetic ions included in the multi-fluid system by a self-consistent treatment of their charge density and current component effects on the electromagnetic fields.
By considering the reduction of the multi-fluid plasma system to the special case of hot electrons and cold ions in the steady state, the existence of two types of electrostatic solitary waves (ESWs) were shown, characterized and typed here as upper and lower solutions. The character of the upper and lower solutions of the coupled set of ODEs are consistent with the two types of ESWs observed by Geotail near Earth's bow shock [11], with one type of ESW associated with electrostatic potential hills and the other with potential wells.
We have shown that setting m e = 0 (a typical approximation for PIC as well as other simulations) recovers classical solitons, labeled here as lower solutions, such as presented in the well known text book on plasma physics by Chen [19].
A type of soliton, labeled as upper solutions, was found that depends critically on a nonzero electron mass to prevent unbounded acceleration of the electron fluid velocity. These so-called upper solutions are consistent with the second type of ESW observed by Geotail [11] associated with 'positive' electrostatic potential hills and components of two-fluid plasma flow returning upstream from the ESWs towards the bow shock. The classical type solitons, or lower solutions, were found to have an electric field strengths that grow linearly with increasing Mach number. Lower solutions are smooth functions and straight-forward to evaluate over the Mach number range 1 < M ep < 1.8 (for γ = 5/3). However for all upper solutions and for lower solutions where M ep ≥ 1.9, the insertion of a unique solution is required.
A method for constructing unique solutions was outlined, and applied to both upper and lower types, where at the insertion point on the x-axis, the charge density is viewed as a delta like function where dE x /dx → ∞ and the electric field is forced to zero by insertion of a unique solution, with the result that a bounded, single-valued solution is formed that is fully compatible with the original governing system of ODEs.