Slowdown and Heating of Interstellar Neutral Helium by Elastic Collisions beyond the Heliopause

Neutral components of the interstellar matter penetrate heliospheric boundaries and are used to diagnose the very local interstellar medium (VLISM) in front of the heliosphere. Comprehensive observations of interstellar neutral (ISN) helium with the IBEX-Lo sensor (Fuselier et al. 2009) on the Interstellar Boundary Explorer (IBEX; McComas et al. 2009) mission reveal the VLISM temperature and its relative velocity to the Sun (Möbius et al. 2009a, 2009b). Analyses of these observations have assumed that the main population of ISN helium is not modified outside the heliopause, and therefore the atom fluxes measured at 1 au can be directly related to their distribution function in the pristine VLISM, i.e., before it is affected by the heliosphere, using Liouville's theorem (Lee et al. 2012, 2015; Sokół et al. 2015). Additionally, atoms are ionized in the inner heliosphere due to charge exchange with solar wind ions, electron impact ionization, and photoionization (Sokół et al. 2019). These analyses showed that the ISN helium temperature is 7500 K and flows with a speed of 25.4 km s⁻¹ relative to the Sun toward direction (757, −51) in the ecliptic coordinates (Bzowski et al. 2015; McComas et al. 2015; Schwadron et al. 2015). These parameters are often interpreted as the pristine relative velocity and temperature of the VLISM because the primary ISN helium so far has not been expected to be significantly modified beyond the heliopause in the region perturbed by the heliosphere (the outer heliosheath).

The Warm Breeze discovery indicated that another ISN helium population also exists in the outer heliosheath (Kubiak et al. 2014, 2016) with a higher temperature and lower speed. Bzowski et al. (2017, 2019) showed that this population is created from He⁺ ions in the outer heliosheath that charge exchanged with the primary population of the ISN helium. Both populations are detected in the three lowest energy steps of IBEX-Lo (Swaczyna et al. 2018). Still, even with a comprehensive uncertainties system (Swaczyna et al. 2015), the observations show statistically significant differences from the predictions. This discrepancy may suggest that the model is too simple to fully represent the observations. Swaczyna et al. (2019a) suggested that the primary ISN helium is not fully equilibrated in the VLISM and, instead of a Maxwellian, follows a kappa distribution. Additionally, Wood et al. (2019) showed that a bi-Maxwellian distribution with different perpendicular and parallel temperatures (T_⊥/T_∥ = 0.62) relative to the inflow direction gives a significant improvement in the goodness of fit.

In the outer heliosheath, charge exchange collisions between helium atoms and He⁺ ions dominate all other ionization processes of helium atoms. The resonant character of this process is often used to claim that the momentum of the newly created atom is the same as the momentum of the parent ion. Swaczyna et al. (2019b) demonstrated that this simplification is not justified for typical collision speeds in the outer heliosheath, and angular scattering of colliding particles impacts the distribution function of the created neutral population. However, this angular scattering can be neglected in studies of the global heliosphere because the amount of momentum exchange is not significantly changed (Heerikhuisen et al. 2009; Izmodenov et al. 2000).

In this study, we focus on elastic collisions of helium atoms that are rarely considered in the heliospheric studies. While elastic collisions do not change charge state of particles, they lead to momentum exchange. The heating of interstellar atoms due to such collisions with the solar wind has been considered since the 1970s. Many of these studies (Wallis 1975; Holzer 1977; Fahr 1978; Kunc et al. 1983; Fahr et al. 1985) applied a continuous momentum transfer approach and concluded that this process might increase the temperature of ISN atoms at 1 au by a few hundred to a few thousand kelvins from the interstellar value. However, Gruntman (1986) argued that this approach should not be used as each atom undergoes a small number of binary collisions inside the heliosphere, and thus the continuous transfer approximation is not valid. Additionally, Chassefière et al. (1986) pointed out that Fahr et al. (1985) significantly overestimated the elastic cross section, consequently overemphasizing heating by elastic collisions with the solar wind. Recently, Gruntman (2013, 2018) used a "one or none" collision approximation for this process and showed that they increase ISN helium temperature by ∼200 K and slightly enhance the tail of the distribution due to rare but large scattering angles.

