The Role of Electrons and Helium Atoms in Global Modeling of the Heliosphere

We present a new three-dimensional, MHD-plasma/kinetic-neutrals model of the solar wind (SW) interaction with the local interstellar medium (LISM), which self-consistently includes neutral hydrogen and helium atoms. This new model also treats electrons as a separate fluid and includes the effect of Coulomb collisions. While the properties of electrons in the distant SW and in the LISM are mostly unknown due to the lack of in situ observations, a common assumption for any global, single-ion model is to assume that electrons have the temperature of the ion mixture, which includes pickup ions. In the new model, electrons in the SW are colder, which results in a better agreement with New Horizons observations in the supersonic SW. In the LISM, however, ions and electrons are almost in thermal equilibrium. As for the plasma mixture, the major differences between the models are in the inner heliosheath, where the new model predicts a charge-exchange-driven cooling and a decrease of the heliosheath thickness. The filtration of interstellar neutral atoms at the heliospheric interface is discussed. The new model predicts an increase in the H density by ∼2% at 1 au. However, the fraction of pristine H atoms decreases by ∼12%, while the density of atoms born in the outer and inner heliosheath increases by 5% and ∼35%, respectively. While at 1 au the density of He atoms remains unchanged, the contribution from the “warm breeze” increases by ∼3%.


