The Effect of Suprathermal Protons in the Heliosheath on the Global Structure of the Heliosphere and Heliotail

In the interaction between the solar wind (SW) and the local interstellar medium, various processes create ions with energies up to ∼10 keV that are out of thermal equilibrium with the “core” population. Wave–particle interactions tend to isotropize the velocity distributions, but the collisionless nature of the SW precludes thermalization. Suprathermal protons can charge-exchange with interstellar hydrogen, producing energetic neutral atoms that are seen by the Interstellar Boundary EXplorer spacecraft. We have developed a model for the presence of several suprathermal populations in the SW downstream of the heliospheric termination shock. The model uses magnetohydrodynamics to satisfy the first three moments of the total ion distribution, and couples these through charge-exchange to neutral hydrogen, conserving mass, momentum, and energy in the combined system. The proton population is separated into a cool core and three suprathermal populations, and hydrogen atoms may charge-exchange with protons from any of those four populations. The phase-space properties of the pick-up ions are selected based on data and theoretical considerations. In this paper we quantify the impact of suprathermal protons on the global structure of the heliosphere by comparing our new model to a traditional Maxwellian fluid model, and a kappa-distribution model. We find that the differences in momentum and energy transfer rates from the protons onto neutral hydrogen between the models leads to different plasma properties in the heliotail, and also changes the size of the heliosphere. Including the energy-dependent charge-exchange cross section into the collision integrals reduces the magnitude of these differences.