Chassefière & Bertaux (1987) considered elastic collisions outside the heliopause and concluded that about 10%–30% of ISN helium atoms are elastically scattered to form an additional slowed (by ∼3 km s⁻¹) and heated (by ∼4000 K) population. This population should not be confused with the secondary population created from the charge exchange collisions. They obtained this result using a momentum transfer cross section that significantly underestimates the number of collisions in the outer heliosheath but simultaneously increases momentum exchange in these collisions. In the present study, we revisit this problem and present how these collisions change the interpretation of IBEX observations.

Far from the heliopause, charged and neutral populations in the pristine VLISM are in thermal equilibrium and thus can be described by Maxwell distributions. However, in the outer heliosheath, the interstellar plasma diverges to flow around the heliopause, while ISN atoms can still penetrate the heliosphere (Parker 1961; Axford 1972; Zank 1999). Consequently, flows of ions and atoms are separated, so they are no longer in thermal equilibrium. While elastic collisions are not frequent enough to restore equilibrium, they affect the neutral components passing the outer heliosheath. In this study, we estimate this effect using Monte Carlo calculations.

2.1. Elastic Collision Scattering

Calculations of elastic collision scatterings require using relevant differential cross sections. The order of magnitude of integral cross sections for various elastic scatterings is similar. Hence, we consider the four most abundant species in the VLISM: hydrogen atoms, protons, helium atoms, and He⁺ ions. Collisions with less abundant species can be safely neglected due to orders of magnitude lower abundances in the VLISM (Frisch et al. 2011).

Separation of charge exchange and elastic scatterings for collisions of ions and atoms of the same species (e.g., He⁺–He⁰) is not possible from quantum mechanical considerations because the colliding nuclei are not distinguishable. However, differential cross sections calculated in the indistinguishable approach show two distinctive maxima at scattering angles between directions of the incoming and outgoing atom close to 0° and 180° (Krstić & Schultz 1998; Barata & Conde 2010; Schultz et al. 2016). The scattering amplitudes contributing to each maximum can be separated and used to calculate two separate elastic and charge exchange cross sections. This approximation does not account for interference of scattering amplitudes, but this effect can be neglected for collision speeds considered in this study. In this study, we follow the distinguishable particle approach (Appendix A). Moreover, we do not use the momentum transfer cross sections that can be used as an effective term for a continuous momentum transfer between populations because such cross section cannot be used to track individual collisions between particles.

Comprehensive sources of these cross sections are limited. Krstić & Schultz (1998) calculated differential cross sections for collisions of helium atoms with protons, which are available in the ALADDIN database ⁴ in a tabulated form for 24 center-of-mass energies in the range 0.1–100 eV. We do not find any other comprehensive set of the needed differential cross sections. Therefore, we calculate these cross sections using the Jeffreys–Wentzel–Kramers–Brillouin (JWKB) method as described by Nitz et al. (1987), utilizing interaction potentials for collisions of helium atoms with hydrogen atoms (Gao et al. 1989), protons (Helbig et al. 1970), helium atoms (Ceperley & Partridge 1986), and He⁺ ions (Marchi & Smith 1965; Barata & Conde 2010). The calculations are performed for a wide range of energies from 10⁻⁴ to 10⁴ eV, as described in Appendix A. We check that our calculations for collisions with protons are consistent with the differential cross sections obtained by Krstić & Schultz (1998). Figure 1 shows the obtained elastic collision cross sections of helium atoms for energies of 0.01, 0.1, 1, and 10 eV to span over the collision energy range expected in the outer heliosheath. Atom–ion differential cross sections show characteristic oscillations related to more partial wave contributions, as phase shifts decrease much more slowly with the order of the partial wave than in the case of atom–atom collisions. We check the obtained cross sections against previously obtained sources for the available energies and collision pairs (Gao et al. 1989; Krstić & Schultz 1998; Barata & Conde 2010). As the collision energy increases, the scattering angles in elastic collisions are smaller.

Zoom In Zoom Out Reset image size

