Transport of Interstellar Neutral Helium throughout the Heliosphere

A number of physical processes accompanying the solar wind interaction with the local interstellar medium (LISM) are governed by charge exchange between ions and neutral atoms of interstellar origin. A new, 3D, MHD-plasma/kinetic-neutral model is developed that self-consistently includes both neutral hydrogen and helium atoms, and their feedback on the plasma, through charge exchange and photoionization. Focusing on the transport of interstellar neutral helium, quantitative estimates are provided for bulk properties, deflection angles, and velocity distribution functions (VDFs) along the upwind direction. It is shown that the average deflection of secondary He atoms born in the outer heliosheath (OHS) from their original direction in the LISM is ∼12° in front of the heliopause, and occurs in the directions parallel to the plane formed by the velocity and magnetic field vectors in the unperturbed LISM. While these properties are consistent with Interstellar Boundary Explorer observations of the “warm breeze,” we show that charge exchange in the OHS leads to remarkable deviations of their VDF from the Maxwellian distribution. He atom filtration in the OHS results in a significant temperature anisotropy and VDF asymmetries, even for the primary helium atoms that experience no charge exchange at all. This is an entirely kinetic phenomenon that shows that primary He atoms observed at 1 au have distributions substantially different from those in the LISM.


Introduction
Charge exchange between ions and neutral atoms has long been acknowledged to be of fundamental importance in the solar wind (SW) interaction with the partially ionized local interstellar medium (LISM) and the resulting formation of the heliosphere (Blum & Fahr 1969;Wallis 1971Wallis , 1975. Interstellar neutral (ISN) helium is of particular interest, since its large mean free path for charge-exchange collisions with other elements, relatively high abundance at 1 au, and ionization properties make it ideal for inferring LISM properties (Möbius et al. 2004(Möbius et al. , 2009Witte 2004;McComas et al. 2015). It was recently shown that due to temperature anisotropy ISN He cannot be described by a single Maxwellian distribution function at 150 au (Wood et al. 2019).
Interstellar Boundary Explorer (IBEX) observations resulted in the identification of a second population of He atoms ). This population, dubbed the "Warm Breeze" (WB), is estimated to be 30% hotter and 50% slower, at 150 au, than the primary ISN He flow, and deflected by ∼8°from the original LISM flow direction. Although four possible explanations were proposed for the WB origin by Kubiak et al. (2014), recent studies (e.g., Swaczyna et al. 2015;Kubiak et al. 2016Kubiak et al. , 2019Bzowski et al. 2017) provided strong evidence that the WB originates from charge exchange between the pristine ISN He and He + ions occurring in the outer heliosheath (OHS), the region of interstellar matter between the heliopause (HP) and the bow shock or, more likely, the bow wave (BW) (e.g., Pogorelov et al. 2008;Zank et al. 2012). The WB's deflection occurs predominantly in the plane passing through the Sun, formed by the interstellar magnetic field (ISMF) and velocity vectors (the B-V plane) (Kubiak et al. 2016). It was suggested by Lallement et al. (2005Lallement et al. ( , 2010, and confirmed in Pogorelov et al. (2008), that this is the plane of dominant deflection of H atoms. One approach to model the He transport is based on Müller & Cohen (2012), Sokół et al. (2015), Bzowski et al. (2017Bzowski et al. ( , 2019, and Kubiak et al. (2019). It maps the local phase-space distribution to a Maxwellian distribution at infinity and includes He atom production and losses and gravity. The properties of He + ions were estimated using the proton properties from different global heliosphere models, but helium populations were not self-consistent components of the model.
In this Letter, we present the first results from our 3D MHDplasma/kinetic-neutral model, which includes the effect of ISN He on the plasma. Our simulations are performed by adding He ion and atom transport to the Multi-Scale Fluid-Kinetic Simulation Suite (e.g., Pogorelov et al. 2014Pogorelov et al. , 2021. We investigate the bulk properties of the pristine and secondary populations of ISN He, including their deflection and velocity distribution functions (VDFs) along the upwind direction. Our simulations demonstrate the deviation of the primary, and especially of the secondary, He atom distributions from the isotropic Maxwellian distribution. We show that the filtration process results in the temperature anisotropy and asymmetries in the VDFs.
Our results demonstrate the importance of secondary He atoms for the SW-LISM interaction, and help improve the accuracy of derivation of LISM properties from IBEX data. We also show that selective exclusion of He atoms from their subset that experiences no charge exchange affects the resulting distribution function at IBEX. This is an entirely kinetic phenomenon that was omitted in the recent identification of the effect of elastic collisions of He atoms with protons in the OHS ).

Model
We model the SW-LISM interaction in the presence of H and He neutral atoms using a self-consistent, fully 3D, MHD-plasma/ kinetic-neutral model. The plasma description is based on the ideal MHD treatment for the mixture of charged particles (Pogorelov et al. 2004(Pogorelov et al. , 2017. Our approach is to solve the conservation laws for the mixture, now involving H + , He + , and electrons. Pickup ions exist, but are not treated separately at this point. In addition, we solve everywhere a set of two equations for the He + ions number density, + n He , and pressure, He He where γ = 5/3, u is the plasma velocity, and + S n He and + S p He are the source terms due to the He + + He charge exchange and photoionization. These sources also enter the mixture equation, in addition to the source terms describing the H + + H charge exchange. We assume that He + , protons, and electrons are comoving. Here Q C describes the thermal equilibration process due to Coulomb collisions between protons and He + ions:  (Malama 1991;Baranov & Malama 1993;Lipatov 2002;Izmodenov et al. 2003;Heerikhuisen et al. 2006Heerikhuisen et al. , 2008. Typical He + He + collisions occur in the OHS in the energy range 0.1-10 eV. According to Scherer et al. (2014), the process He + He + → He + + He at such low energies is the dominant for He in the OHS. Therefore, we consider only these collisions. Barnett et al. (1990) provide us with the charge-exchange cross sections, σ cx (He + He + ), in the energy range from 0.16 eV (0.04 eV amu −1 ) to 2 MeV. 4 The fit to Barnett's data yields to where E is the collision energy in keV. We fit the data in the range of 0.16 eV-2 keV. The maximum relative error of our fit is within 1% for E < 10 eV, and within 1.5% for E < 3 keV.
Since the cross sections for He + H + collisions are two orders of magnitude smaller, we neglect such collisions. We use the equation of state for any ion species, p s = n s k B T s , to obtain the temperature of protons and helium ions. The source terms for density, momentum, and energy (S e , S m , and S ρ , respectively) enter the mixture equations, while the He + pressure source term in Equation (2) is computed as in Pogorelov et al. (2016) (see also DeStefano & Heerikhuisen 2020): He He He For the moment, newborn neutral atoms are derived using the Maxwell-Boltzmann distributions for both H + and He + ions. Helium atoms are affected by the Sun's gravity, whereas an exact balance of gravity and radiation pressure is assumed for the H atoms. Bzowski et al. (2013) showed that photoionization is the dominant loss mechanism for He atoms in the SW, the electronimpact ionization rate being nearly one order of magnitude smaller. The photoionization rate is The ISN H and He densities are chosen to be n ∞,H = 0.17 cm −3 and n ∞,He = 1.53 × 10 −2 cm −3 , respectively. Notice that the relatively large hydrogen density is in agreement with the recent estimate by Swaczyna et al. (2020). The temperatures of both plasma and neutral gas are 7500 K in the unperturbed LISM, and we assume the Maxwell-Boltzmann distribution of the neutrals there.
We restrict ourselves to steady-state solutions corresponding to a spherically symmetric SW. Consistent with Zirnstein et al. (2016), we choose the following quantities at 1 au: proton density = + n 5.74 H cm −3 , speed V = 450 km s −1 , temperature T = 80,000 K, and magnetic field B = 37.5 μG. The He + density at the inner boundary is negligible as compared to protons. We do not include alpha particles in the model; therefore, the production of neutral He in the SW is essentially absent. The MHD equations are solved on a Cartesian grid with the resolution of 1.25 au in each direction near the heliospheric nose and 0.625 au within 20 au from the Sun The computational domain is a cube, extending 1000 au into the tail direction and 680 au in the upwind direction.

Simulation Results
To present our results in the B-V plane, we choose a reference frame with the x-axis in the direction antiparallel to V ∞ , the y-axis parallel to V ∞ × B ∞ , i.e., perpendicular to the B-V plane. The z-axis is therefore in the B-V plane. Figure 1 shows 2D distributions of number density (top panels), the in-plane deflection (middle panels), and bulk velocity magnitude (bottom panels) of the neutral He. The left, middle, and right panels correspond to the total (He), primary (He 0 ), and secondary (He 1 ) populations, respectively. It should be noticed that a bow shock is absent in this simulation. Instead, a bow wave exists in the form of weak discontinuity, i.e., at the BW the gradient of all plasma quantities increases abruptly. This condition disappears for larger interstellar magnetic fields (e.g., Pogorelov et al. 2017). A neutral atom is tagged as "secondary" if it originates in the LISM where > ¥ + T T 1.027 He . The outer boundary of the He 1 production region lies then in the BW precursor. The heliocentric distance of the HP is equal to 110 ± 2 au in the x-direction and 120 ± 2 au along the Voyager 1 2 ) are shown. From the left to the right, we show the results for the total, pristine, and secondary helium. The orientation of the LISM velocity and ISMF vectors is shown in panel (a). The adopted definition of positive deflection for a generic particle with speed v He is indicated in panel (d). The distributions are computed using 256 2 bins in the B-V plane, each of them with the thickness of 10 au in the y-direction and Δx = Δz = 6 au. Note. These statistics are computed using truncated conical volumes (20°aperture, Δx = 50 au) centered over the x-direction. Extrapolated quantities at 1 au are obtained from a linear fit to the statistics of Figure 2, in the range 8-20 au.
trajectory. The particle deflection in the B-V plane, θ P , is determined with respect to the −x-direction (Figure 1(d)). The direction toward zero average deflection of He 1 atoms at the HP is ∼26°.5 off the x-axis, in the direction (238°, 25°.4) (Figure 1(f)). This does not match the direction toward the minimum bulk speed of He atoms, which makes an angle of about 13°.5 (panel (i)) with the x-axis. This direction is also close to that of the He + pressure maximum in the OHS (245°.5, 25°.6). Table 1 summarizes the average helium properties in four regions between 10 and 450 au in the upwind direction, and extrapolated statistics at 1 au. Figure 2 shows the linear distributions in the direction of the LISM flow. The maximum density in the helium wall is observed at ∼25 au from the HP, where n He ≈ 1.065 n He,∞ and the density of the secondary He is » ¥ n n 0.15 He He, 1 (Figure 2(a)). The density of the primary He in the OHS decreases linearly from the value at infinity to » ¥ n n 0.92

He
He, 0 at the HP. In the SW, the filtration process is such that the densities of both He 0 and He 1 atoms decrease with decreasing heliocentric distance. Both the Sun's gravity and selective filtration are responsible for the absolute increase of the velocity component parallel to V ∞ , v P ≡ v x (Figure 2(b)). The average sunward velocity component of secondary He reaches the minimum of ∼12 km s −1 in the He wall. Panel (c) shows the parallel temperature, computed as where the angular brackets mean averaging over an ensemble of particles inside the volume, and V x is the bulk speed. Panel (d) shows the perpendicular temperature, In the SW, our simulations show that both the primary and secondary ISN He experience parallel cooling and perpendicular heating with decreasing heliocentric distance. The secondary He shows full temperature anisotropy (see panels (c) and (d)). In the upwind direction, the maximum parallel temperature in the OHS is  T ,He 1 ≈ 45,000 K. The maxima of perpendicular temperatures are T z,He 1 ≈ 15,200 K and T y,He 1 ≈ 13,000 K, achieved at ∼25 au and ∼40 au from the HP. The perpendicular anisotropy grows in the OHS with decreasing heliocentric distance, so that up to HP, and decreases to ∼0.8 at 10 au. The decrease in the parallel temperature is related to the narrowing and distortion of the v x distribution function (Figure 4(c)). This cooling effect is due to both the selective filtration in the OHS and gravitational focusing.
It is worth noting that our model allows for different temperatures for H + and He + ions. In the absence of thermalization, the temperature ratio at the HP would be » + +  (4)) and ∼10 −13 cm −3 , respectively, and the mean free path is in the range of 0.8-10 au. Consequently, charge-exchange collisions do not introduce temperature differences. Figure 3 shows the 2D probability densities of ISN He deflections in the SW and OHS regions. Our calculation confirms that the average deflection of the ISN He gas occurs in the B-V plane. The average out-of-plane deflection is negligible for both populations at all radial distances (see Table 1). Moreover, the filtration process of the primary He does not produce any noticeable average deflection. On average, the secondary helium instead is deflected by  q »  12 ,He 1 in the OHS at the heliospheric nose. As a result of the filtration process, in the SW the average deflection is smaller, becoming ∼6°. The distribution functions of secondary He show a remarkable asymmetry in the B-V plane, being skewed toward positive deflections. Particles with deflections up to ±180°in Figure 3 belong to the indirect beams (Fahr 1978) of secondary particles with sufficiently small v x and/or small impact parameter to be deflected by the Sun back into the upwind direction. The distributions in Figure 3 are in qualitative This is clearly seen in the VDFs in Figure 4. Panels (a) and (b) and panels (e) and (f) show the joint (v x , v z ) distributions for the total and the secondary He. These distributions are normalized so that ( ) x z x z He 2D , and fits to the Maxwellian distribution are shown with gray curves. Panels (c) and (g) show the v z -integrated distributions, while panels (d) and (h) present v x -integrated distributions. The pristine He is well described by a Maxwellian distribution in the narrow range of speeds ±15 km s −1 from the peak. In the SW, the filtration process results in a depletion of the VDF tails (panels (c) and (d)). Indirect beams of particles generate a tail in the antisunward side of VDF[v x ] and in the perpendicular VDFs. Although their density is negligible at R > 10 au, they may become important at smaller distances. We disregarded these particles in the temperature computation. The perpendicular temperature anisotropy of He 1 atoms can be deduced from . Such large indices make these distributions close to Maxwellian, but only near the peak.

Discussion and Conclusions
We have extended our MHD-plasma/kinetic-neutral model to reproduce features of the He atom transport, such as the He wall and focusing cone in the downwind direction. As shown by Zirnstein et al. (2021) the inclusion of He atoms may substantially affect the choice of the densities of protons and H atoms.
The average ISN He deflection is shown to occur in the B-V plane. Between 110 and 160 au, the average deflection of secondary atoms is ∼12°.0 from the upwind direction, and its median is ∼8°.0. The incoming direction (248°.1, 14°.5) differs by 3°.5 from that estimated by Kubiak et al. (2016) at 150 au. The average sunward velocity component of He 1 at 150 au is ∼12.5 km s −1 , which is within 1.5 km s −1 from the Kubiak et al. (2016) estimate of ∼11 km s −1 . The largest discrepancy from the previous publications is found in the temperature. The total He 1 atoms temperature at 150 au is He ,He ,He 1 1 1 ≈ 25,500 K, significantly larger than the WB estimates of ∼9500 K by Kubiak et al. (2016) and closer to the earlier estimate of 15,000 -+ 8000 6000 K (Kubiak et al. 2014). The temperature anisotropy of secondary ions, T P /T ⊥ ≈ 3.5, is 2.2 times the best fit of IBEX data obtained by Wood et al. (2019). Such a high temperature is due the atoms that originate in different regions of the OHS and arrive at any specific location from different directions, and it is not an artifact of the relatively large volumes used for collecting statistics. This kinetic effect is related to the large mean free paths (∼1500 au) of ISN He, which allow the neutral atom distributions to be significantly different from those of the plasma (Figure 2). It is of interest that even in the case of an initially isotropic, Maxwellian helium distribution in the LISM and Maxwellian OHS plasma, the distribution of the secondary He is not Maxwellian and becomes fully anisotropic. The distribution of primary He remains close to Maxwellian in a narrow range from the peak, and the total temperature remains essentially unchanged, at least up to 10 au from the Sun. However, their distribution in the SW is anisotropic, exhibiting parallel cooling and perpendicular heating. Moreover, all distributions become asymmetric with decreasing distance. Due to the small volume available for collecting statistics, further improvements will be necessary to obtain the VDFs at 1 au and strictly compare our simulation with IBEX observations. Finally, we emphasize that these kinetic effects are entirely related to charge exchange as a quantum mechanical process involving indistinguishable particles with the cross section strongly peaked at zero pitch angle. This makes charge exchange of the like particles indistinguishable from head-on collisions (McNutt et al. 1998). This process is dominant from the viewpoint of the SW-LISM interaction, since it gives birth to energetic WB. Swaczyna et al. (2021) have recently shown that elastic collisions between ISN He and protons may increase the antisunward parallel temperature of the primary ISN He by ∼1000 K, and reduce its speed by ∼0.45 km s −1 . Interestingly, the former effect is opposite to the parallel cooling described here. The the semiclassical JWKB approach of Swaczyna et al. (2021) does not include the charge exchange and therefore misses the filtration process of He atoms described in this Letter. Since both elastic and charge-exchange collisions are important for the interpretation of IBEX data, these results complement each other. In summary, it ought to be concluded that large-scale simulations are necessary to infer the temperature and velocity of the unperturbed LISM from observations made at 1 au.