On the Energization of Pickup Ions Downstream of the Heliospheric Termination Shock by Comparing 0.52–55 keV Observed Energetic Neutral Atom Spectra to Ones Inferred from Proton Hybrid Simulations

We present an unprecedented comparison of ∼0.52–55 keV energetic neutral atom (ENA) heliosheath measurements, remotely sensed by the Interstellar Boundary Explorer (IBEX) mission and the Ion and Neutral Camera (INCA) on the Cassini mission, with modeled ENAs inferred from interstellar pickup protons that have been accelerated at the termination shock, using hybrid simulations, to assess the pickup ion energetics within the heliosheath. This is the first study to use hybrid simulations that are able to accurately model the acceleration of ions to tens of keV energies, which is essential in order to model ENA fluxes in the heliosheath, covering the full energy range observed by IBEX and CASSINI/INCA. The observed ENA intensities are an average value over the time period from 2009 to the end of 2012, along the Voyager 2 (V2) trajectory. The hybrid simulations upstream of the termination shock, where V2 crossed, are constrained by observations. We report an energy-dependent discrepancy between observed and simulated ENA fluxes, with the observed ENA fluxes being persistently higher than the simulated ones. Our analysis reveals that the termination shock may not accelerate pickup ions to sufficient energies to account for the observed ENA fluxes. We, thus, suggest that the further acceleration of these pickup ions is most likely occurring within the heliosheath, via additional physical processes like turbulence or magnetic reconnection. However, the redistribution of energy inside the heliosheath remains an open question.