We estimate the scattering effect on the relative motion of helium atoms colliding with these species, as presented in Figure 2. Panel (a) shows the integral cross section for these collisions as a function of relative speed. Collisions with protons and He⁺ ions show stronger dependence on the relative speed, but all are comparable for speeds relevant for this problem. Mean scattering angles, presented in panel (b), are significantly larger at lower speeds. Panel (c) displays the mean relative change of the speed in the helium atom initial motion direction. This change is between ∼0.05 and 1.5 km s⁻¹ per collision and depends on colliding species and relative speed. Such changes are not negligible compared to the bulk flow of the VLISM relative to the Sun. Moreover, all considered elastic cross sections are higher than the charge exchange cross section with He⁺ ions, thus we expect that these collisions affect a significant portion of the helium atom population. For this Letter, we only consider the consequences of more frequent elastic collisions. However, charge exchange collisions, mostly with He⁺ ions, result in creating the Warm Breeze and remove some atoms from the primary population (Bzowski et al. 2019), but this effect is not considered in this study. Consequences for the primary ISN helium should be much smaller than found for the primary ISN hydrogen by Pogorelov et al. (2008, 2009) due to lower ionization rate and narrower distribution of helium atoms in the velocity space.

Zoom In Zoom Out Reset image size

2.2. Outer Heliosheath Plasma and Neutrals

2.3. Transport of Neutral Helium Atoms

Following the method presented in Section 2.3, we investigate changes in the ISN helium distribution function in the outer heliosheath. The pristine VLISM populations follow the Maxwell distribution with the speed u_VLISM = 25.4 km s⁻¹ relative to the Sun, and the temperature T_VLISM = 7500 K. Figure 3 shows the resulting distribution function of the primary atoms at the heliopause. The probability distribution function is widened because of elastic collisions in the outer heliosheath but shows significant asymmetry in parallel and perpendicular directions to the inflow direction (z-axis). Moreover, the parallel component shows an extended tail in the antisunward direction. The average parallel component of the atoms' velocities is u_HP = 24.85 km s⁻¹, and their temperature is T_HP = 8850 K. The temperature describes the mean energy of the particles and is related to the mean square speed $\left\langle {v}^{2}\right\rangle$ as $T=\tfrac{m\left\langle {v}^{2}\right\rangle }{3{k}_{{\rm{B}}}}$, where m is the atomic mass, and k_B is the Boltzmann constant.

Zoom In Zoom Out Reset image size

The fitted Maxwell distribution has the relative speed u_M = 24.95 km s⁻¹, and the temperature is T_M = 8600 K. This function well describes the perpendicular component of the distribution, with only small discrepancies in the tails, but fails to describe the asymmetric parallel component. The peak of the Maxwell fit is shifted to lower relative speed than the peak, and the extended antisunward tail is not captured. However, the asymmetric kappa distribution introduced in Section 2.3 reproduces all these features. The peak position is at u_κ = 25.15 km s⁻¹. The perpendicular component has temperature T_κ⊥ = 8650 K, and kappa index κ_⊥ = 30. This finding shows that this component is very close to the Maxwell distribution, and only small differences are expected in the tails of the distribution (see Figure 3). The parallel part has different parameters for the sunward (T_κ∥S = 7950 K, κ_∥S = ∞) and antisunward (T_κ∥A = 10,600 K, κ_∥A = 6.5) parts of the distribution. Note that the sunward part is consistent with the Maxwell distribution, which is represented by κ_∥S = ∞ .

The evolution of the distribution parameters in the outer heliosheath is presented in Figure 4. The distribution of the primary ISN helium starts to change properties around 240 au from the Sun, i.e., where significant deviations of the plasma and neutrals flow from the pristine VLISM parameters begin in this model. The properties progressively change across the outer heliosheath, but the change is more rapid closer to the heliopause, where the discrepancy between parameters of the ISN helium and the plasma components is the highest.

Zoom In Zoom Out Reset image size

The mean numbers of collisions for one helium atom traveling through the outer heliosheath based on the Monte Carlo calculations and their standard deviations are 1.58 ± 1.30, 0.20 ± 0.45, 0.88 ± 1.05, and 1.70 ± 1.32 for collisions with protons, He⁺ ions, primary ISN hydrogen atoms, and secondary ISN hydrogen atoms, respectively. The average number of all elastic collisions is 4.35, and only 1.5% of primary ISN helium atoms are not elastically scattered in the outer heliosheath at all.