Introduction
The heliosphere is formed due to the interaction between the expanding solar wind (SW) and the local interstellar medium (LISM; Parker 1961). Interstellar neutral (ISN) atoms, among which hydrogen and helium atoms constitute the most abundant species, play a key role in this interaction (e.g., Blum & Fahr 1969;Wallis 1971Wallis , 1975Baranov et al. 1979Baranov et al. , 1991Baranov & Malama 1993;Williams et al. 1995;Zank et al. 1996;Pauls & Zank 1997;Pogorelov et al. 2006Pogorelov et al. , 2007Frisch et al. 2015;Zank 2015;Bzowski et al. 2019). The heliopause (HP) is a tangential discontinuity that separates the ionized component of solar origin from the interstellar plasma. However, neutral atoms can enter the heliosphere and a fraction of them is ionized along the path, which provides a fundamental coupling between the SW and and LISM. Resonant charge exchange is the dominant ionization process in the outer heliosphere. In addition, ionization of ISN atoms gives birth to pickup ions (PUIs), suprathermal particles that constitute the energetically dominant component of the SW plasma in the outer heliosphere beyond ∼10 au from the Sun (Möbius et al. 1985). On the other hand, fast neutral atoms born in the supersonic SW and slower, but hotter, atoms born in the inner heliosheath (IHS, the region between the termination shock, TS, and the HP) are able to reach the LISM up to distances of ∼1000 au before being ionized. In this way, the LISM is heated and decelerated by the interaction with the neutral SW, which affects the formation of a bow shock (BS), the magnetic and velocity fields, the distribution and the nature of turbulence in that region (e.g., Fraternale et al. 2022), cosmic rays of energy as large as 10 TeV (Zhang & Pogorelov 2016;Zhang et al. 2020), etc. The part of the LISM modified by the presence of the heliosphere, regardless of what physical processes are responsible for this modification and which quantities are affected, is now commonly called the very local interstellar medium (VLISM; see the evolution of this term in Zank 2015 and. In particular, the LISM plasma is affected by the heliosphere in two different ways: (i) by the streams of neutral atoms born in the supersonic SW and in the IHS, and propagating outwards beyond the HP; and (ii) by its direct collision with the HP. The most recent extension of the term VLISM covers the case of subfastmagnetosonic LISM, when the region of the HP-LISM interaction can extend to very large distances. It also accounts for all radial directions, including the polar and heliotail directions. Thus, the region of space occupied by the VLISM plasma depends both on the heliospheric processes and on the properties of the unperturbed LISM. To clarify, the LISM-VLISM boundary is not necessarily sharp and the VLISM properties are neither homogeneous nor isotropic. The pristine LISM is expected to be found a few thousands of astronomical units in the direction of the heliospheric "nose" (upstream) and clearly extends beyond 10,000 au in the heliotail and polar directions (Zhang et al. 2020).
Because of the LISM plasma being mediated by the secondary neutral atoms, a BS can only exist if the flow is superfast magnetosonic in front of it. One cannot determine whether a BS exists or not without performing a numerical simulation, even if the unperturbed LISM is superfast. In the absence of the BS, the LISM plasma is modified not only by secondary neutrals but also by a giant rarefaction that extends into the upstream direction. Relatively strong gradients of all quantities exist between the BS and the HP, therefore it is sometimes useful to introduce an interstellar region called the outer heliosheath (OHS), which is filled by the VLISM. The terminology is not consistent throughout the literature, as discussed in the review paper of Kleimann et al. (2022). In the present simulations, we use"OHS" to indicate the dynamic interstellar region relatively close to the heliosphere, where the plasma temperature exceeds the unperturbed value by an arbitrary threshold set to 3%. The region identified in such a way is also used to identify the secondary ISN atoms, as described in Section 4.4.
Numerical modeling is necessary to understand the global, three-dimensional features of the heliosphere and interpret in situ measurements, performed in the outer heliosphere by Voyager and New Horizons (NH), and remote observations by the Interstellar Boundary Explorer (IBEX; see reviews of Pogorelov et al. 2017a andKleimann et al. 2022). Moreover, more sophisticated models need to be developed as new, higher-resolution data will be provided by upcoming missions, such as the Interstellar Mapping and Acceleration Probe (IMAP; McComas et al. 2018).
MHD-plasma/kinetic-neutrals models couple an MHD description of the ionized components to a kinetic description of neutral atoms. Unlike multifluid models, where neutral fluids that originate in the thermodynamically distinct regions of space are described by the separate systems of Euler gas dynamics equations, kinetic models of the neutral atom transport make it possible to determine the actual distribution function of neutral atoms. The source terms that couple the plasma and the neutrals are collected in this case statistically using a Monte Carlo method (Malama 1991;Baranov & Malama 1993;Lipatov 2002;Izmodenov et al. 2003;Heerikhuisen et al. 2006).
Recently, we have developed the first global, threedimensional, MHD-plasma/kinetic-neutrals model that treats both hydrogen and helium atoms and ions self-consistently . Our modeling framework is the Multi-Scale Fluid-Kinetic Simulation Suite (MS-FLUKSS; Pogorelov et al. 2010Pogorelov et al. , 2014Pogorelov et al. , 2017b, and references therein), an adaptive mesh refinement, hybrid parallelized software (Borovikov et al. 2013). Our kinetic model includes charge-exchange ionization (no elastic scattering), photoionization, and energy-transfer effects due to Coulomb collisions. In this paper, we extend the model by adding electrons as a separate fluid which is assumed to be comoving with the ions.
The properties of electrons in the distant supersonic SW, IHS, and VLISM are mostly unknown due to the lack of in situ observations. Nonetheless, electrons are critically important. They affect the pressure balance at the HP, IHS thickness, ionization of neutral atoms and dust particles, generation of radio emissions associated with shocks in the VLISM, etc. Moreover, the development of such a model is necessary to incorporate the effects of electron impact ionization. The latter physical process is known to be of importance in the IHS. The electron thermal speed in the IHS may be large enough to make the electron impact ionization rate comparable to the rate of ionization due to H+p charge exchange (Scherer et al. 2014).
The distribution of temperature in the IHS is still a highly debated topic. Some studies argue that electrons downstream of the TS are significantly hotter than thermal protons, with pressure similar to that of the proton mixture which includes energetic PUIs. For example, Chalov & Fahr (2013), Fahr et al. (2017), andChalov (2019) suggest that the ratio of electron to thermal proton temperatures, T e /T p,th , can be as large as 10 downstream of the TS. Fahr & Verscharen (2016) have proposed an explanation for this, based on their multi-ion model. However, such heating requires a remarkable acceleration of electrons at the TS or inside the IHS. The results mentioned above were obtained by fitting the Voyager/Low-Energy Charged Particle (LECP) electron observations at energies 20 keV using a single isotropic kappa distribution function. Further theoretical analysis of kappa distributions has been carried out recently by Fahr & Dutta-Roy (2019), Scherer et al. (2019), and Fahr & Heyl (2020). Since electrons below 20 keV were not measured, another possible scenario is that the bulk electrons in the IHS have properties similar to those of thermal protons, and that suprathermal electrons do not contribute significantly to the pressure. This would be more consistent with observational data analyses (e.g., Schwartz et al. 1988), numerical and theoretical results (Gedalin et al. 1995), particle-in-cell (PIC) simulations of collisionless shocks, including the TS (e.g., Matsukiyo 2010;Matsukiyo et al. 2019;Gedalin et al. 2020), and near-Earth shock observations by Cluster or the Magnetospheric Multiscale Mission (e.g., Schwartz et al. 2011;Cohen et al. 2019). All the above do not show any preferential heating of electrons.
In addition, Baranov & Ruderman (2013) emphasized the importance of parallel electron thermal conduction in the heliosheath. Izmodenov et al. (2014) demonstrated that the IHS thickness may be significantly reduced when the heat capacity ratio of the plasma mixture in the IHS is tuned toward isothermal values (γ = 1). With the exception of the fully kinetic model presented by Malama et al. (2006), MHD-plasma/kinetic-neutrals models with separate treatment of electrons have not been developed yet. Multifluid models with electrons do exist (e.g., Glocer et al. 2009;Usmanov et al. 2016), but their treatment of neutral atoms is less accurate than the corresponding kinetic models. All current MHD-kinetic models have to incorporate some assumptions about the electrons' pressure. If PUIs are not treated separately (and therefore the properties of the thermal SW are not available), a possible assumption would be that the electron pressure is a fraction of the pressure of the ion mixture, which includes PUIs. It may seem surprising that this fraction is typically set to 1, which potentially can lead to substantial overestimation of the electron pressure and underestimation of the proton pressure, at least in some regions of space, as will be shown below. In fact, "pickup" electrons are not produced by charge exchange. When PUIs are treated separately, it is conventionally assumed that electrons carry the pressure of thermal protons (e.g., Pogorelov et al. 2016). This seems to be another extreme, since Coulomb collisions or wave-particle interactions should affect the ultimate electron temperature.Using the analytic expressions to describe the energy partition between protons and electrons in the heliosheath, Heerikhuisen et al. (2019) found a noticeably smaller heliosphere when the electron temperature was~0.1 of the ion mixture, with additional SW cooling in the heliotail.
Here we investigate the implications of using different models for electrons treated as an additional fluid. We describe different treatments for electrons and ions. In particular, we discuss the process of filtration of the ISN H and He atoms at the HP, the properties of plasma in different heliospheric regions, and the effect of electrons on the heliospheric structure, e.g., the IHS thickness. We further elaborate on the effects of He atoms in the global heliosphere .
Using the most recent consensus values for the LISM parameters (Zirnstein et al. 2016;Bzowski et al. 2019;Swaczyna et al. 2020), we compare and validate our results with NH observations at distances between 20 and 50 au from the Sun (McComas et al. 2021). As the first step in our analysis, we assume that the SW is stationary and spherically symmetric at the inner boundary.
The paper is organized as follows. The physical model is presented in Section 2. In Section 3, we describe the simulation setup and boundary conditions. The results are presented and discussed in Section 4. Concluding remarks are presented in Section 5.

Physical Model
To address the outlined problem, we further improve our three-dimensional, MHD-plasma/kinetic-neutrals model of the SW-LISM interaction (Borovikov et al. 2008;Pogorelov et al. 2008). The current version includes a kinetic treatment of He atoms and a fluid description of He + ions . The transport of neutral atoms is modeled with a Monte Carlo method (e.g., Heerikhuisen et al. 2006Heerikhuisen et al. , 2008Borovikov et al. 2008), which now has two independent kinetic modules (for H and He atoms).

Description of Plasma
The plasma description is based on the solution of the set of conservation laws for the mixture of thermal protons, pickup protons, electrons, and He + ions, hereinafter referred to as the plasma mixture. We will adopt the term ion mixture to refer to the mixture of thermal protons, PUIs, and He + ions, and the term proton mixture to indicate all protons (thermal and PUIs). Note that PUIs are present in the mixture but they are not treated separately. We solve the governing equations for the plasma mixture (Equations (1)-(3)) in the ideal-MHD conservation-law form in the presence of source terms describing charge exchange between ions and neutral atoms. This ensures the exact MHD description of discontinuities. All ion species are assumed to be comoving. This assumption is reasonable, as far as the motion perpendicular to magnetic field lines is concerned, because of the expected influence of the motional electric field (Parker 1963). The difference in velocities of the constituent charged particles parallel to B is possible, but is assumed to be suppressed via the isotropization of newly created PUIs. In addition, we solve the continuity equation for the He + ion density (r + He ) and the pressure equation for + p He (Equations (4)-(5)). We also solve the pressure equation for electrons (Equation (6)).
The temperature of each plasma component, s, is then calculated using the equation of state, p s = n s k B T s , where s = {p th , PUI, e, He + } is used to denote thermal protons, pickup protons, electrons and He + ions, respectively. Hereinafter, the subscript (p) will be used for the proton mixture. Ion temperature, or their thermal speed, is required to describe charge exchange between ions and neutral atoms. The plasma velocity distribution function is also required. Although both the Maxwellian and kappa distributions are currently implemented in Borovikov et al. (2008) and Heerikhuisen et al. (2008), we assume the Maxwellian distribution function for ions when calculating the source terms. Future simulations with kappa-distributed plasma and time-dependent boundary conditions will allow us to carry out a comparison with observationbased studies on the thermodynamics of the IHS that assume a kappa-distributed SW (e.g., Livadiotis et al. 2011).
Since the LISM is collisional on relatively small scales (e.g., the p-p mean free path for Coulomb collisions is 0.3-4 au for temperatures in the range between 7500 and 3 × 10 4 K; see Baranov & Ruderman 2013;Mostafavi & Zank 2018;, we implemented the source terms in the He + and electron pressure equations that describe the thermal equilibration of protons, He + ions, and electrons by Coulomb collisions (Equation (8)). These terms are computed using the local plasma properties. Our equations read Here ρ s , p s , and T s , are the mass density, thermal pressure, and temperature for each plasma component, u and B are the bulk velocity and magnetic field vectors, respectively, and u and B their magnitudes. The total energy density is E = ρ(ò + u 2 /2 + B 2 /8π), where ò is the specific internal energy (ρò = p/(γ − 1)) and p * = p + B 2 /8π is the total pressure. We assume that the adiabatic index γ is equal to 5/3 for all species. The terms r S s E p ,m, , represent the sources in the mass, energy, momentum, energy, and pressure equations due to charge exchange and photoionization, for each species' s. These term are computed via two Monte Carlo kinetic solvers, whose general description can be found, for example, in Lipatov (2002, see their Ch. 2, Equations (2.175)-(2.180)). S p e for electrons includes photoionization only. The source terms Q s C represent the energy-transfer rate to a species s due to Coulomb collisions with the other species (e-p, e-He + , p-He + ). Expression (9) defines the relaxation rate, n , associated with the process of thermal equilibration. It is valid under the assumption that all charged particles are comoving and have Maxwellian distributions (see Braginskii 1965;Richardson 2019). All quantities in Equation (9) are expressed in CGS units; Λ sj is the Coulomb logarithm (see Equations (10) and (11)), q = 4.803 × 10 −10 statC is the elementary charge, k B the Boltzmann constant (=1.6 × 10 −12 erg eV −1 if temperatures are expressed in eV), and m j and m s are the particle masses. Note that n n =s j js n n j s , where n denotes the number density. Looking at the individual components of the sum in Equation (8), the contribution to the internal energy-transfer rate from species j to species s is equal by magnitude and opposite in sign to that from s to j (note ).

Description of Neutral Atoms
Two independent kinetic modules (for H and He) compute the source terms in the mass, momentum, energy, and pressure equations due to charge exchange and photoionization. Statistics are gathered using a Monte Carlo method instead of solving the six-dimensional Boltzmann equation directly. The particle splitting/recombination strategy adopted here is described in Baranov & Malama (1993) and Heerikhuisen et al. (2009). It is implemented to improve statistics in regions where the plasma grid is fine, while keeping the computational cost relatively low. The current model does not include elastic collisions between ions and neutral atoms with oblique scattering. We model charge-exchange collisions for the processes H + H + → H + + H and He + He + → He + + He. The latter is the dominant process involving He, as shown by Scherer et al. (2014). The source terms S p in the pressure equations are computed from those in the density, momentum, and energy equations through Equation (7) DeStefano & Heerikhuisen 2020).
The following charge-exchange cross sections are used for H+p and He+He + : where  is the collision energy in eV. The expression for hydrogen was provided by Lindsay & Stebbings (2005), while the one for helium is a fit to the Barnett et al. (1990) data in the energy range from 0.16 eV (0.04 eV amu −1 ) to 2 MeV . Note that the He+He + cross sections in Janev et al. (1987) are ∼8 to 16% higher than Barnett's in the energy range we are considering (https://www-amdis.iaea. org/ALADDIN/).

Simulation Setup and Boundary Conditions
In this study we present the simulation results obtained with three physical models. In the first model, Model 0 (simulation S0), we do not include Equation (6) and assume = + + p p p e p H e . Electrons in this model are therefore "hot," their temperature being equal to the temperature of the ion mixture. Coulomb collisions between He + and protons are also included in Model 0. This model was used earlier by . Model 1 (S1) includes electrons as a separate fluid, but the Coulomb collision source term Q C (also for He + ) is ignored. In contrast to S1, Model 2 (S2) additionally includes the collisional heating term. Table 1 shows the boundary conditions that are adopted in all simulations presented here. We assume a spherically symmetric, steady SW. The boundary values at 1 au are based on the averaged OMNI data (in the ecliptic plane) for the period from 2011.5 to 2022. This is done to compare our solutions with NH observations of SW thermal protons and PUIs between 22 and 47 au performed during the same interval of time. The NH PUI density measurements have been recently revised by McComas et al. (2021), which resulted in higher values by a factor of 2 as compared to the previous estimates (McComas et al. 2017). We expect this to pose a number of challenges for modelers. Based on the recent NH data, Swaczyna et al. (2020) obtained a new estimate for the ISN H density of 0.195 ± 0.033 cm −3 , which is ∼26% larger than the previous estimates (e.g., Zirnstein et al. 2016) and closer to the LIC H density suggested by Slavin & Frisch (2008) and more recently recalculated by Bzowski et al. (2019). Here we investigate the effect of such increased ISN H density, while keeping the rest of the LISM parameters the same as in . These parameters are summarized in Table 1 together with the references in the literature. We use the LISM ions parameters, in particular the total nucleon density of ≈0.09 cm   (Zirnstein et al. 2016) bases of IBEX-Lo observations. The uncertainties of the LISM properties are discussed in detail in Bzowski et al. (2019), and cannot be avoided in the present study. For example, the LISM ionization ratio derived from soft X-ray measurements can suffer from their contamination by charge exchange occurring in the heliosphere (e.g., Koutroumpa et al. 2009). The inner boundary of our simulations is placed at 10 au. The solution inside this inner radius is analytical. The proton density at 1 au is 8 cm −3 , exceeding the average OMNI data 4 but within the standard deviation (see the discussion in Section 4). We assume that the electron temperature (T e ) inside the inner sphere evolves nonadiabatically as T e ∝ R −0.85 , similarly to the thermal proton temperature observed by Voyager 2 (V2). An extensive statistical study of the SW properties based on Wind observations spanning the time interval from 1994 to 2005 (Wilson et al. 2018) has shown that the median T e /T p ratio is 0.49 for the fast wind and 1.42 for the slow wind (1.28 for all data), while the ratio of the mean temperatures is 0.96. Previously, Cranmer et al. (2009) presented observational data (from Ulysses and Helios) for the proton and electron heating and discussed the radial evolution of the temperature within 5 au in the fast SW.
The simulations presented in this study were obtained on a Cartesian computational grid with four level of refinement. The simulation domain is a cube with a side of 1680 au. The finest grid size is 0.625 au in all directions. We ran the simulations for at least 1600 yr of physical time. The global coordinate system has the z-axis parallel to the Sun's spin axis. The x-axis belongs to the plane containing the z-axis and V LISM and is directed upstream into the LISM. The Sun's coordinates are (x S , y S , z S ) = (1000, 850, 850), therefore the simulation domain extends 1000 au into the tail direction.
In addition to Equations (1)-(6) we solve the level-set equation to correctly track the position of the HP, as described by Borovikov et al. (2011).

Supersonic Solar Wind
In Figure 1, we show the distributions of pressure (left panels) and temperature (right panels) of each plasma component (from the top to the bottom). The density is shown in Figure 2. We present the quantities in the radial direction intersecting the position of the NH spacecraft at 45 au (ecliptic J2000 longitude and latitude are 286°.6 and 1°.94, respectively). NH observations of thermal SW protons (gray squares), PUIs (green circles), and proton mixture (black circles) are also provided. We can only compare the distributions for the proton mixture (black, continuous curves), since PUIs are not tracked separately in our simulations. The left panels show the total ram pressure, and thermal pressure. The total pressure p+p mag is shown with the cyan line. The positions of the TS, HP, and BS are also indicated. The bulk speed and density agree with observations. To achieve such agreement, we set the density at 1 au to be 8 cm −3 , which is about 23% larger than the average OMNI data. In fact, the proton density measured at NH is found to fall as ∼R −1.83 , which is significantly slower than the R −2 trend for a strictly spherical expansion (Elliott et al. 2019). This behavior is not fully understood and has not been reproduced by global simulations so far. In our simulations, the radial exponent in the total proton density distribution is −1.98, when the fit is computed for R ä [1 au, R TS ]. However, the exponent computed at heliocentric distances R ä [30, 70] au is −1.85, which is caused by the mass-loading effect associated with the birth of PUIs. The exponent further decreases to −1.68 if the fit is done at at a distance of 10 au ahead of the TS. Though timedependent effects, e.g., fast/slow wind stream interactions, may potentially explain the discrepancy (see Richardson et al. 2022), it is interesting to note that reproducing the proton density observed by NH seems to be challenging also for datadriven simulations (e.g., Kim et al. 2018). As far as the radial velocity component is concerned, it decreases by 9.6% at 45 au and by 17% at the TS with respect to its value at 1 au (415 kms −1 ). The decrease of 7.4% is observed between 20 and 45 au. Data-driven models like Kim et al. (2016) are able to reproduce the NH velocity including its low-frequency variations. Our solutions are in a good agreement with the average trend (see the ram pressure in Figure 1(a)-(c)). In the supersonic SW, we do not see significant differences in the properties of the plasma mixture in all three models. The model differences are evident, however, in the pressure and temperature of each plasma components, and in the properties of neutral atoms. When the partial pressure of electrons is assumed to be equal to the pressure of the ion mixture, electrons in the distant SW become hot, since their temperature is equal to that of the ion mixture, which is energetically dominated by PUIs beyond ∼10 au: He This is shown by simulation S0 in the top panels of Figure 1. It is clear that the proton mixture pressure and temperature in S0 disagree with observations (compare the black curve and the black points in Figure 1). For R ä [44, 47] au, the simulation underestimates the proton pressure by a factor ∼2.65. In Model 1 and Model 2, this factor is significantly reduced to 1.35 and 1.38, respectively. As for the proton temperature, Model 0 underestimates it by a factor of 2.3, while the factor reduces to 1.1 for Model 1 and Model 2. For R ä [30, 40] au, the latter two models underestimate the temperature by a factor of 1.45. The origin of this discrepancy may be associated with the rapid increase of dynamic pressure observed at 1 au in 2015, which is a time-dependent effect. During the same period of time, the observed SW speed also exceeds the simulated values.
Our simulations S1 and S2 show that the proton mixture pressure decreases with distance as ∼R −0.83±0.01 (R ä [30, 70] au), while the temperature increases monotonically as ∼R 1.05±0.01 reaching T p ≈ 600,000 K in front of the TS. In S0, instead, T p ≈ 310,000 K. The latter value seem unrealistically low. In fact, assuming that the thermal proton temperature in front of the TS is T p,th ≈ 20,000 K (based on V2/Plasma Science, PLS, observations of Richardson et al. 2008), T PUI ≈ 3 × 10 6 K, and n PUI /n p ≈ 0.25, the proton mixture temperature would be ∼750,000 K, which is more in accord with our simulations S1 and S2. This suggests that the electron pressure and temperature in the supersonic SW cannot not be significantly higher than that of the thermal protons.
In simulation S1, the evolution of electrons from the boundary (at 10 au) is adiabatic and they become extremely cold in front of the TS (magenta curve in Figure 1(d)). Naturally, this is not a realistic behavior. Instead, panels (c) and (f) show the nonadiabatic behavior of electrons in simulation S2 due to the proton-electron energy transfer associated with Coulomb collisions (magenta curves). The electron temperature has a minimum at ∼30 au and increases as ∼R 0.095±0.006 (R ä [30, 70] au). As a result, electrons are heated up to ∼6000 K in front of the TS. This value is probably still lower than expected, since electrons may be heated by the other physical processes, e.g., turbulence/wave-particle interactions, kinetic instabilities, collisionless-shock processes, etc. Nonetheless, it is interesting that the thermal proton temperature observed at NH at its current location at ∼50 au has similar values.
To show a possible contribution of PUI waves to electron heating, we performed an additional run, labeled S3 (Figure 1(f), orange curve), where the right-hand side of Equation (6) includes the heating source term This expression has been commonly used in turbulence transport equations to model the source of turbulence due to Figure 1. Simulated and measured distributions of pressure and temperature along the NH trajectory. From the top to the bottom, the results are shown for simulations S0, S1, and S2. The 25 days averaged NH measurements are shown with the gray squares and green circles for the thermal SW protons and PUIs, respectively. We use black circles to show NH data for the proton mixture (thermal SW and PUIs). The left panels show the magnetic pressure, ram pressure, and thermal pressure of each species. The right panels show the temperature distributions. The bottom panels additionally show the results of simulation S3, which includes a heating source in the supersonic SW due to PUI waves (orange curves). The vertical gray lines indicate the TS, HP, and the BS, which reveals itself as a a weak discontinuity (WD), as discussed in Section 4.3. isotropization of PUIs (e.g., Breech et al. 2009;Isenberg et al. 2010;Usmanov et al. 2016). Here S n PUI is the rate of production of PUIs (discussed later in Section 4.4; see also the left panels of Figure 5), α = 0.4 is the fraction of energy that goes to electrons, and f D = 0.25 is a constant governing the fraction of energy that goes to PUI waves during the isotropization process, which is still the subject of theoretical discussions. In principle, Q turb should affect thermal protons and electrons through the turbulence transport equations which, in turn, provide the heating source terms according to the Kármán-Taylor phenomenology. In the simplified approach presented here (no turbulence transport equations), only the PUI waves generated in situ are considered. It is also assumed that the wave power is fully converted into heat instantaneously. The results of simulation S3 show that the the minimum of T e shifts to 25 au and the radial trend becomes ∼R 0.4±0.02 (R ä [30, 70] au). The combined heating due to Coulomb collisions and PUI waves may increase the electron temperature in front of the TS to ∼10,000 K.

Inner Heliosheath
Immediately downstream of the TS, the proton pressure and temperature are higher in Model 1 and Model 2 as compared to Model 0, similarly to the supersonic SW. The total proton pressure reaches ∼0.12 pPa in S0, and 0.25 in S1 and S2. The latter value is admittedly consistent with the estimates of Dialynas et al. (2020) based on IBEX-Lo observations. It is worth noting, however, that the highest-energy neutrals may not exceed ∼1.1 keV in the simulations presented here. The shock crossing by electrons in our simulations is not adiabatic, because Equation (6) (without sources) ensures the adiabatic behavior only in the smooth flow regions. While that equation is not as accurate for the description of shock transitions as the full set of conservation laws, there is no alternative for us while describing the three-dimensional flow of such complexity. Indeed, the downstream-to-upstream electron temperature ratio (T e,d /T e,u ) in the simulations is ∼25%-30% larger than the adiabatic jump r r = g- . It is known from the theoretical analyses of collisionless shocks and observations that the ratio of the electron temperature downstream of the TS to the adiabatic value depends on the the electron thermal speed and upstream Mach number (Schwartz et al. 1988;Gedalin et al. 1995;Matsukiyo 2010). For upstream electron temperature values ranging from 10,000 to 20,000 K (thermal speeds v th,e = 550-1100 km s −1 ) and Alfvén Mach number M A ∼ 7.5, the abovementioned papers suggest that the electron heating across the shock cannot exceed the adiabatic heating considerably. Assuming a compression ratio at the TS of ∼2.5-3 and an upstream bulk flow speed of 320 km s −1 (Burlaga et al. 2008;Richardson et al. 2008), one obtainsT e,d ad [20,000, 50,000] K for adiabatic heating. The PIC simulations of Matsukiyo & Scholer (2014) show that for the upstream electron beta β e = 0.1 (as in S2) the relative gain of thermal energy acquired by electrons is between 1% and 3%, leading to T e,d /T p,th,d ∼ 0.2−0.5 downstream of the shock. According to V2 observations (Richardson et al. 2008) the thermal proton temperature behind the TS was ∼180,000 K, therefore we estimate that the downstream electron temperature may fall into the range between ∼35,000 and ∼100,000 K. These values are significantly smaller than the 10 6 K obtained from S0 and seem to be in a better agreement with the results of S2 and S3.
The  (2019), as mentioned in Section 1. Based on V2/LECP data, those studies conclude that the electron temperature in the IHS may be 6.7-10 times higher than the thermal proton temperature. This would imply an overadiabatic heating of electrons at the TS by a factor ∼15-50. The physical mechanism that may explain so large an electron heating is still debated, as it disagrees with the studies mentioned in the previous paragraph.
An extremely high peak of magnetic pressure at the HP and the corresponding minimum of thermal pressure are wellknown consequences of simulations that assume a unipolar heliospheric magnetic field (HMF), as discussed in detail by Pogorelov et al. (2015Pogorelov et al. ( , 2017aPogorelov et al. ( , 2021. We still follow this approach for the moment, to increase the complexity of modeling gradually. As far as the proton mixture is concerned, its temperature in the IHS is essentially in accord with the estimates based on Figure 2. Simulated and measured density distribution of different plasma components along the NH trajectory. The cyan curve shows the total nucleons' density. From the top to the bottom, simulations S0 to S1 are shown. The thermal SW proton and PUI density data from NH are shown with gray squares and green circles, respectively. The magenta points show the electron density derived from Voyager 1/Plasma Wave Subsystem plasma wave observations in the VLISM. IBEX data. For instance, Livadiotis et al. (2011) find that the IHS temperature is mostly distributed around ∼10 6 K and has a maximum of ∼2 × 10 6 K. In the present calculations the temperature of the mixture of protons immediately behind the TS ranges from ∼1.5 × 10 6 to ∼3 × 10 6 K and decreases with the radial distance toward the HP, as shown in Figure 3 (black curves). The highest proton temperature (and the largest radial gradient) is found in the models with colder electrons.
In simulations S1 and S2 (featuring cold electrons), we observe a reduction of the IHS thickness, as compared with simulation S0 (see Figure 1). This is associated with a decrease in pressure and temperature of the plasma mixture in the IHS, as shown in Figure 3 (shaded region). This is one of the important results of this study, which will be discussed in more detail in Section 4.4 below. It is in qualitative agreement with the results obtained by Heerikhuisen et al. (2019), who tested the effects of colder electrons on the structure of the heliosphere assuming that the electron to proton temperature ratio in the IHS is T e /T p = 0.1. It is worth noticing that this ratio obtained from our Model 2 is noticeably smaller.The decrease in the IHS thickness seems to be in agreement with Voyager observational data that suggest the width of the IHS along the Voyager 1 (V1) and V2 directions to be ∼30 au. However, it is important to emphasize that the actual IHS thickness is time dependent and was not actually measured by the Voyagers since the spacecraft crossed the TS and the HP at different moments of time. The IHS width in the present simulations is probably still larger than expected, but it is known that a separate treatment of PUIs would lead to a decrease of the IHS width (e.g., Pogorelov et al. 2016). In the present simulations, there is a small difference between the IHS width in the V1 and V2 directions. In the time-dependent, datadriven model of Kim et al. (2017, where PUIs were not treated separately), the IHS width is very variable, ranging from 30 to 45 au, with the width being often ∼5 au larger in the V1 direction as compared to the V2 direction. The latter result holds in the present simulations, as well as in recent energetic-neutral-atoms-based analysis of IBEX data by Reisenfeld et al. (2021), where the results of different models for the IHS thickness are also compared (see Table 1 therein). This topic was recently discussed in Pogorelov (2023) and highlighted in the review paper of Kleimann et al. (2022, see their Section 9.2). The radial stand-off distance of the HP in the V1 direction is ∼125.5 au in S0, and ∼121.8 au in S1 and S2. Along the V2 direction, the distance to the HP is smaller by ∼0.6-1 au. In reality, V1 and V2 crossed the HP at 121.58 au and 119.0 au, respectively. It should be recalled, however, that the HP position is not stationary. V1 crossed the HP in 2012 August, approximately 2 yr before solar maximum, while the V2 crossed it in 2018 November, or ∼2 yr before solar minimum.

Outer Heliosheath
All simulations presented in this study show the presence of a BS in the upstream direction, which is so weak that it practically degenerates into a weak discontinuity (WD). In contrast to MHD shocks, WDs exhibit only jumps in the quantity derivatives but not in the quantities themselves.
With the parameters from Table 1, the fast-magnetosonic Mach number in the unperturbed LISM is M f ∼ 1.1, which represents a necessary but not sufficient condition for the existence of a fast-mode shock in front of the HP. Figure 4 shows the radial Alfvén and fast-magnetosonic Mach numbers (M Ar and M fr , respectively) along a particular radial direction that intersects the BS at 90°. It is seen that M Ar > M fr > 1 far upstream of the shock. However, M fr becomes close to 1 in front of the BS. The quantities ahead and behind of the BS are clearly affected by the secondary ions created by charge exchange with ISNs and with the secondary neutrals born in the SW. The mediation effect is stronger in S1 and S2. It is also of interest to note that when He + ions are not included in the simulation, while all other boundary parameters of the LISM are kept unchanged, M f drops below 1, becoming ∼0.85. In this case the BS disappears completely (see the gray curves in Figure 4 from simulation S0_noHe_nH0195 presented in the Appendix).   In S0 and S2, relatively small mean free paths for Coulomb collisions make the plasma components reach almost a thermalequilibrium state. In these cases, the temperatures in front of the HP along the V1 and V2 trajectories range from ∼32,000 K (T e ) to ∼40,500 K (T p ). In S1, without Coulomb collisions, the proton temperature reaches ∼47,000 K, electrons are only heated to ∼10,000 K due to their compression, and He + ions reach temperature of ∼15,500 K. The proton temperatures in front of the HP in our simulations are in agreement with thermal plasma measurements from the PLS instrument on board V2, which recorded T p ∼ 30,000-50,000 K (Richardson et al. 2019). Unfortunately, large error bars in observations do not allow for a more precise comparison between the observed and modeled temperatures.
It is worth pointing out that inclusion of He + ions and He atoms requires modification of the LISM proton density, to keep the HP at the observed position (e.g., Bzowski et al. 2019). One puzzling aspect concerns the electron density in front of the HP observed by the Voyager/Plasma Wave Subsystem (Gurnett et al. 2020). As shown with the magenta symbols in Figure 2, the simulations presented here are underestimating the electron density. This does not happen in simulations without He when the LISM proton density is set equal to ∼0.09 cm −3 to compensate for the lower dynamic pressure caused by the absence of He + ions; see Figure 12 in the Appendix. This remains an open issue that will be addressed in the future. In particular, it is worth investigating the effects of larger p/He + density ratios in the LISM. It cannot be excluded also that involving solar-cycle-related variability may resolve this puzzle.

Charge-exchange-driven Cooling Mechanism and Shrinking of the Inner Heliosheath
To explain the effect of cooling described in Section 4.2 on the plasma mixture in the IHS, we look at the source terms in Equations (1)-(3) obtained with the Monte Carlo method. The left panels of Figure 5 show the source term S n PUI in the meridional plane (x-z). This term represents the rate of production of new ions per unit volume due to charge exchange or photoionization, and has units of cubic centimeters per second. Although we do not track PUIs as a separate component in this study, we can still estimate their production rate. For the sake of simplicity, we assume that any chargeexchange event creates a PUI, regardless of the region of space where it occurs. Otherwise, different populations of PUIs should be taken into account (Malama et al. 2006), some of them with characteristics closer to the parent ions (see, e.g., Pogorelov et al. 2016).
All three models are compared in Figure 5. The right panels show the pressure source term, S p , due to charge exchange for the the plasma mixture. Since there is no net production of ions (S n = 0 except near the Sun due to photoionization), S p is proportional to the heating/cooling rate of the plasma mixture due to charge exchange.
We find that the charge-exchange rate (i.e., S n PUI ) increases significantly in the IHS when electrons are treated separately (simulations S1 and S2), being enhanced up to ∼50% in the upstream direction. As a consequence, the pressure source term S p is found to be negative throughout the IHS, but its absolute value increases both in S1 and S2 as compared to S0.
To understand the reasons for these changes, we have inspected the expressions for the analytical source terms presented in Pogorelov et al. (2016). A comparison between our simulation results and the analytical model is shown in Figure 6 at a fixed location in the IHS. This figure shows the dimensionless source terms as functions of the proton mixture temperature only. The reference values are =´-S 9.14 10 n ref 9 cm −3 s −1 and = S 0.0583 p ref pPa, respectively. The qualitative behavior can be summarized as follows: 1. Treating electrons as a separate fluid results in colder electrons, because hot electrons are not created by charge exchange. Then, as a direct consequence of the conservation laws, the proton mixture becomes hotter. This is shown in Figure 1 (black curves). 2. Hotter protons result in higher average p-H collision energy and, consequently, higher probability of charge exchange. This increases S n PUI , as shown in Figure 6 (magenta curve and symbols) and in the left column of Figure 5. 3. For plasma conditions typical of the IHS, we find that the pressure source term has a negative sign ( Figure 5, right panels). Since ∂S p /∂T p is always negative, an increase in the proton temperature leads to an enhanced pressure drain from the IHS and associated cooling of the mixture (green curve in Figure 6). 4. Interestingly, there exists a special proton temperature T p,c such that S p > 0 for T p < T p,c . In this case, chargeexchange collisions heat the plasma mixture. This threshold value, while dependent on the other parameters, is about 75,000 K for the conditions used in Figure 6, which is significantly lower than realistic temperatures of the proton mixture in the IHS. In contrast, the condition S p > 0 is satisfied in the supersonic SW and in the OHS ( Figure 5, right panels), where charge exchange heats the plasma mixture. In the OHS, when electrons are colder, as in Model 1 without Coulomb collisions, hotter protons result again in cooling of the OHS plasma mixture. This explains the lower temperature in front of the HP in simulation S1 (compare the black and green curves in Figure 3). In the OHS, the presence of helium and Coulomb collisions must be taken into account to explain the profiles of simulation S2 (magenta curve in Figure 3). 5. Squares and rhombi in Figure 6 show the source terms obtained in our simulations with Model 0 and Model 1 at a fixed location in the IHS (R = 100 au, upstream direction). They are close to the analytical predictions. However, the variations δS n, p /δT p are enhanced in our simulations. This can be understood on realizing that that the analytic expressions are calculated assuming Maxwellian distribution functions for both ions and neutrals. In our kinetic simulations, however, only the plasma mixture is assumed to be Maxwellian, whereas the neutral atoms have anisotropic distributions self-consistently obtained from the Monte Carlo simulations. 6. It is not straightforward to link the abovementioned source term behavior to the actual steady-state values of temperature and pressure shown in Figure 1. Lower bulk flow speeds in the IHS, especially near the HP, leave more time for the source terms to act effectively, which explains why the major changes in the bulk quantities are obtained in the IHS. As a result, a global shrinking of the IHS is observed. This is also seen in the downwind region, which is not steady, as shown in Section 4.5. The Figure 5. Two-dimensional distributions of the PUI density source terms, S n PUI , in the meridional plane (left panels) and of the pressure source terms for plasma, S p (right panels). From the top to the bottom, we show the distributions obtained from models S0, S1, and S2, respectively. The black dashed line shows the position of the HP tracked using the level-set method. IHS thickness reduces by ∼3-5 au in the nose direction and up to ∼60 au in the tail, when measured in the B-V plane, along the direction perpendicular to axis defined by V LISM . This can be seen in both Figure 9, which shows the two-dimensional distributions of the temperature (top panels) and current density (bottom panels) in the B-V plane, and in Figure 10, which highlights some features with three-dimensional plots.
The distributions of H and He atoms are also model dependent. In Figures 7 and 8 we show the densities of different populations of neutral H and He atoms, respectively. Their definitions are provided below. We have derived such distributions using virtual spherical probes placed at chosen points belonging to the simulation domain. This allowed us to obtain good statistics in relatively small regions of space, which is a challenge for this kind of kinetic simulations. To ensure a meaningful comparison with near-Earth observations, one should avoid averaging over excessively large volumes. The radii of such spherical probes are set to 2 au for the leftmost point shown in the figures and increases nonlinearly to the largest radius of 10 au, at 680 au.
Population 0 of neutral atoms (H (0) , He (0) ) includes the ISN atoms that originated in the unperturbed LISM. Those of them that arrived at the inner boundary of the computational region obviously experienced no charge exchange at all, except for the trivial events in the uniform, unperturbed LISM. Since the LISM plasma is perturbed by the presence of the HP and the BS is not universally present, the definition of Population 0 requires some additional criteria to distinguish them from Population 1. Population 1 neutrals (H (1) , He (1) ) include atoms born in the VLISM where the plasma mixture temperature is T > 1.03 T LISM . Note that, with the above threshold, the region where He (1) atoms (the so-called "warm breeze") are created extends upstream of the BS ). Population 2 (H (2) , He (2) ) includes atoms born in the IHS, and Population 3 (H (3) , He (3) ) refers to those born in the supersonic SW. Since the density of He + ions is assumed to be very small at the inner boundary and alpha particles are not included, He (2) and He (3) are negligible and therefore they are not shown in Figure 8. Figure 7 shows the effect of filtration of ISN H atoms. The distribution of different H populations for our three simulations are presented. The left panels show the density profiles. The right panels show the relative differences of the densities obtained in S1 and S2 with respect to those in S0. For hydrogen, n H /n H,LISM = 0.6 at 50 au in S0, and decreases by ∼2.3% in S2. The qualitative behavior of these curves is in agreement with the hydrogen density distributions shown in Pogorelov et al. (2009), where the comparison with the multifluid model results was also provided.
The distributions are compared to the fit of NH observations carried out by Swaczyna et al. (2020). The analytical formula employed in that study is , where λ ≈ 4 au is the H ionization cavity size, θ is the angle between the radial direction toward the probe and the H inflow velocity vector, andñ TS is assumed to be the density at the TS. In Figure 7(a), we present these analytical estimates and a few NH data points, which are shown with the yellow curve and square symbols. The H density at the TS derived by Swaczyna et al. (2020) is 0.127 ± 0.019 cm −3 , which agrees well with our simulations. However, the radial gradient of the H density is larger in our solutions. We obtain a power-law fit: . It incorporates the changes due to variations in the TS and HP stand-off distances from simulations S0-015 and S0-022 presented in the Appendix. As expected, the largest differences are seen in the density of H (2) atoms. The fraction of them flying in the supersonic SW region and reaching 1 au increases up ∼35% (Figures 7(c) and (h)). Despite their low density, H 2 neutrals are significantly hotter and affect the tail of the velocity distribution functions observed at 1 au. As far as the H 1 atoms created in the OHS are concerned, their density in the supersonic SW is similar in S0 and S2. The small differences below 5% are likely due to the variation of the HP position, since Coulomb collisions make the OHS plasma properties similar in these two simulations. An increase of up to 13% is seen in S1. Interestingly, the total density of all all H atoms at 1 au is larger in the models with separate electrons (Figure 7(f)). However, the fraction of pristine ISN H atoms reaching 1 au decreases by 8% and 17% in S1 and S2, respectively.
We observe a small increase ∼2.5% of He (1) atoms reaching 1 au in simulations with separate electrons (Figures 8 (c) and (d)). We were not able to detect significant differences in the total neutral He density.

The Heliotail Configuration
We found the heliotail to be unsteady in our simulations. This was verified by running them for more than 1800 yr of physical time. Similar instabilities were obtained in a similar statement of the problem (unipolar HMF) by Pogorelov et al.  (2015), both in the multifluid and MHD-kinetic statements of the problem. In agreement with the earlier theoretical studies of Yu (1974), two columns of increased plasma density are created in the heliotail. At the current state of knowledge, it is clear that their presence is entirely due to the absence of solarcycle effects, and they are especially prominent when the HMF strength is exaggerated by the assumption of unipolar HMF (Pogorelov et al. 2017a). Also in agreement with Yu (1974) and Pogorelov et al. (2015), these columns become unstable due to some kind of kink instability, which eliminates the reason for the flow collimation and excites instability of the HP. We observe that these instabilities are enhanced when electrons are treated separately. This is well seen in the current density distributions shown in Figures 9(c) and (d) and in the three-dimensional plots shown in Figure 10; see also Figure 1 in Pogorelov et al. (2015). The simulations of Michael et al. (2021), performed in the absence of helium and without electrons treated separately, also show large-scale oscillations in the tail. The onset point of instabilities in the Michael et al. (2021) simulations appears to be significantly closer to the zaxis if compared to our simulations, especially S0. A Rayleigh-Taylor mechanism proposed by Opher et al. (2021) seems to be applicable only to the initial stage of their simulations where the uniform, high-density distribution of H atoms is suddenly introduced into the steady-state MHD solution . The nature of the instability is different in the present study, and its onset occurs where the SW flow vector and the ISN flow vectors are nearly aligned. Thus, the flow pattern is typical of a kink instability.
We suggest that the cooling mechanism described in Section 4.4 plays an important role in the enhancement of this instability. The onset point moves closer to the z-axis in the simulations with separate electrons. It is worth pointing out that the momentum transfer due to charge exchange acts against the plasma flow, being directed toward the positive direction of x. This drag is mostly aligned-and opposite-to the bulk flow, and is stronger inside the lobes where the charge-exchange rate is higher (Figure 5), which may affect the stability of the plasma. The heliotail splitting does not occur in the simulations presented here (at least, not within 1000 au from the Sun). The two-lobe flow collimation persists because of the assumption of unipolar HMF, but clearly the heliosphere is not tailless, as stated by Michael et al. (2021).
Note, however, that the heliotail structure obtained in our simulations does not necessarily represent the one in the realistic heliosphere. The inclusion of the solar cycle, especially the variable boundary between the slow (dense) and fast (rarefied) SW, affects the solution in the heliotail drastically, even under the assumption of unipolar HMF. Interestingly, the solar cycle forces the system at a frequency ∼2.9 × 10 −9 Hz (11 yr), which is similar or even higher than the typical frequency of fluctuations found here. An estimate of the wavenumber along the x-direction yields k x ∼ 2.8 × 10 −13 m −1 (∼150 au wavelength), corresponding to the frequency ∼2.2 × 10 −9 Hz (14 yr period) obtained using an average bulk flow speed of 50 km s −1 . Eventually, including the realistic HMF leads to remarkable changes of the tail's topology. A detailed discussion is provided by Pogorelov et al. (2015Pogorelov et al. ( , 2017a. Nonetheless, the global cooling effect of the IHS plasma described in this study is applicable to timedependent simulations and plays some role in the onset and evolution of instabilities in more realistic heliotail solutions. We also emphasize that an ad hoc strategy to compute the charge-exchange source terms should be implemented to ensure low noise in them and to obtain improved statistics in the presence of unsteady structures and relatively short-scale fluctuations, as discussed in Heerikhuisen et al. (2013). The occurrence of spurious instabilities driven by poorly resolved kinetic source terms has to be carefully evaluated, as well as the effect of grid resolution and domain size. Here we set a short cadence for the kinetic neutral solvers (equal to 3-5 times the time step used for the plasma), to prevent the occurrence of mirror-like spurious fluctuations that can develop in the IHS and near the poles, as well as grow in the tail. A new version of our kinetic code more suitable for the study of time-dependent heliosphere will be presented in a separate publication.

Role of Coulomb Collisions
We have investigated the possible effect of Coulomb collisions on the process of thermal equilibration of all charged particles. The top panels of Figure 11 show the inverse relaxation times (collision frequencies) for p-e collisions (left panels) and p-He + collisions (right panels). We note here that they represent the relaxation times of the thermal-equilibration process (energy exchange) that are used in the pressure equation. Therefore, they are not the basic Coulomb collision frequencies (see, e.g., Zank 2014;. As shown in the previous sections, our simulations demonstrate that Coulomb collisions are mostly effective in the VLISM, where the mean free paths are ∼0.3-10 au for temperatures from ∼7500 to ∼50,000 K. This determines the These distributions demonstrate the charge-exchange-driven cooling effect caused by the inclusion of fluid electrons, and the resulting global shrinking of the heliosheath. The structures resembling detached magnetic vortices are actually two-dimensional cuts of the unstable, collimated plasma columns undergoing a kinklike instability (see Figure 10). quasi-thermalization of electrons, protons, and He + ions. We have determined that the proton and He + temperatures in front of the HP experience significant changes due to Coulomb collisions. In their absence, the plasma components have noticeably different temperatures. This has direct implications for the properties of secondary H and He atoms in the OHS, which are ultimately observed by IBEX at 1 au. In our model, the secondary and primary protons are treated in the OHS as a unique population of thermal protons. Here we label as "secondary" the protons born in the OHS due to charge exchange between the ISN H atoms and thermal protons. Assuming a typical collision energy of 3.4 eV, the initial temperature of the secondary protons would be ∼39,500 K. The mean free paths for Coulomb collisions between these secondary protons and primary thermal protons, the temperature of which is T ∼ 10,000-25,000 K because of the adiabatic compression, is about 5-15 au, which is relatively small. Moreover, the density of primary protons that never experience charge exchange in the OHS is very small as compared to the density of the secondaries. This justifies our treatment of the secondary protons as a part of the thermal population as a whole.
As shown in Figure 11, the energy-transfer mean free paths for p-e collisions is λ ep ∼ 4.5 au in the VLISM adjacent to the HP, whereas it decreases to ∼0.3 au at infinity. For He + -p collisions, the mean free path is l+ 0.8 pHe au at the HP and ∼0.04 au at infinity. In the supersonic SW at 70 au from the Sun, λ ep ∼ 27 au. This is sufficient to heat electrons to ∼6000 K in front of the TS, which is in good agreement with the results of Usmanov et al. . We show four surfaces: the iso-β surface at β = 100, the iso-J surface at J = 6 × 7 × 10 −17 A m -2 , the TS (dark gray), and the HP (translucent gray). The first two surfaces can be distinguished by the different color maps, both showing the temperature of the mixture. A few velocity streamlines that originate in the LISM and reach the vicinity of the stagnation point at the HP are also shown. especially in the presence of Coulomb collisions. In our simulations, Coulomb collisions result in the thermal equilibrium of charged particles in the OHS, making the solutions of S0 and S2 similar in this region.

Concluding Remarks
Building on the results of , we have developed a new MHD-plasma/kinetic-neutrals model of the global heliosphere, which in addition to the self-consistent treatment of hydrogen and helium atoms and ions includes a separate pressure equation for electrons. We have performed three main steady-state simulations (S0, S1, and S2) using the current consensus properties of the unperturbed LISM and 11 yr averaged OMNI data as boundary conditions.
In the supersonic SW, our results show that treating electrons separately allows us to better reproduce the NH proton observations between 20 and 47 au. In particular, we are able to reproduce the observed proton mixture temperature and pressure only if the electrons have a temperature similar to that of thermal protons.
One of the nontrivial consequences of treating electrons as a separate fluid is the reduction of thermal energy of the plasma mixture in the IHS. The physical mechanism underlying this cooling, described in Section 4.4, resides in the increased rate of p+H charge exchange and associated absolute increase of the (negative) pressure source term. Consequently, the IHS Figure 11. Coulomb collision relaxation rates for thermal equilibration (left panels) and the corresponding mean free paths (right panels) are shown in the meridional plane. The top panels show the quantities associated with the p-e collisions, and the bottom panels show the He + -p collisions.
Here we present the results of four additional simulations that we have carried out to investigate the effect of modifying the densities of LISM neutral H, protons, and He + ions. Table 2 shows the values of density for all ion and neutral atom species in the unperturbed LISM that we used in the simulations presented in this paper, both in the Appendix and in the main text. The quantities not shown in Table 2 are identical to those in Table 1. For Model 0, we have additionally run simulations for a hydrogen density of 0.15 cm −3 (simulation S0-015) and 0.22 cm −3 (S0-015). We have also compared simulations with and without helium ions and atoms. In the first case (simulation S0-noHe-0054-0195), we preserve the parameters of S0 but remove helium ions. In particular, the LISM proton density is 0.054 cm −3 . In the second case (simulation S0-noHe-0089-0195), we increase the LISM proton density to 0.899 cm −3 in order to achieve the same ion mass density as in S0, therefore the same LISM ion dynamic pressure. As compared with the simulations involving helium, in S0-noHe-0054-0195 we find a ∼15% decrease of total pressure (thermal+magnetic) in front of the HP, a ∼30% decrease in electron density, and a ∼10% increase in temperature. The HP moves outward by ∼12 au, the TS moves outward by 6 au, and the BS disappears.
In S0-noHe-0089-0195, we find a ∼10% decrease in total pressure, a ∼33% increase in electron density, and a ∼15% decrease in temperature. The HP moves inward by ∼5 au, the TS moves outward by ∼1 au, and the BS moves inward by ∼5 au 8.98 × 10 −3 5.40 × 10 −2 2.20 × 10 −1 1.53 × 10 −2 Model 0, increased LISM H density S0-noHe-0054-0195 0 5.40 × 10 −2 1.95 × 10 −1 0 A s S0, but no helium S0-noHe-0089-0195 0 8.99 × 10 −2 1.95 × 10 −1 0 A s S0, but no helium and increased n p,LISM (same p ram,LISM ) Note. All densities are in units of cubic centimeters. All other parameters as in Table 1. and becomes more pronounced (see the black and red symbols in Figure 12). Figure 12 also shows the simulation results obtained with the models treating electrons as a separate fluid (blue and magenta symbols). We show the results along the V1 trajectory direction and in the direction toward the stagnation point on the LISM side of the HP. The stagnation point in our simulations is in the direction of (237°.4, 11°.3), when expressed in ecliptic J2000 coordinates. There is an uncertainty of ±2°.5 for both angles.
He + ions carry about 40% of the total plasma ram pressure in the LISM. This is the major factoraffecting the global heliosphere, especially because of the uncertainty in the He ionization degree in the LISM. Numerical models are also dependent on this uncertainty. The momentum transfer due to charge-exchange collisions is a second-order effect. It is of interest to see that the HP position moves outward by ∼1 au if the He+He + charge exchange is suppressed in a model where the same He + density is preserved.
In the original publication, the Faraday Equation (1) was inadvertently omitted from the system of equations presented in the paper.
We have also realized that the cyan curves describing the total pressure in panels (b) and (c) of Figure 1 are incorrect due to a compilation error. The corrected panels are provided in Figure 1 below.
In addition, we have corrected a misprint in the legend of Figure 12. The legends of simulations S1_0195 and S2_0195 were mixed up. The blue points actually represent the results from simulation S2_0195, while the magenta points represent the results from simulation S1_0195 (with the exception of panel (c)).
These changes do not affect the simulation results and conclusions.