Introduction
As the solar system and its surrounding heliosphere move through the local interstellar medium (LISM), interstellar neutral (ISN) atoms, mostly atomic hydrogen with densities of >0.12 /cm 3 (Dialynas et al. 2019;Swaczyna et al. 2020), enter the heliosphere and undergo charge-exchange collisions with the continuously flowing solar wind (SW) protons. During this process, SW protons gain an electron and become neutral hydrogen, still flowing outward at the SW velocity. Newly created ions from the ISN population are advected outward with the SW under the force of the V × B electric field, typically as a ring-beam distribution that is isotropized by scattering from self-excited and preexisting magnetic fluctuations, forming a population that is commonly known as pickup ions (PUIs). Recent observations from the New Horizons spacecraft at ∼38 au (McComas et al. 2017) showed that PUIs are heated in the frame of the SW with increasing distance (McComas et al. 2021), before reaching the termination shock (TS). At the shock, they are further heated, with a fraction of their distribution being reflected off the shock surface and undergoing additional heating, (Zank et al. 1996.
The two Voyager spacecraft (V1 and V2) reached the TS in 2004 and 2007 at distances of ∼94 (Decker et al. 2005;Stone et al. 2005) and ∼84 au (Decker et al. 2008), respectively. Voyager observations revealed that the TS was mediated by PUIs (in agreement with the theoretical prediction of Zank et al. 1996Zank et al. , 2010, where roughly 80% of the solar wind flow energy was transferred to PUIs in the heliosheath (HS), the region between the TS and the heliopause (i.e., the interface between our heliosphere and the very local interstellar medium (VLISM)), including a substantial part (>15%) that went into >28 keV protons (Decker et al. 2008;Richardson et al. 2008). The heliopause was observed at distances of ∼122 au Gurnett et al. 2013;Krimigis et al. 2013;Stone et al. 2013) and ∼119  The heated PUIs that populate the HS charge-exchange with the interstellar neutrals and are measured remotely by the Interstellar Boundary Explorer (IBEX; 0.01-6 keV) and Cassini/Ion and Neutral Camera (INCA; 5.2-55 keV), providing energetic neutral atom (ENA) full-sky maps Krimigis et al. 2009). These images showed the existence of a bright and narrow ribbon (Schwadron et al. 2009) of ENA emissions that is measured by IBEX-Hi and is thought to lie beyond the HP, formed through a secondary ENA process (e.g., Heerikhuisen et al. 2010;McComas et al. 2017) and the globally distributed flux (GDF; Schwadron et al. 2011), a "background" ENA flux, with the IBEX ribbon, which is largely produced in the HS, removed (Dayeh et al. 2011;Zank et al. 2010;McComas et al. 2020;Zirnstein et al. 2020). Spectral agreements between Cassini/ INCA ENA measurements and the Voyager ion measurements suggest that the source of ENA emissions at energies above 5.2 keV are most likely produced by charge-exchange interactions inside the HS (Dialynas et al. 2013(Dialynas et al. , 2017. Therefore, understanding the PUI distribution downstream of the TS is essential to study the pressure balance and acceleration mechanisms inside the HS (Dialynas et al. 2019(Dialynas et al. , 2020. This understanding is needed to determine the emission of ENAs from the HS because these ENAs are used to remotely sense the boundaries of our heliosphere and its interaction with the VLISM. This letter provides an unprecedented comparison between observed ENAs in the HS (Section 3), combining IBEX-Hi measurements of energies from 0.52 to 6 keV with higherenergy ENAs from Cassini/INCA from 5.2 to 55 keV, and modeled ENA spectra inferred from hybrid simulations of interstellar pickup protons accelerated at the TS using realistic parameters (Section 2). Previous studies have investigated the PUI heating downstream of the TS using self-consistent fluid models (Fahr & Chalov 2008;Fahr & Siewert 2010;Zieger et al. 2015;Opher et al. 2020), or particle-in-cell (PIC) simulations (Lembège and Yang 2018). Global fluid simulations only capture adiabatic processes and not kinetic processes occurring at the TS. On the other hand, PIC simulations are constrained to provide accurate results only to lower energies due to the size of numerical boxes. This is the first study incorporating a self-consistent hybrid simulation, which is able to accurately predict the acceleration of PUIs downstream of the TS to energies up to tens of keV. Our comparison shows that the observed ENA fluxes are higher than the modeled ones throughout the 0.52-55 keV energy range and that the identified discrepancy is energy-dependent. We conclude (Section 4) with the suggestion that there is an acceleration region for these ions that lies within the HS. We finally discuss additional physical processes that could result in energy redistribution within the HS, which could guide future modeling efforts.

Simulated PUI Distribution
We use the hybrid model by Giacalone et al. (2021) to produce interstellar pickup protons downstream of the TS, which have been reflected and some further accelerated at the TS (Zank et al. 1996). The hybrid simulation is twodimensional, however, vectors, such as average proton velocity, magnetic field, electric field, etc., are not confined to the simulation plane, and they point in three directions. The initial magnetic field and bulk plasma velocity consist of an average component and a turbulent component. The simulation is a self-consistent, kinetic treatment of SW protons and pickup protons, as well as massless, charge-neutralizing fluid SW electrons. Giacalone et al. (2021) determined the intensity of accelerated interstellar PUIs at three different locations: the V2 crossing, the flank, and the tail. For the purposes of this study, we use the results at the V2-crossing location, where simulation parameters are based on V2 observations just before the TS crossing (see the second column of Table 1 in Giacalone et al. 2021). Figure 1 shows energy spectra from the hybrid simulation for two proton populations, the SW (green dotted lines) and PUI (purple dotted lines), in the shock frame. These spectra are averaged over nearly the entire downstream region of the simulation. The simulations in Figure 1 show high-energy tails that are formed for both the SW (starting at ∼3 keV) and PUI (starting at ∼5 keV) distributions, with the intensity of accelerated SW protons being significantly lower than that of the PUIs. The simulated PUI intensity decreases significantly above ∼50 keV, due to artificial losses arising from the small simulation domain and run time (Giacalone et al. 2021). As a result, there is a substantial mismatch between the simulated spectra and the ion measurements from the Low Energy Charged Particle (LECP) detector on V2 (e.g., Decker et al. 2008;Dialynas et al. 2019). It should be noted, however, that the simulated PUI tail intensities in Giacalone et al. (2021) right before the abrupt decrease, i.e., up to ∼80 keV, agree well with the corresponding LECP intensities as measured downstream of the TS.
We then use these PUI intensities produced by the hybrid model downstream of the TS to infer ENA intensities throughout the HS. We use the ENA model described in Kornbleuth et al. (2021a), where bulk plasma velocity streamlines are extracted from the BU MHD model simulations, described in Section 3, in order to simulate ions crossing the TS and transiting through the heliosphere until they chargeexchange with interstellar neutrals and become ENAs.
In the ENA model, similarly to the methodology followed in Zank et al. (2010) and Zirnstein et al. (2017), PUIs and SW protons are treated as one fluid, and this requires partitioning of the plasma energy among three different ion populations, namely, (i) thermal SW protons, (ii) PUIs created in the supersonic SW and being adiabatically transmitted across the TS, and (iii) PUIs being reflected at the TS until they have sufficient energy to overcome the cross-shock potential barrier (Zank et al. 1996). We partition the total thermal energy of the plasma via where n p and T p are the density and temperature of the plasma, n i is the density of the respective ion population, and Γ i is the temperature fraction, T i /T p , where T i is the temperature of the respective ion population. We model downstream plasma conditions in the V2 direction by fitting the ion spectra produced by the hybrid simulation with a Maxwellian distribution in the case of cold SW protons, with n Sw = 0.74 * n p and thermal energy E SW = 0.07 * E p , and for the transmitted (energized) PUIs with parameter κ = 10 (1.57), density = , where n p is the plasma density and E p is the thermal energy of the plasma extracted from the BU MHD model of Kornbleuth et al. (2021b). These fits are depicted with dashed lines in Figure 1 (green for SW, blue for transmitted PUIs, and red for the reflected and energized PUIs). The purple solid line shows the total PUI population spectra, combining the fits of the transmitted and energized populations. By fitting these distribution functions, we are able to determine the density and temperature ratios for each population downstream of the TS, which is used in ENA modeling, as described. The use of the hybrid simulation from Giacalone et al. (2021) ensures a reasonable agreement of the simulated fluxes with the measured LECP ions (within a limited energy range of ∼28 to <80 keV), while the separate treatment of the thermal solar wind protons and PUIs allows us to constrain high-energy ion properties downstream of the TS. In our simulations we assume quasi-neutrality, i.e., = = + + n n n n n e p Sw PUI PUI trans energ , and that electrons and solar wind protons have the same temperature. It should be noted that the latter assumption, used in all global models of the heliosphere, is not necessarily accurate. Theoretical studies have shown that, in principle, electrons can be hot ahead of the TS, with the same temperature as the PUIs (Chalov & Fahr 2013;Chashei & Fahr 2013. Such observations have not been possible yet because of a measurement gap between ∼6 keV and ∼30 keV in measuring plasma properties on board of Voyager spacecraft (Richardson & Decker 2014). If that is the case, electrons can gain substantial energy across the TS (Zieger et al. 2015;Fahr & Siewert 2015;.  axis (Funsten et al. 2009, see also   These observed ENA fluxes are directly contrasted with simulated ENA fluxes produced by the same ENA model used in Kornbleuth et al. (2021a). More specifically, we fit the PUI densities and temperatures just downstream of the TS from the BU MHD simulations with the hybrid results, as described in Section 2, based on ratios relative to plasma density and temperature. Then, those PUIs are transported farther into the HS, following bulk plasma velocity streamlines, extracted from the BU MHD model simulations, until they charge-exchange with interstellar neutrals and become ENAs.

Comparison of Modeled and Observed ENA Spectra
At the outer boundary of the BU MHD simulation (1500 au from the Sun), we assume an interstellar proton density of n p,LISM = 0.04 cm −3 and interstellar neutral hydrogen density of n H,LISM = 0.14 cm −3 (lower than those reported by Swaczyna et al. 2020, i.e., ∼0.19 cm −3 ), based on the study of Izmodenov & Alexashov (2015) who used these parameters to best fit the TS distance with respect to Voyager measurements. The neutral and ionized populations in the interstellar medium are assumed to have the same bulk velocity and temperature at the outer boundary in the pristine ISM, given by v ISM = 26.4 km s −1 (longitude = 75°.4, latitude = −5°.2 in the ecliptic J2000 coordinate system) and T ISM = 6530 K, respectively. Based on the work of Izmodenov & Alexashov (2020), the BU MHD simulation uses 22 yr averaged solar cycle conditions from the years 1995 to 2017, with heliolatitudinal variations of the SW speed and density based on interplanetary scintillation data (Tokumaru et al. 2012) and SOHO/SWAM full-sky maps of backscattered Lyα intensities (QuéMerais et al. 2006;Lallement et al. 2010;Katushkina et al. 2013Katushkina et al. , 2019. The proton distributions predicted by the hybrid model (SW and transmitted and energized PUIs) are transported within the HS following the MHD streamlines, even though that assumption might not be exactly accurate for the most energetic ones, and they charge-exchange with interstellar neutral hydrogen producing ENAs. The ENA flux along a radial line of sight (LOS) is given by where ¢ r is the vector along a particular LOS as a function of θ and j and s is the distance along the vector, m p is the mass of a proton, and f p is the phase-space velocity distribution, which is treated as a Maxwellian for SW protons and as a kappa for PUIs, as described in Section 2. The TS of the MHD model is at 82 ± 4 au, and the heliopause is at 126 ± 1 au. The integral only captures 44 au in the HS along the V2 LOS, which is larger than the width of the HS as measured by V2, which crossed the TS at ∼84 and the heliopause at ∼119 au. The velocity of the protons in the frame of the plasma is given by where v p and v i are the velocities of the bulk plasma and the parent proton, respectively. For the density and temperature of the given proton population, we use n i and T i , respectively, as described in the previous section. ¢ ( ( )) r n s H is the neutral H density along the LOS, s ( ) E ex is the chargeexchange cross section from Lindsay & Stebbings (2005), and ( ) S E is the survival probability, which represents the likelihood that an ENA created in the HS will charge-exchange prior to being observed at 1 au. The survival probability of an ENA is calculated along the radial LOS as where dr is the radial element over which we are integrating, υ ENA is the speed of the ENA, and υ rel is the relative velocity between the ENA and the bulk plasma given by Pauls et al. (1995), where v ENA is the velocity of the ENA, u p is the bulk-averaged plasma velocity, and υ th,p is the thermal speed of the plasma. It should be noted that Equation (4) assumes a Maxwellian distribution of protons, even though a kappa distribution, with a nonthermal tail, could have been a better approximation. Nonetheless, the relative velocity between the ENAs and the plasma is sufficiently large in the realm where the kappa is applied, such that the Maxwellian assumption has only a small impact. Heerikhuisen et al. (2019) showed that, although assuming a kappa distribution for the HS protons, compared to Maxwellian one, leads to different plasma properties in the heliotail, including the energy-dependent charge-exchange cross section into the collision integrals reduces the magnitude of these differences. In calculating the survival probability, we assume the observer to be located at 100 au because the ENA observations are corrected for survivability from 100 au to the inner solar system. The modeled and observed ENA fluxes in Figure 2(a) show a clear energy-dependent discrepancy throughout the 0.52-55 keV energy range (blue shaded area), with the observations being consistently higher. The ratios of the observed ENA fluxes over the modeled ones, as a function of energy, are shown in Figure 2(b) (and embedded table). Those ratios increase with energy for ENAs with energies from 0.71 and higher, peaking at about 8.38 keV, and it starts decreasing with energy beyond this point. The modeled ENA spectra also exhibit different slopes than the measured ENA spectra throughout the 0.52-55 keV energy range. The resulting spectra from the model are softer (harder) than the observed ones in the 0.52-8 keV energy range (beyond ∼8 keV), i.e., the ratio between the measurements and the model within that energy range is increasing (decreasing). The overall best comparison between the model and the measurements occurs for >18 keV, which is not surprising considering that (i) the PUIs from the hybrid simulation that we use (Giacalone et al. 2021) agree well with the ∼28 keV measurements from LECP in V2, and (ii) the conversions of the in situ ∼28 to 540 keV LECP ions to ENAs (Dialynas et al. 2020) result in ENA spectra that retain power-law slopes similar to those of the measured 5.2-55 keV ENA spectra from Cassini/INCA (Krimigis et al. 2009;Dialynas et al. 2020).

Discussion and Conclusions
We have compared observed ENAs from IBEX-Hi (0.52-6 keV) and Cassini/INCA (5.2-55 keV) with modeled ENA spectra inferred from interstellar pickup protons that have been accelerated at the TS using hybrid simulations with realistic parameters. Our comparison shows that there is an energydependent discrepancy between the modeled and observed ENAs throughout the 0.52-55 keV energy range, with the observed ENAs being consistently higher than the modeled ones. The disagreement with observations potentially indicates that an additional physical process, energizing further the PUIs, has yet to be identified. Such a process would also require a source of energy. In the hybrid simulation, energy is conserved, and it is the upstream solar wind ram energy that is accelerating the ions. Even though the hybrid model is local, it captures the main components of turbulence (that is, <0.01 of the ram pressure), so the distribution of energy across the TS is captured correctly. How the redistribution of energies farther downstream in the HS (including SW ram pressure, thermal pressure, etc.), as well as any other source of energy yet to be identified, could further accelerate the PUIs is still an open question. Chalov et al. (2003) pointed out the potential importance of high levels of turbulence in the inner HS in order to produce observed fluxes of hydrogen ENAs > 60 keV and argued that combined observations of ENA fluxes at energies near 1 keV and at energies of several tens of keV are necessary to constrain the parameters of the SW turbulence at the TS and inner HS. Interaction of PUIs with turbulence both before they cross the TS and in the HS can lead to their stochastic acceleration due to wave-particle interactions (Chalov & Fahr 2000;Chalov et al. 2007). Furthermore, interplanetary shocks transmitted through the TS have been shown to result in further heating of the PUIs in the HS (Mostafavi et al. 2019). Nonetheless, it should be noted that for ENAs of energies <1 keV, an alternative possible reason for the data-model discrepancy could be the source of these ENAs not being the HS but the VLISM, as discussed in Fuselier et al. (2021).
Using test-particle simulations, Zirnstein et al. (2021) showed that although a moderate level of turbulence with a power ratio of @ D ( ) 0.011

B B
2 o (based on V2/MAG measurements from Burlaga et al. 2008 at the TS) was sufficient to reproduce a superthermal PUI tail downstream of the TS, the derived intensities were lower than those observed by IBEX-Hi. In fact, to produce a proton distribution consistent with IBEX observations, a 10 times larger turbulence power ratio must be applied at the shock foot. The authors speculated that this enhanced level of turbulence in the PUI foot may be related to shock self-re-formation (Lembège et al. 2004), shock-front ripples (Umeda et al. 2014) originating from instabilities induced by ion temperature anisotropies (Winske & Quest1988), cross-field currents (Lembège and Savoini 1992), or interplanetary SW turbulence transported through the TS (Giacalone et al. 2021;Zank et al. 2021).
In this study, we have used a wider range of ENA energies, including also higher-energy ENAs from the Cassini/INCA instrument (∼5.2-55 keV), and we also used a self-consistent hybrid simulation to produce the PUI population in the HS. Similarly to Zirnstein et al. (2021), we find that even selfconsistent hybrid simulations do not accelerate PUIs enough at the TS in order to agree with the observed spectra. Furthermore, we find that the data-model discrepancy exhibits a clear energy dependence, where the modeled ENA spectra do not capture the observed spectral slopes. Thus, we infer that an additional acceleration region for these low-energy PUIs is most likely the HS.
Models have shown that processes like turbulence or magnetic reconnection (Drake et al. 2010;Opher et al. 2011) may occur throughout the HS, as opposed to only close to the TS. Zhao et al. (2019) analyzed V2 data and found evidence of magnetic flux ropes/islands accompanied by energetic protons, suggesting reconnection within the HS. The source of that HS turbulence has been investigated in several recent studies. Zieger et al. (2020) found that dispersive fast magnetosonic waves can be responsible for energy dissipation downstream of the TS. Zank et al. (2018) explored the transmission of nearly incompressible turbulence across the TS. Compressive turbulence within the HS may also further heat the PUIs after crossing the TS (Fisk & Gloeckler 2017). Zank et al. (2021) studied the interaction and transmission of quasi-2D turbulence through a collisionless perpendicular shock wave in the large beta regime and found that the downstream spectral amplitude is increased significantly. Simulations also demonstrate that a Rayleigh-Taylor-like instability (Opher et al. 2021), which mixes the HS and LISM plasmas, can form turbulent heliospheric jets with scales of order 100 au and a turnover timescale of years . How far downstream the turbulence ensues, whether the above processes could further accelerate the PUIs in sufficient numbers, depending on the source of energy, and how these accelerated PUIs get transported further into the HS should be investigated in future studies.
Although exploring each of the processes mentioned above, which could potentially act to further accelerate PUIs within the HS, lies beyond the scope of the present Letter, we anticipate that the reported discrepancy will guide relevant future modeling efforts. Such investigations will be especially crucial when data from the upcoming Interstellar Mapping and Acceleration Probe mission (IMAP; McComas et al. 2018) will become available, covering ENAs of energies from 0.005-300 keV taken from instruments on the same platform in L1 orbit, offering major advancements compared to previous observations, such as a significant increase in the collection power of all the ENA cameras on board.