So far, analyses of IBEX observations assumed that the primary ISN helium follows the Maxwell distribution in the outer heliosheath (e.g., Bzowski et al. 2015; Schwadron et al. 2015; Swaczyna et al. 2018). These studies neglected elastic scatterings; therefore, the speed and temperature are different in the pristine VLISM compared to the parameters found previously. To estimate first-order corrections to the parameters of the VLISM derived from the IBEX observations, we assume that the absolute differences in the parameters are the same as in our simulations. Our analysis shows that the speed and temperature at the heliopause following the Maxwell fit change by Δu = u_M − u_VLISM = − 0.45 km s⁻¹ and ΔT = T_M − T_VLISM = 1100 K compared to the pristine VLISM. Assuming that previously derived values correspond to conditions at the heliopause with the Maxwell distribution function u_HP,IBEX = 25.4 km s⁻¹ and T_HP,IBEX = 7500 K, the pristine VLISM parameters are u_VLISM,IBEX = u_HP,IBEX − Δu = 25.85 km s⁻¹ and T_VLISM,IBEX = T_HP,IBEX − ΔT = 6400 K. Furthermore, we expect the ISN helium distribution at the heliopause is better represented by the asymmetric kappa distribution with parameters u_HP,κ = 25.6 km s⁻¹, T_HP,κ⊥ = 7550 K, κ_HP,⊥ = 30, T_HP,κ∥S = 6800 K, κ_HP,∥S = ∞ , T_HP,κ∥A = 9500 K, and κ_HP,∥A = 6.5. To verify if these changes to initial conditions can change our findings, we repeat our calculations with an ad hoc model in which flow speeds are proportionally scaled to match u_VLISM,IBEX at the outer boundary, and all temperatures beyond the heliopause are decreased by ΔT. The obtained relative changes Δu and ΔT agree with these found above.

The changes in the outer heliosheath significantly exceed their uncertainties from previous IBEX studies (Bzowski et al. 2015; Swaczyna et al. 2018). To illustrate how they affect the IBEX-Lo observations, we calculate the expected count rates for three orbits observed in 2009: 14, 16, and 18 for the three distributions characterized above. These rates are calculated using the analytic full integration model (Schwadron et al. 2015). Figure 5 shows results compared with the observed rates on these orbits. The Maxwell distribution with the pristine VLISM parameters shows a significantly narrower signal in spin angle, representing a lower ratio of thermal speed to bulk speed in this case. The asymmetric kappa distribution and the Maxwell fit at the heliopause give similar rates. However, small changes are present in the distribution wings.

Zoom In Zoom Out Reset image size

The vast majority of primary ISN helium atoms from the VLISM are elastically scattered at least once in the outer heliosheath before entering the heliosphere. Consequently, ISN helium atoms are significantly heated (∼1100 K) and slowed (∼0.45 km s⁻¹), thus their distribution at the heliopause does not represent the pristine VLISM parameters, as is commonly assumed (Möbius et al. 2004). This result shows that the pristine VLISM parameters obtained from IBEX observations need to be adjusted for this effect. Using the current "working values" provided by McComas et al. (2015), we estimate the pristine VLISM relative speed is u_VLISM ≈ 25.4 + 0.45 = 25.85 km s⁻¹, and the temperature is T_VLISM ≈ 7500–1100 = 6400 K. Note that this change is across the IBEX correlation tube of these two parameters (McComas et al. 2012, 2015). Additionally, the estimated heating is about ∼5 times stronger than heating due to collisions with the solar wind protons and α-articles, as estimated by Gruntman (2018).

Our results differ significantly from the findings by Chassefière & Bertaux (1987), who estimated that only 10%–30% of atoms are elastically scattered and produce an additional component of ISN helium slower by ∼3 km s⁻¹ and hotter by ∼4000 K compared to the primary population. However, the average changes of both populations in their study would yield a speed 20% × 3 = 0.6 km s⁻¹ and temperature 20% × 4000 = 800 K, which are close to our results. We note that those authors used momentum transfer cross sections that should give correct average changes in momentum but underestimate the number of elastic collisions (Schultz et al. 2016).

The elastic collisions cause a significant asymmetry in the ISN helium distribution at the heliopause naturally aligned with the interstellar bulk flow (see Figures 3 and 4). This alignment is a result of collisions occurring at different positions in the outer heliosheath with slower and heated populations. The magnitude of the temperature change is smaller than the asymmetry found by Wood et al. (2019), who used a bi-Maxwell distribution to model the ISN helium distribution. However, this difference may be caused by two vastly different distribution functions used in the respective studies.