Introduction
The heliosphere is created from the collision of two plasma flows-the solar wind (SW) that emanates from the Sun and the local interstellar medium (LISM) that permeates the local galactic environment surrounding our solar system. The need for improvements in the fidelity of models of the SW-LISM interaction has grown in step with the increasing availability of data from the boundaries of the heliosphere. In particular, NASA's Voyager 1 (V1) and Voyager 2 (V2) probes have sent back invaluable in situ data (e.g., Burlaga et al. 2005Burlaga et al. , 2008Stone et al. 2008Stone et al. , 2013, including encounters with the heliospheric termination shock (TS, in 2004 by V1 and in 2008 by V2) and heliopause (HP, in 2012 by V1 and in 2018 by V2). Remote sensing data of the plasma environment in the heliospheric interface comes mostly from energetic neutral atoms (ENAs) that are born in charge-exchange collisions between energetic protons and the cold slow-moving interstellar hydrogen atoms that pass through the SW-LISM region. NASA's Interstellar Boundary EXplorer (IBEX; e.g., McComas et al. 2009McComas et al. , 2017 spacecraft collects such ENAs from Earth orbit and produces all-sky maps of ENA flux every six months. This combination of in situ and remote sensing data has been the primary driver for model development over the past 15 yr. It has long been recognized that interstellar neutrals play a critical role in shaping the heliosphere (Wallis 1975). Early simulations often utilized axially symmetric geometry without magnetic fields, while modeling the neutrals as either particles (Baranov & Malama 1993) or, more computationally expedient, as a fluid (Pauls et al. 1995;Liewer et al. 1996;Zank et al. 1996b;Fahr et al. 2000). Once the asymmetry of the heliosphere was recognized, when V2 encountered the TS in the southern hemisphere much closer than V1 had done in the northern hemisphere, three-dimensional (3D) simulations that included the LISM magnetic field became the standard (Izmodenov et al. 2005;Opher et al. 2006;Pogorelov et al. 2006;Washimi et al. 2007;Ratkiewicz & Grygorczuk 2008;Heerikhuisen et al. 2008).
Charge-exchange between interstellar hydrogen and SW protons creates ENAs, some of which can be seen by IBEX, but also injects new protons into the plasma that are out of thermal equilibrium. Due to the way in which they are picked up by the moving plasma, such suprathermal ions are referred to as pick-up ions (PUIs). Most heliospheric simulations employ magnetohydrodynamics (MHD) to model the ions, and the way in which PUIs enter the energy balance is through the chargeexchange source terms that are applied to the MHD equations as a result of ion-neutral coupling. This explicitly assumes a complete equilibrium of PUIs, protons, and electrons. While MHD captures the energetics of PUIs, it does so only through the conservation laws of the moments of the proton distribution function. Going beyond this simplistic approach for PUIs in global simulations of the heliosphere has been a continuing goal of the modeling community. An early attempt to model PUIs as separate from, but coupled to, the supersonic SW was developed by Isenberg (1986). Zank et al. (2014) used these ideas to develop a more general PUI mediated plasma model that can be applied to the supersonic SW, the inner heliosheath, and the very local interstellar medium (VLISM). Indeed, the region where the plasma outside the heliosphere is mediated by PUIs can be used to define the VLISM (Zank 2015). Recently, Zank et al. (2018) extended their model to include the transport of low frequency MHD turbulence and the associated heating of the thermal plasma. Comparison of the New Horizons SWAP thermal plasma and PUI data and the Voyager 2 magnetometer and plasma data with the model reveals good agreement. In global simulations of axially symmetric heliospheres and without velocity diffusion, Fahr et al. (2000) included PUIs using a fluid approach, while Malama et al. (2006) solved for the PUI distribution function using a Monte Carlo particle approach. The obvious next step is to include PUIs in 3D simulations of the SW-LISM interaction.
The importance of suprathermal tails in the proton velocity distribution generated by PUIs was highlighted by the first IBEX maps, which showed significant flux of ENAs in the 0.5-6 keV energy range. Taking 3D MHD-neutral simulations of the heliosphere computed without a separate PUI component, Prested et al. (2008) demonstrated that realistic levels of ENA flux could be generated using a post-processing of the solution assuming the protons in the subsonic SW follow a generalized Lorentzian (or "kappa" distribution, e.g., Livadiotis & McComas 2013). Heerikhuisen et al. (2008) then assumed a kappa-distribution for heliosheath protons during the chargeexchange process, allowing the suprathermal tails to affect the source terms. The feedback of charge-exchange with suprathermal protons decreased the width of the heliosheath by several astronomical units, an effect also noticed by Malama et al. (2006) in an axially symmetric heliospheric simulation with PUIs.
In an attempt to include the physics of shock dissipation into determining how suprathermal protons are added to the simulations, Zank et al. (2010) proposed a model with three populations of protons with different temperatures immediately downstream of the TS. The Zank et al. model recognizes that thermal protons and PUIs are transmitted differently through the quasi-perpendicular TS. By virtue of their filled-shell distribution (Vasyliunas & Siscoe 1976) in phase space, a subset of the PUIs are reflected preferentially by the crossshock electrostatic potential, which forms the primary dissipative mechanism at the TS (Zank et al. 1996a). Consequently, the reflected PUIs are strongly heated at the TS, while the remaining transmitted PUIs and thermal protons are only heated adiabatically. Furthermore, some PUIs can even experience multiple reflections (Lee et al. 1996;Zank et al. 1996a) to gain even more energy than singly reflected PUIs (see also Lembège & Yang 2016;Kumar et al. 2018;Lembège & Yang 2018, for recent particle in cell simulations that analyze the V2 crossing of the TS). Recent observations of preferential heating of PUIs at an interplanetary shock about 34 au from the Sun by the SWAP instrument on board New Horizons  further validates this approach. The ENA flux resulting from an implementation of the approach by Zank et al. (2010) was tested in Desai et al. (2012), who found agreement to within 50% of IBEX data. In an improved version of the model that also included ENAs from PUIs outside the heliopause the flux levels almost exactly matched the spectral data from IBEX (Desai et al. 2014;Zirnstein et al. 2014) along a particular viewing direction.
In this paper we use the PUI model for the subsonic SW from Zirnstein et al. (2014), itself an extension of the model from Zank et al. (2010), to construct PUI distributions at all points downstream of the TS. We then use these distributions in the computation of the charge-exchange source terms in the neutral particle code, which are subsequently applied to the MHD solver to obtain global heliospheric solutions that include the feedback of charge-exchange with PUIs onto the plasma. In Section 2 we describe the multipopulation PUI model and how it is implemented. The results are presented in Section 3, where we compare four simulations, each with a different way to include charge-exchange between suprathermal protons and neutral hydrogen. A discussion of the results and our conclusion follow in Section 4.

Model of the SW-LISM Interaction
Our model for the 3D steady-state SW-LISM interaction couples a Monte Carlo particle code for neutral hydrogen (Heerikhuisen et al. 2006(Heerikhuisen et al. , 2008 to an MHD code for the ions (Pogorelov et al. 2006(Pogorelov et al. , 2009a and forms part of the MSFLUKSS code (Pogorelov et al. 2011). The particle code tracks H-atoms of interstellar origin through the domain on ballistic trajectories that include stochastic charge-exchange collisions with protons. Each collision results in the H-atom losing its electron and becoming a proton that joins the plasma, while the electron joins with the proton to create a new H-atom with the properties of the ambient plasma. The changes in momentum and energy associated with each charge-exchange collision are stored on a grid and passed to the MHD code to be used as source terms for the next iteration. The ion and neutral codes are iterated until a steady-state heliosphere is achieved, which typically requires about 40 iterations. As we show in Section 3.3, the rate of charge-exchange depends on the relative velocities of potential partners, so that the shape of the proton velocity distribution affects the overall solution we obtain.

Boundary Conditions
We consider four runs of our global heliosphere code, that differ only in the way that neutral hydrogen couples to the protons from the plasma (see Section 2.2). The boundary conditions we apply are the same as the 3 μG case of Zirnstein et al. (2016). On the LISM boundary, we assume an ion density of 0.09 cm −3 , a neutral hydrogen density of 0.154 cm −3 , both at a temperature of 7500 K. The flow speed for both species is 25.4 km s −1 from ecliptic J2000 coordinates (255°.7, 5°.1), and the magnetic field with strength 3 μG is directed toward ecliptic J2000 coordinates (226°. 99, 34°.82). On the inner boundary of our domain (a sphere of radius 10 au), we inject SW with parameters that correspond to the following 1 au values convected by adiabatic expansion to 10 au: 5.74 cm −3 density; 51,100 K temperature; 450 km s −1 velocity; 37.5 μG radial component of magnetic field. The density, velocity, and temperature are uniform over the boundary, while the magnetic field at the boundary follows the Parker spiral. To avoid the presence of an artificial flat current sheet while maintaining the presence of magnetic field in the heliosphere, we adopt a unipolar field for the SW with magnetic field pointing away from the Sun in both the northern and southern hemispheres.
Our simulation naturally creates a "neutral SW" as chargeexchange in the supersonic SW gives rise to a steady output of new hydrogen atoms that move out radially at the SW speed. To generate a realistic spread of neutral SW speeds in the simulation, we make use of historic SW speeds as a function of latitude. The charge-exchange rate for H-atoms traveling through the supersonic SW is computed using the SW speed and density supplied by the MHD code, but when a chargeexchange occurs in the supersonic SW, we set the new hydrogen atom to have a speed that is selected from a 22 yr database (based on interplanetary scintillation observations) of yearly averaged SW speed (Sokół et al. 2013), corresponding to its latitude and a random point in time. Using the MHD density and speed for the charge-exchange rate helps keep the source terms consistent with the system of equations we are solving.

The Multicomponent Suprathermal Ion Model
The multicomponent ion model is based on the decomposition of the IHS proton distribution as described in Zirnstein et al. (2014), itself an extension of the three population ion model presented in Zank et al. (2010). We assume that there are four different types of protons that make up the SW downstream of the TS: (1) the "core" (i.e., thermal) population that came from the Sun and was mildly heated (less than the equivalent adiabatic compression) by crossing the TS; (2) a population of PUIs created in the supersonic SW that were "transmitted" through the TS with moderate heating (approximately adiabatic); (3) a population of PUIs energized at the TS by shock "reflection"; and (4) a population of PUIs "injected" into the IHS plasma when one of the first three populations charge-exchanges with an H-atom.
Assuming the populations are comoving, and that the total velocity distribution is isotropic, the conservation laws for the first three moments (number density, momentum, and energy) of the ion distribution are simply the ideal MHD equations. We ignore the consequences of moments 4 and above (the third moment is zero by isotropy), and instead assume that the four populations can each be approximated by a Maxwellian. Our earlier work without self-consistent feedback showed this to be a reasonable approximation Desai et al. 2012Desai et al. , 2014, consistent with other simulations of PUIs at the TS by Wu et al. (2010). Under this approximation, the total velocity distribution function evolves due to the different charge-exchange rates experienced by each population (see below where n i is the number density and temperature of population i. This leaves us with two constraints and eight unknowns (the density and temperature for each of the populations). In Section 2.3 we explain how we solve for these at each location in the subsonic SW downstream of the TS. The MHD solver mentioned in Section 2 is used to update the total ion distribution, and the presence of the four separate proton populations only affects the coupling to hydrogen. For the multicomponent model, charge-exchange for a given H-atom can occur with any of the four populations. Chargeexchange results in a transfer of energy and momentum between the plasma and neutral populations that is included into the MHD system of equations through source terms. These source terms can be computed by integrating over all chargeexchange partners at a given location. Equation (3) illustrates the energy source term, but the integral is similar for the momentum source.
Here the charge-exchange cross section depends on the collision energy according to (Lindsay & Stebbings 2005) where E is the kinetic energy (in keV) in the frame of the collision and σ ex is in units of 10 −16 cm 2 . For fluid neutrals that follow a Maxwellian velocity distribution coupled to a Maxwellian proton distribution, the source term integrals can be analytically computed provided that the cross section is assumed to be constant (e.g., Pauls et al. 1995;McNutt et al. 1998;DeStefano & Heerikhuisen 2017). For the particle code, we need the charge-exchange rate β for a single hydrogen atom moving through a proton distribution, which is derived by integrating over all possible proton partners: If the cross section is assumed to be constant over the range of interaction energies (i.e., combined bulk motion and thermal spread) between the H-atom and ambient protons, v rel can be computed analytically for either Maxwellian (e.g., Ripken & Fahr 1983;Malama 1991;Heerikhuisen et al. 2006) or kappadistributed protons (Heerikhuisen et al. 2008). However, Equation (4) shows that the cross section becomes small at large energies, such that v rel σ ex (v rel ) decreases above ∼10 keV, reducing the charge-exchange rate with high energy protons. As such, the analytic forms are not valid for very hot Maxwellians (T>20 MK), or kappa-distributions with significant tails (e.g., κ<2). For this reason the charge-exchange rate (5) for H-atoms in our Monte Carlo particle code is computed numerically, and the resulting charge-exchange partner is selected from a numerically precomputed cumulative distribution that includes the effect of the energy-dependent cross section. See Heerikhuisen et al. (2015) and DeStefano & Heerikhuisen (2017) for a more detailed analysis of the role of the cross section in the charge-exchange integrals.

Solving for the PUI Populations
Equations (1) and (2) form the basis of our multicomponent PUI model for the inner heliosheath plasma (see also Zirnstein et al. 2017). To find the four density and four temperature values needed to compute the charge-exchange rate in the H-atom solver at a given location, we first trace the flow back to the TS. At the TS n inj =0, and we need only consider the remaining six unknown quantities-the MHD code providing n MHD,TS and T MHD,TS . The total PUI density just upstream of the TS can be computed by integrating up the production of PUIs along a radial line in the supersonic SW, according to where δ is sufficiently small (∼0.2 au) so that n H is effectively zero initially, dr=u SW dt, and ν ph and ν e are the photoionization and electron impact ionization rates for neutral hydrogen that scale as r −2 with distance from the Sun. The PUI density just downstream of the TS can then be inferred from the relative PUI density just upstream of the TS by using the compression ratio obtained from the MHD solver, so that the PUI fraction is continuous across the shock.
To prescribe the total proton distribution at the TS, we utilize two parameters: (1) the ratio of PUI density to the total density α=n PUI,TS /n MHD,TS ; and (2) the ratio of PUIs that are significantly energized at the TS β=n refl,TS /n PUI,TS , where n PUI,TS =n trans,TS +n refl,TS (since = n 0 inj,TS ). These parameters can be adjusted in response to comparisons between simulation results and observational data, or other theoretical considerations. For the simulations in this paper we fix these ratios to be the same at all points on the TS. Hence we have the three densities just downstream of the TS from a a b ab In the simulations considered here we obtain α from the integral for n PUI,TS using Equation (6), combined with n MHD,TS from the MHD code. In addition we fix β=0.1 everywhere on the TS. This value is based on previous calculations of PUI acceleration at the TS , and the PUI energy implied by the measurements of the TS crossing by Voyager 2 (Richardson et al. 2008). The temperatures at the TS for the three proton populations are obtained by satisfying TS . The three constants in (8) are somewhat arbitrary, but have been chosen to provide the best match to Voyager 2 data at the TS (Richardson et al. 2008), and data of ENA flux measured by IBEX (Desai et al. 2012(Desai et al. , 2014Zirnstein et al. 2014). With these choices of constants, and with α=0.2 and β=0.1, we have * = T T TS MHD,TS , and for a typical MHD (i.e., total) temperature of T MHD =1.2 MK near the TS, we obtain T core =120,000 K, in approximate agreement with measurement of the core SW temperature by Voyager 2 just downstream of the TS (Richardson et al. 2008). Note that in the simulation 0.2α0.25, depending on the location on the TS, so that in general * ¹ T T TS MHD,TS . The density and temperature of each of the three components at all locations just downstream of the TS is computed first. Then for each grid location in the inner heliosheath the plasma streamline is traced back to the TS. Along the way we accumulate the losses due to charge-exchange for each population according to where the loss rate along each streamline segment = s u d dt depends on the relative speed given by (Pauls et al. 1995 From (10) it can be seen that a more energetic (i.e., hotter) population has a higher relative speed for the interaction. Below ∼10 keV, σ ex v rel increases with energy, so that from Equation (9) hotter populations suffer greater losses as plasma propagates through the heliosheath than colder populations (see also Zirnstein et al. 2014). This effect causes a change in the relative abundance of each population as the plasma moves away from the TS, with the more energetic populations being removed more quickly. Note that the loss for each species can be computed in this way because the temperature fractions from Equation (8)  where * T TS depends on which location on the TS the streamline connects to. Note that for particular values of α, β, and the constants in Equation (8) chosen here, * T TS is equal to T MHD,TS . However, this will not be true in general, and hence we retain * T TS in the above description of the method.
All charge-exchange collisions in the inner heliosheath result in the loss of one of the populations present at the TS, with the corresponding gain of a newly "injected" population. Hence from Equation (1 Finally, the temperature of the injected population at an arbitrary location in the heliosheath can be obtained from (2) = ---

Simulations with Four Different Pick-up Ion Models
To understand the impact of the multi-ion model on the global structure of the heliosphere, we ran four versions of our coupled MHD-plasma kinetic-neutral code (see Section 2).
Each version of the code assumes something different about the nature of the PUIs that exist in the inner heliosheath plasma. The assumptions for each of the four cases considered can be summarized as follows: case I The PUIs in the subsonic SW are assumed to thermalize with the core SW, resulting in a proton distribution function that is a single Maxwellian. case II The PUIs form a power-law tail with a cool core. In this case we assume the proton distribution function follows a kappa-distribution, and for the runs considered here we choose κ=1.63 everywhere. case III The proton distribution function downstream of the TS is assumed to be comprised of three distinct populations of PUIs, plus the core SW, as described in Sections 2.2 and 2.3. case IV Here we take case III, but additionally allow the electron pressure to be less than the total ion pressure, in places where the latter is dominated by energetic PUIs. This case is explained in more detail in Section 3.2.
These different forms of the proton distribution enter the calculation only during the charge-exchange process. The proton distribution is assumed to be isotropic, and charge-exchange is the only collisional process where the shape of the proton distribution becomes important. The charge-exchange integrals themselves are solved in a Monte Carlo way by gathering the rates of energy and momentum removed from (or deposited into) the plasma during charge-exchange collisions computed in the particle solver for neutral hydrogen . So even though we solve the same MHD equations in all four cases, the reason the solutions to the SW-LISM boundary value problem are different, is that the energy and momentum source terms due to charge-exchange between protons and hydrogen atoms are different in each case. Figure 1 illustrates how the multicomponent model (cases III and IV) evolves with distance from the TS, compared to the single Maxwellian and single kappa models that maintain their shape in velocity space. This evolution occurs because energetic protons charge-exchange at a higher rate than low energy protons from the core, as explained in Section 2.2. The figure shows that close to the TS, the multicomponent model has significantly more protons with velocities between 2 and 10 thermal speeds than either the single Maxwellian or single kappa-distribution models, while deeper into the heliosheath most of these high energy protons have been replaced by a lower energy "injected" population.

Modeling the Electron Pressure
In MHD, the total pressure P is given by P=P p +P e =2P MHD , where P p is the proton pressure, P e is the electron pressure, and = P n k T MHD MHD B MHD is the pressure obtained from the MHD model, assuming a proton-electron plasma. In most implementations, the number density, velocity, pressure, and temperature of electrons are assumed to be equal to the corresponding values for protons, so that P p = P e =P MHD . Cases I-III described in Section 3.1 follow this convention. However, because the PUIs are injected through charge-exchange, rather than being created via some continuous Liouville-type process, the electrons should not be expected to have the same suprathermal components. So for case IV we set the electron pressure immediately downstream of the TS to be equal to a small fraction of the total (electron + proton) pressure, and then allow this fraction to grow to become equal to the proton pressure at some distance in the heliotail. In other words, we allow > > P P P p M H D e downstream of the TS, while maintaining the MHD assumption = + = P P P P 2 p e M H D . While different values could be justified on the basis of different physical motivations, for case IV we use this approach to set the electron temperature at the TS to 10% of the effective temperature T MHD obtained from the MHD code, and we increase the temperature linearly along a streamline until T e =T MHD at a distance of 800 au. Hence Equation (2) becomes g g min , 1 . 14 core core trans trans refl refl inj inj

MHD MHD
Here l is the distance along the streamline from the TS, L is the distance beyond which we assume an equipartition between proton and electron temperatures, and γ is the fraction of the MHD temperature the electrons have at the TS. For case IV we use L=800 au and set γ=0.1.

Differences in the Energy and Momentum Transfer Rates
One of the key questions we seek to answer in this paper is how the global MHD solution responds to the differences in the shape of the velocity distributions for the four cases we consider. As described in Section 2, the shape of the proton velocity distribution impacts the charge-exchange process, and exerts its effect on the MHD solution through source terms in the energy and momentum equations. The main regions where the source terms impact the solution are: (i) in the supersonic SW, mostly resulting in a slowing and heating due to the injection of suprathermal PUIs; (ii) in the subsonic SW, where ENAs are created by the removal of energetic ions, resulting in cooling of the plasma; and (iii) in the region ahead of the heliosphere between the heliopause and bow shock (or "bow wave," McComas et al. 2012;Zank et al. 2013) where interstellar neutrals and ENAs created inside the heliosphere both deposit energy and momentum. The source terms vary in magnitude by many orders of magnitude, depending on location, so to aid comparisons we focus on the ratio between sources obtained from the different cases we consider.
In Figures 2 and 3 we show the ratio of the source terms in the models where PUIs alter the velocity distribution, and the case where the protons and PUIs thermalize into a single Maxwellian (i.e., case I). The energy source term ratios plotted in Figure 2 show that the supersonic SW is effectively unchanged, as expected given that our PUI models only operate in the subsonic SW downstream of the TS. The ratios also show that the differences outside the heliopause are relatively small, except perhaps case IV which creates significantly more heliospheric ENAs than the other cases. The subsonic SW shows the greatest differences between the source terms due to the differences in how we model the PUIs in that region.
From Figure 1 we see that the kappa-distribution has "wings" of high energy particles and an elevated narrow "core," at the expense of particles with energies between 1 and 3 times the thermal speed when compared to a Maxwellian Figure 2. Ecliptic plane slices of the energy source term from cases II-IV, divided by the energy source term from charge-exchange computed assuming the protons and PUIs thermalize into a single Maxwellian distribution (case I). Note that the actual magnitude of the energy source term varies by 4-5 orders of magnitude in regions where it is significant-i.e., inside, and immediately in front of, the heliosphere. distribution with the same density and pressure. One might expect the high energy wings to yield enhanced energy transfer, but even though such particles are highly mobile, their chargeexchange rate is low due to the energy-dependent cross section (see Equation (4)). As noted in Section 2.2, because our Monte Carlo H-atom code includes the cross section inside the collision integral, the rate of charge-exchange with protons declines precipitously above a few tens of keV. The result is that close to the TS, where the effective plasma temperature is greatest, the kappa-distribution has a lower rate of energy transfer to the neutral hydrogen population than an equivalent Maxwellian. Farther down the heliotail, the wings of the distribution decrease in energy as the plasma cools, to the point where the cross section no longer limits the energy transfer, resulting in a greater energy transfer rate compared to a Maxwellian. Figure 2(B) shows that in the multicomponent model (case III), the subsonic SW loses energy more rapidly than the Maxwellian case within about 200 au of the TS, while farther down the heliotail the loss rate is similar to the Maxwellian case. The main reason is that this model, in contrast to the kappa-distribution model, has most of its high energy distribution below 10 keV, above which the cross section limits the charge-exchange rate. More specifically, from Figure 1 we can see that the multicomponent model has significantly more particles with 8-10 times the thermal speed compared to the kappa-distribution, and close to the TS where the total proton temperature is just over 1 MK this corresponds to protons in the 6-10 keV range.
As expected, Case IV has significantly stronger source terms in the heliosheath, compared to the other three cases. As explained in Section 3.2, in this case the electron temperature close to the TS is low, allowing for a PUI pressure that is enhanced compared to the other three cases in which the electron and proton pressures are assumed equal. Figure 2(C) shows that case IV is about twice as efficient at removing energy from the plasma as case I (note the larger range on the color scale than in Figures 2(A) and (B)). This has a dramatic effect on the size of the heliosphere, so that the heliopause for case IV is approximately the inner edge of the saturated blue region (see also Figure 5) where the source term ratio is in fact negative because the energy sources change sign at the heliopause. Figure 3 shows the ratio of the momentum source term magnitude to the magnitude of the momentum source term for a single Maxwellian. While the MHD code works with the momentum source term vector, the magnitude used in the Figure 3 can be thought of as indicating where the plasma is slowed down by charge-exchange. For the single kappadistribution (Figure 3(A)), the source term ratio is lower close to the TS, indicating less momentum transfer away from the plasma. The reasoning is similar to that discussed above for the energy transfer, and is a consequence of the energy-dependent cross section reducing the charge-exchange rate with the high energy wings of the distribution. Interestingly, as the plasma flows toward the heliotail, it eventually cools sufficiently for the wings of the distribution to participate in the chargeexchange process. At this point we can see an abrupt increase in the momentum source term that corresponds to a significant slowing of the plasma in the tail (this slowing of the heliotail flow is discussed in more detail in Section 3.4). In the more distant heliotail the source terms are similar (see also Figure 4), indicating that the flow speeds are similar in both the Maxwellian and kappa-distribution cases. Figure 3(B) shows that momentum source for the multicomponent model (case III) is slightly lower than the Maxwellian case close to the TS, but similar elsewhere. Figure 3(C) shows that in the region where the energy source for case IV greatly exceeds the corresponding source in the Maxwellian case (Figure 2(C)), the momentum sources for case I and case IV are actually quite similar. The saturated yellow region farther down the heliotail indicates an enhanced momentum transfer rate for case IV, which turns out (see Section 3.4 and Figure 4) to be driven mostly by the higher plasma density (and hence enhanced charge-exchange rate) that develops farther down the heliotail in this case.

Plasma and Neutral Properties in the Heliotail
In Figure 4 we show the bulk properties of the protons and the hydrogen atoms along an axis parallel to the interstellar flow for each of the four cases considered. The single Maxwellian case (blue) acts as a baseline for the comparison with the other cases that model the presence of PUIs in various ways. The most striking features of this figure are the similarity between the plasma properties of case III and case I, and that case II and case IV are generally on opposite sides of the baseline (case I).
The plot of plasma and neutral densities in Figure 4(A) shows a region of almost uniform proton density ∼0.002 cm −3 for about 300 au downstream of the TS. Beyond this region, all cases show an approximately linear increase in proton density to the outer boundary of the simulation domain at 1000 au from the Sun. The neutral hydrogen ionization cavity extends to about the TS, beyond which the hydrogen density remains approximately constant at ∼0.05 cm −3 at least to ∼1000 au into the heliotail. While case I and case III produce relatively similar proton density profiles, case II is noticeably lower farther down the tail, while the density in case IV rises more rapidly to become about double that of the other cases in the distant heliotail.  Figure 4(B) shows the plasma and neutral flow speeds in the heliotail. All cases include a two-step slow-down of the plasma, first a sudden slowing to about 60-70 km s −1 about 100 au past the TS, then a slight acceleration over the next 100 au, followed by another deceleration to half the speed over about 100 au. Beyond this, from about 400 to 500 au from the Sun, the plasma flow speed remains relatively steady around 30-40 km s −1 . In general, case I and case III have similar profiles, while case II slows down later and settles into a slightly higher speed, and case IV slows earlier and is slowest in the distant heliotail. The neutral flow accelerates in all cases from about 20 km s −1 at the TS to about 30 km s −1 near the outer boundary. The neutral flow near the TS in the tail direction is mostly made up of secondary interstellar atoms (i.e., those created in the disturbed LISM around the HP) rather than primary atoms (i.e., those from the pristine LISM), because the latter are colder and thus have straighter trajectories that intersect the ionization cavity around the Sun. This results in a bulk neutral flow that is somewhat slower than the LISM value, because the secondary interstellar atoms are slower and warmer (∼20,000 K), and more easily reached this location without crossing the ionization cavity around the Sun. As charge-exchange proceeds in the heliotail, slow neutrals are replaced with new ones that move with the plasma flow on average, which accelerates the neutral flow to become closer to the plasma speed.
The plasma and neutral temperatures in the heliotail are shown in Figure 4(C). It is clear that even though all models start with almost the same plasma temperature at the TS, the rate at which the plasma cools varies significantly. Note that the proton and electron temperatures for case IV are different, in accordance with the assumptions of the model, as described in Section 3.2. As with the density and speed, the plasma profiles of case I and case III track each other closely, while case II and case IV cool slower and faster respectively. Some of the cooling is due to a transfer of energy from the plasma to the neutrals, so that case II neutrals are coolest, while case IV neutrals are hottest. However, for cases I and III this trend is reversed, likely due to differences in the outer heliosheath properties and the way in which incoming neutrals are filtered. The slightly higher speed of case 2 neutrals near the TS provides additional support for this idea.

Shape and Size of the Heliosphere
Differences among the four cases in the energy transfer rate from protons to neutral hydrogen, due to the different assumptions about the shape of the proton velocity distribution function, results in four different heliospheres. Figure 5 shows the approximate location of the heliopause in the ecliptic plane. The location of the heliopauses shown in the figure are determined using a temperature condition that varies linearly from T;200,000 K near the nose to T;30,000 K where the heliotail exits the domain. This approach ensures that the HP is located in the steepest part of the temperature gradient associated with it. Despite differences in the coupling source terms (Section 3.3) that result in different plasma and neutral properties in the heliotail (Section 3.4), the overall size of the heliosphere for cases I, II, and III is remarkably similar. Only case IV, which has the biggest departure in source terms from the baseline, shows a significantly smaller heliotail.
To quantify the differences among the four cases on the front of the heliosphere, Table 1 shows the distance from the TS to the HP in both the Voyager 1 and Voyager 2 directions. Here we see that the models with suprathermal ions present in their velocity distributions (cases II, III, and IV) all have a slightly narrower heliosheath than the single Maxwellian case I. The differences are less than 10%, even for case IV, which almost doubles the rate at which energy is removed from the plasma. For comparison, Voyager 1 crossed the TS at ∼94 au, and the HP at ∼122 au, for a heliosheath thickness of ∼28 au. This observed heliosheath thickness cannot be explained by traditional steady-state MHDneutral models. Adding suprathermal protons can decrease the heliosheath thickness, but as we have shown, with the energydependent cross section included in the collision integrals, the effect of PUIs on the solution is moderate. Voyager 2 crossed the TS at ∼84 au, and crossed the HP at ∼119 au, for a heliosheath thickness of ∼35 au, which is much closer to our predictions. Overall, our steady-state results suggest that the increased removal of plasma pressure from the heliosheath due to charge-exchanging PUIs has a smaller affect on the thickness of the heliosheath than time-dependence due to the solar cycle and solar activity (e.g., Pogorelov et al. 2009bPogorelov et al. , 2013Pogorelov et al. , 2017Washimi et al. 2011Washimi et al. , 2017.

Discussion and Conclusions
We have performed simulations of the SW-LISM interaction in which we included the feedback of the presence of PUIs onto the plasma through their charge-exchange coupling to neutral hydrogen. We considered three models for PUIs in the subsonic SW, and compared these to a baseline case where the proton distribution is a single Maxwellian. We have previously demonstrated the importance of allowing neutral hydrogen to charge-exchange with suprathermal protons in order to obtain sufficient ENA flux to explain the observations by the IBEX spacecraft (e.g., Heerikhuisen et al. 2008). The calculations presented here represent an improvement over the multicomponent PUI models we employed previously (e.g., Zank et al. 2010;Zirnstein et al. 2014), because here we compute the charge-exchange of neutrals with the PUI components and apply the resulting energy and momentum source terms to the MHD equations.
Another important aspect to simulations with suprathermal ions comes from the energy-dependent charge-exchange cross section. Suprathermal protons are much more mobile than thermal protons, and as such would be expected to participate in more charge-exchange collisions (Heerikhuisen et al. 2015). However, above ∼10 keV the exponential factor in the chargeexchange cross section (Equation (4)) significantly decreases the collision rate, provided that the cross section is not assumed constant as it is in the analytic approximations generally used in global simulations of the heliosphere. The impact of including the energy-dependence of the cross section into the simulations is perhaps best seen by comparing the source terms obtained from a kappa-distribution to those obtained using a Maxwellian velocity distribution for protons. Figure 2 shows that the energy transfer rates from protons to hydrogen atoms are in fact lower in the kappa-distribution case, even though this case has a power-law tail to the velocity distribution. The model heliospheres that result from our simulations with the four different proton distributions turn out to be remarkably similar. Figures 4 and 5 show that the size of the heliosphere and the bulk properties of the plasma and neutral populations do not differ much for the models where we assume an equipartition between electron and proton temperatures (cases I, II, and III). Only the case where we allow the protons to carry more thermal energy than electrons (case IV) shows a noticeably smaller heliosphere with the heliotail cooling much more quickly and the plasma density about double that of the other cases in the distant heliotail.
The model with PUIs and unequal electron and proton temperatures represents the most realistic of the four cases we have considered, given that the electrons likely have a temperature similar to the core SW near the TS, with wave/ turbulent heating resulting in eventual equilibration farther down the heliotail. Such a model has not been implemented self-consistently before, and represents an interesting new approach to including the enhanced transport of energy from the ions to the neutrals in a simple way. We showed that allowing for a greater proton temperature as compared to the electron temperature roughly doubles the charge-exchange energy transfer rate (Figure 2). By extension, the flux of ENAs from the heliosheath would also increase by a comparable factor, perhaps going some way to make up for the lack of ENA flux from the heliotail in previous simulations. However, Zirnstein et al. (2017) implemented an approach where electrons were cooler than protons in their post-processing routine, and found that the ENA flux was still underestimated by a factor of two or three. We anticipate future simulations to develop this approach further.
The heliosheath thickness inferred from the Voyager 1 data (i.e., 28 au) is significantly smaller than any traditional steadystate MHD-neutral simulation of the heliosphere is able to produce. This has lead to a search for mechanisms that could remove additional pressure from the heliosheath. Other than the novel idea to look at a high thermal conduction limit (Izmodenov et al. 2014), most approaches involve the charge-exchange between neutral hydrogen and highly mobile suprathermal protons (e.g., Malama et al. 2006;Heerikhuisen et al. 2008;Pogorelov et al. 2016;Opher et al. 2018). Our results indicate, however, that the effect of enhanced chargeexchange with PUIs is greatly reduced once the energydependent cross section is taken into account, and the thickness of the heliosheath shrinks by at most 10% (Table 1), when all other factors remain unchanged. We conclude therefore, that time-dependence is likely the dominant factor in the apparent narrowness of the heliosheath in the V1 direction, as has also been demonstrated in recent time-dependent simulations (e.g., Pogorelov et al. 2017;Washimi et al. 2017). The fact that the heliosheath thickness in the V2 direction is at least 34 au (much closer to the model prediction of ∼40 au), lends additional support to this notion because the most likely interstellar magnetic field direction compresses the heliosphere more in the V2 than the V1 direction on average. Future work will focus on assigning a different breakdown of the PUI populations as a function of location on the TS. It is likely that the properties of the TS differ in response to changes in the radial distance, SW conditions, and the fraction of PUIs.
So in the future we will apply observational (e.g., ENA) and theoretical considerations to constrain the density and energy partitioning (Equations (1) and (2)) at the TS.
J.H. acknowledges support from a UAH IIDR grant. This work was supported by NASA grants NNX16AG83G, NNX14AP24G, NNX14AJ53G, 80NSSC18K1649NS, NNX14AF43G, DOE grant DE-SC0008334, and IBEX subaward SUB 0000167/NASA 80NSSC18K0237. G.P.Z. has been supported in part by the NSF EPSCoR RII-Track-1 Cooperative Agreement OIA-1655280. Computational resources were provided on Blue waters through NSF PRAC award OAC-1811176, on Comet and Stampede through XSEDE project MCA07S033, and on Pleiades through NASA project SMD-16-7570. This work was also partially supported by the IBEX mission as a part of NASA's Explorer program.