This investigation gives an estimation of the elastic collision effects on the distribution of the ISN helium atoms. Nevertheless, a full quantitative analysis using IBEX observations to determine the pristine VLISM parameters is necessary and requires tracking the elastic scattering across the outer heliosheath using a 3D global model of the heliosphere. Such a study would allow for the determination of possible deflection of the inflow speed in the outer heliosheath that cannot be estimated using the method presented in this Letter. Finally, this Letter also serves as a source of the cross sections for helium atom elastic collisions in the outer heliosheath.

This material is based upon work supported by the National Aeronautics and Space Administration under grant No. 80NSSC20K0781 issued through the Outer Heliosphere Guest Investigators Program. This work is also partially funded by the IBEX mission as part of the NASA Explorer Program (80NSSC20K0719).

Differential cross sections for elastic collisions are calculated using the semiclassical JWKB approximation, using known interaction potentials for the considered pairs of particles. Phase shifts (δ_l ) are computed with the integral form given in Nitz et al. (1987, Equations (A1) and (A4)). We calculate phase shifts for each energy until their magnitudes for ten consecutive orbital angular momenta l are smaller than 10⁻⁵ radians. Scattering amplitudes are calculated for 768 values of center-of-mass scattering angles ${\theta }_{m}=\tfrac{\pi }{2}(1+{x}_{m})$, where x_m is the mth root of the Legendre polynomial of degree 768:

$$\begin{eqnarray}&&f\left(\theta \right)=\displaystyle \sum _{l}\displaystyle \frac{\left(2l+1\right)}{2{ik}}\left({e}^{2i{\delta }_{l}}-1\right){P}_{l}\left(\cos \theta \right),\end{eqnarray} \tag{ 1 }$$

where k is the wavenumber and P_l are the Legendre polynomials. This selection of angles constitutes abscissas for the Gauss–Legendre quadrature and is the same as used in previous comprehensive studies of differential cross sections (e.g., Krstić & Schultz 1998; Schultz et al. 2016). Finally, differential cross sections are given as $d\sigma (\theta )/d{\rm{\Omega }}={\left|f\left(\theta \right)\right|}^{2}$.

In the case of a collision between He⁺ ions and He atoms, we calculate separate scattering amplitudes f^g(θ) and f^u(θ) for the gerade and ungerade ground states, respectively. Moreover, the nuclei are not distinguishable and thus the charge transfer amplitude in angle π − θ and elastic scattering amplitude in angle θ can interfere. Consequently, the total differential cross section is

$$\begin{eqnarray}&&\displaystyle \frac{d\sigma }{d{\rm{\Omega }}}=\displaystyle \frac{{\left|{f}^{{\rm{g}}}\left(\theta \right)+{f}^{{\rm{u}}}\left(\theta \right)+{f}^{{\rm{g}}}\left(\pi -\theta \right)-{f}^{{\rm{u}}}\left(\pi -\theta \right)\right|}^{2}}{4},\end{eqnarray} \tag{ 2 }$$

but elastic collisions and charge exchange (transfer) cross section can be approximately separated and expressed as

$$\begin{eqnarray}&&\displaystyle \frac{d{\sigma }_{\mathrm{el}}}{d{\rm{\Omega }}}=\displaystyle \frac{{\left|{f}^{{\rm{g}}}\left(\theta \right)+{f}^{{\rm{u}}}\left(\theta \right)\right|}^{2}}{4},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d{\sigma }_{\mathrm{ct}}}{d{\rm{\Omega }}}=\displaystyle \frac{{\left|{f}^{{\rm{g}}}\left(\pi -\theta \right)-{f}^{{\rm{u}}}\left(\pi -\theta \right)\right|}^{2}}{4},\end{eqnarray} \tag{ 4 }$$

respectively. Since the elastic and charge exchange cross sections show significant peaks for angles close to 0° and 180°, respectively, the interference term accounted for in Equation (2) but neglected with the approximation in Equations (3)–(4) is negligible.

All considered cross sections are calculated for energies from 10⁻⁴ to 10⁴ eV. As the cross section shows significant changes for low collision energy, we calculated these cross sections with finer energy resolutions at lower energies. For ranges 10⁻⁴–10⁻², 10⁻²–10⁰, 10⁰–10², and 10²–10⁴ eV, we use a logarithmically spaced grid with Δ(log₁₀ E) = 0.005, 0.01, 0.05, and 0.1, respectively. All differential cross sections calculated in this work have been posted on Zenodo under a Creative Commons Attribution license: doi:10.5281/zenodo.4555716.

Monte Carlo calculations of elastic scatterings and their consequences for the primary ISN helium population start with a test particle sample. The sample of n = 10⁶ helium atom velocities ${\left\{{{\boldsymbol{v}}}_{i,0}\right\}}_{i=1,...,n}$ is drawn from the Maxwell distribution with the temperature and relative speed assumed for the pristine VLISM ISN population in the global heliosphere model. Their initial position is ${\left\{{{\boldsymbol{x}}}_{i,0}\right\}}_{i=1,..,n}=\left(0,0,{d}_{0}+{\delta }_{i}\right)$, where d₀ = 400 au is set slightly farther away than the point at which the flows start to diverge from the pristine values (i.e., near the bow wave), and δ_i is randomly chosen from the range [0, Δd], where Δd is the calculation length step. Collisions of each particle from the sample are tracked until the particle position along the z-axis is closer than the heliopause at 110 au.

The probabilities of collisions are calculated with the following formula:

$$\begin{eqnarray}&&{p}_{k}=\displaystyle \frac{{\rm{\Delta }}{{dn}}_{k}{\sigma }_{k}\left({v}_{\mathrm{rel},k}\right){v}_{\mathrm{rel},k}}{{v}_{\mathrm{atom}}},\end{eqnarray} \tag{ 5 }$$

where n_k is the density of population k, v_rel,k is the mean relative speed between the tracked atom and this population, σ_k is the integral cross section for the elastic collision, and v_atom is the tracked atom speed in the Sun frame. The mean relative speed is calculated assuming that the collision partners are distributed following the Maxwell distribution (e.g., Ripken & Fahr 1983):

$$\begin{eqnarray}&&{v}_{\mathrm{rel},k}={u}_{{\rm{T}},k}\left[\left(\omega +\displaystyle \frac{1}{2\omega }\right)\mathrm{erf}\omega +\displaystyle \frac{\exp \left(-{\omega }^{2}\right)}{\sqrt{\pi }}\right]\end{eqnarray} \tag{ 6 }$$

with ${u}_{{\rm{T}},k}=\sqrt{2{k}_{{\rm{B}}}{T}_{k}/{m}_{k}}$, the thermal speed, and $\omega =\left|{{\boldsymbol{v}}}_{\mathrm{atom}}-{{\boldsymbol{u}}}_{k}\right|/{u}_{{\rm{T}},k}$, the ratio of the relative speed of the atom v _atom and the bulk speed of the collision partner population u _k to the thermal speed.

The formula given in Equation (6) is precise if the integral cross section changes linearly with the relative speed. This approximation is not good for helium atom collisions with protons and He⁺ ions, which show significant variations for low collision speeds. However, we verified that this formula could be used with a parabolic fit to the cross section in log–log scale for energies below 100 eV, which is shown with dashed lines in Figure 1(a). We verified that the relative error of this approach is less than 4% for the temperature range considered in this study.

The calculation length step is set as Δd = 3 au, which assures that the probability of elastic collision with each species is much smaller than 1, so we can effectively neglect the likelihood that more than one elastic collision occurs in one length step. Additionally, this step is small enough to ignore changes in the population bulk speed and temperature. We randomly choose whether a collision with any species occurs with the probability calculated with Equation (5). The probability that no collision occurs is given as 1 − ∑_k p_k . In this situation, the atom velocity and position are propagated:

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{i,j+1}={{\boldsymbol{v}}}_{i,j},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{i,j+1}={{\boldsymbol{x}}}_{i,j}+{\rm{\Delta }}d\displaystyle \frac{{{\boldsymbol{v}}}_{i,j}}{\left|{{\boldsymbol{v}}}_{i,j}\right|},\end{eqnarray} \tag{ 8 }$$

where j enumerates positions along the atom trajectory.

If a collision with any of the considered species occurs, we draw 10 possible velocities of the collision partner from the Maxwell distribution with the speed and temperature given by the numerical model in position x _i,j . We verified with a smaller number of modeled atom trajectories that this sample of velocities is enough to represent the difference in the collision probability with different parts of the distribution. Subsequently, we select one of these speeds using probability weights given by the integral cross section for the relative speed of the tracked atom and the collision partner. This two-step random selection provides that we account for differences in the integral cross sections with different parts of the distribution function. Further, we select the scattering angle (θ) using the weighting proportional to $\sin \left(\theta \right)d\sigma /d{\rm{\Omega }}$, where dσ/dΩ is the differential cross section for one of the two closest energy bins for which the cross section is calculated, selected with a weight inversely proportionally to the energy difference between the collision energy and the bin energy. The azimuthal angle (ϕ) is randomly chosen from the range [0, 2π] with uniform probability. Finally, the fractional position of the collision ξ is selected from the range [0, 1]. The post-collision velocity and position given in this case are

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{i,j+1}=\displaystyle \frac{{m}_{\mathrm{atom}}}{{m}_{\mathrm{atom}}+{m}_{k}}{\boldsymbol{r}}\left({{\boldsymbol{v}}}_{i,j}-{{\boldsymbol{v}}}_{\mathrm{CM}},\theta ,\phi \right)+{{\boldsymbol{v}}}_{\mathrm{CM}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{i,j+1}={{\boldsymbol{x}}}_{i}+\xi {\rm{\Delta }}d\displaystyle \frac{{{\boldsymbol{v}}}_{i,j}}{\left|{{\boldsymbol{v}}}_{i,j}\right|}+\left(1-\xi \right){\rm{\Delta }}d\displaystyle \frac{{{\boldsymbol{v}}}_{i,j+1}}{\left|{{\boldsymbol{v}}}_{i,j+1}\right|}.\end{eqnarray} \tag{ 10 }$$

The rotation function ${\boldsymbol{r}}\left({\boldsymbol{v}},\theta ,\phi \right)$ returns scattered velocity in the center-of-mass frame (moving with velocity v _CM) and is defined in Appendix A in Swaczyna et al. (2019b). Each atom is tracked until its position along the z-axis is less than 110 au (i.e., the heliopause distance in this direction). Properties of the atom distribution are tracked in 10 au wide bins along the z-axis, and from each trajectory, we select one speed for each bin.

The distribution function of the primary ISN helium atoms in the outer heliosheath is described using an asymmetric kappa distribution, in which the perpendicular components of velocities are described using a two-dimensional kappa distribution, and the parallel part is composed from two one-dimensional kappa distributions separate for the atoms with speeds below and above the distribution peak. The velocity components perpendicular to the inflow direction are denoted as v_x and v_y , and the parallel component is along the z-axis, v_z . The perpendicular distribution function is (Livadiotis & McComas 2013)

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\perp }\left({v}_{{\rm{x}}},{v}_{{\rm{y}}};{\theta }_{\perp },{\kappa }_{\perp }\right)\\ =\displaystyle \frac{1}{\pi {\theta }_{\perp }^{2}}\displaystyle \frac{{\kappa }_{\perp }\ -\tfrac{1}{2}}{{\kappa }_{\perp }-\tfrac{3}{2}}{\left(1+\displaystyle \frac{1}{{\kappa }_{\perp }-\tfrac{3}{2}}\displaystyle \frac{{v}_{{\rm{x}}}^{2}+{v}_{{\rm{y}}}^{2}}{{\theta }_{\perp }^{2}}\right)}^{-{\kappa }_{\perp }-\displaystyle \frac{1}{2}},\end{array}\end{eqnarray} \tag{ 11 }$$

where θ_⊥ is the speed scale parameter and κ_⊥ is the kappa index in the perpendicular directions. Temperatures are related to the speed scale parameters as $\theta =\sqrt{\tfrac{2{k}_{{\rm{B}}}T}{{m}_{\mathrm{He}0}}}$.

The parallel distribution function is given in the form

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\parallel }\left({v}_{{\rm{z}}};u,{\theta }_{\parallel {\rm{S}}},{\theta }_{\parallel {\rm{A}}},{\kappa }_{\parallel {\rm{S}}},{\kappa }_{\parallel {\rm{A}}}\right)\\ =\left(\displaystyle \frac{{\theta }_{\parallel {\rm{S}}}}{2}\sqrt{\pi \left({\kappa }_{\parallel {\rm{S}}}-\displaystyle \frac{3}{2}\right)}\displaystyle \frac{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{S}}}-\tfrac{1}{2}\right)}{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{S}}}\right)}\right.\\ +{\left.\displaystyle \frac{{\theta }_{\parallel {\rm{A}}}}{2}\sqrt{\pi \left({\kappa }_{\parallel {\rm{A}}}-\displaystyle \frac{3}{2}\right)}\displaystyle \frac{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{A}}}-\tfrac{1}{2}\right)}{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{A}}}\right)}\right)}^{-1}\\ \times \left\{\begin{array}{cc}{\left(1+\displaystyle \frac{1}{{\kappa }_{\parallel {\rm{A}}}-\tfrac{3}{2}}\displaystyle \frac{{\left({v}_{{\rm{z}}}-u\right)}^{2}}{{\theta }_{\parallel {\rm{A}}}^{2}}\right)}^{-{\kappa }_{\parallel {\rm{A}}}} & {v}_{{\rm{z}}}\lt u\\ {\left(1+\displaystyle \frac{1}{{\kappa }_{\parallel {\rm{S}}}-\tfrac{3}{2}}\displaystyle \frac{{\left({v}_{{\rm{z}}}-u\right)}^{2}}{{\theta }_{\parallel {\rm{S}}}^{2}}\right)}^{-{\kappa }_{\parallel {\rm{S}}}} & {v}_{{\rm{z}}}\gt u\end{array}\right.,\end{array}\end{eqnarray} \tag{ 12 }$$

where u is the peak position in the parallel direction, θ_∥A and θ_∥S are the speed scale parameters in the sunward and antisunward directions, and κ_∥A and κ_∥S are the kappa indices in these two directions. The combined distribution is given as a product of these two distribution functions:

$$\begin{eqnarray}&&\begin{array}{l}f\left({v}_{{\rm{x}}},{v}_{{\rm{y}}},{v}_{{\rm{z}}};u,{\theta }_{\perp },{\theta }_{\parallel {\rm{S}}},{\theta }_{\parallel {\rm{A}}},{\kappa }_{\perp },{\kappa }_{\parallel {\rm{S}}},{\kappa }_{\parallel {\rm{A}}}\right)\\ \ \ ={f}_{\perp }\left({v}_{{\rm{x}}},{v}_{{\rm{y}}};{\theta }_{\perp },{\kappa }_{\perp }\right)\times {f}_{\parallel }\left({v}_{{\rm{z}}};u,{\theta }_{\parallel {\rm{S}}},{\theta }_{\parallel {\rm{A}}},{\kappa }_{\parallel {\rm{S}}},{\kappa }_{\parallel {\rm{A}}}\right).\end{array}\end{eqnarray} \tag{ 13 }$$

This distribution function has seven scalar parameters. We find that the sunward part of the parallel component is consistent with the Maxwell distribution, or in other words, that κ_∥S is very large. Therefore, we use the limit κ_∥S → ∞ for the parallel component:

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\parallel }\left({v}_{{\rm{z}}};u,{\theta }_{\parallel {\rm{S}}},{\theta }_{\parallel {\rm{A}}},{\kappa }_{\parallel {\rm{S}}}\to \infty ,{\kappa }_{\parallel {\rm{A}}}\right)\\ \ \ ={\left(\displaystyle \frac{{\theta }_{\parallel {\rm{S}}}}{2}\sqrt{\pi }+\displaystyle \frac{{\theta }_{\parallel {\rm{A}}}}{2}\sqrt{\pi \left({\kappa }_{\parallel {\rm{A}}}-\displaystyle \frac{3}{2}\right)}\displaystyle \frac{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{A}}}-\tfrac{1}{2}\right)}{{\rm{\Gamma }}\left({\kappa }_{\parallel {\rm{A}}}\right)}\right)}^{-1}\\ \quad \times \left\{\begin{array}{cc}{\left(1+\displaystyle \frac{1}{{\kappa }_{\parallel {\rm{A}}}-\tfrac{3}{2}}\displaystyle \frac{{\left({v}_{{\rm{z}}}-u\right)}^{2}}{{\theta }_{\parallel {\rm{A}}}^{2}}\right)}^{-{\kappa }_{\parallel {\rm{A}}}} & {v}_{{\rm{z}}}\lt u\\ \exp \left(-\displaystyle \frac{{\left({v}_{{\rm{z}}}-u\right)}^{2}}{{\theta }_{\parallel {\rm{S}}}^{2}}\right) & {v}_{{\rm{z}}}\gt u\end{array}\right..\end{array}\end{eqnarray} \tag{ 14 }$$

We find the distribution function presented above to effectively match the numerically obtained distribution with the Monte Carlo calculations.