Heliosheath Proton Distribution in the Plasma Reference Frame

Properties of the inner heliosheath (IHS) plasma are inferred from energetic neutral atom (ENA) observations by ∼1 au spacecraft. However, the Compton–Getting effect due to the plasma velocity relative to the spacecraft is rarely taken into account, even though the plasma speed is a significant fraction of the ENA speed. In this study, we transform Interstellar Boundary Explorer (IBEX) ENA spectra to the IHS plasma frame using flow profiles from a 3D heliosphere simulation. We find that proton spectra in the plasma frame are steeper by ∼30% to 5% at ∼0.5 to 6 keV, respectively, compared to ENAs in the spacecraft frame. While radial plasma flows contribute most to the Compton–Getting effect, transverse flows at mid/high latitudes and the heliosphere flanks account for up to ∼30% of the frame transformation for IBEX-Hi at ∼0.7 keV and up to ∼60% for IBEX-Lo at ∼0.1 keV. We determine that the majority of IHS proton fluxes derived from IBEX-Hi measurements in 2009–2016 are statistically consistent with power-law distributions, with mean proton index ∼2.1 and standard deviation ∼0.4. We find significantly fewer spectral breaks in IBEX observations compared to early analyses, which we determine were a product of the “ion gun” background prevalent in ∼2009–2012 before corrections made by McComas et al. in subsequent data releases. We recommend that future analyses of the IHS plasma utilizing ENA measurements take into account the Compton–Getting effect including radial and transverse flows, particularly IBEX and Interstellar Mapping and Acceleration Probe measurements below ∼10 keV.


Introduction
The interaction of the solar wind (SW) with the very local interstellar medium (VLISM) creates the heliosphere as the solar system moves through the interstellar medium. The heliosphere is a region of SW plasma that is separated from the interstellar plasma by a tangential discontinuity known as the heliopause. Because the two plasmas do not mix, their primary interaction is through charge exchange between ions and neutral atoms (e.g., Zank 2015). Inside of the heliopause, the SW emitted from the Sun encounters a termination shock (TS), becoming subsonic in a region known as the inner heliosheath (IHS). The plasma in the IHS has relatively high pressure that provides pressure balance (in a stationary system) on the heliopause, and as it flows away from the TS, it is diverted around the heliopause toward the heliotail. The IHS plasma consists of highly energetic protons (1 keV) whose speeds in the solar inertial frame greatly exceed the bulk plasma speed away from the Sun. If such protons capture an electron during a charge-exchange collision with a passing neutral atom, the resulting hydrogen atom may ballistically travel toward the Sun and can be detected by near-Earth spacecraft. The observation of these "energetic neutral atoms" (ENAs) near Earth allows us to remotely probe the outer heliosphere, understanding its size, structure, and the thermodynamic properties of the IHS plasma and VLISM.
The Interstellar Boundary Explorer (IBEX), a mission dedicated to measuring outer heliospheric ENAs and probing the SW-VLISM interaction (McComas et al. 2009a), has recently completed imaging the outer heliosphere for a full solar cycle (solar cycle 24), revealing the dynamic behavior of the outer heliosphere . Heliospheric ENA observations can be divided into two broad sources: (1) the globally distributed flux (GDF; McComas et al. 2009b;Schwadron et al. 2011) encompassing the entire sky with broad enhancements and depletions in different directions, and (2) the narrow "ribbon" of enhanced ENA emission that encircles the sky (Funsten et al. 2009b;Fuselier et al. 2009;McComas et al. 2009b). The measurements at the highest energies made by IBEX-Hi (Funsten et al. 2009a) observe ENAs in the spacecraft frame of reference between ∼0.5 and 6 keV FWHM. The source of the GDF is believed to be mostly from the IHS where the suprathermal pickup ion (PUI) component was preferentially heated at the TS (e.g., Zank et al. 1996Zank et al. , 2010Richardson et al. 2008;Yang et al. 2015;Kumar et al. 2018) and dominates the ∼1 keV ENA GDF signal. The ribbon, on the other hand, originates from secondary ENAs from outside the heliopause (McComas et al. 2009b;Zirnstein et al. 2015a), originally sourced from the neutralized supersonic SW that propagated outside the heliosphere and ionized in the relatively dense outer heliosheath (OHS) plasma. It is clear that the GDF and ribbon are two, mostly distinct ENA sources (note that part of the GDF may originate from outside the heliopause; Izmodenov et al. 2009;Opher et al. 2013;Desai et al. 2014;Zirnstein et al. 2014), which has been all but proven by the differences in their evolution over time , with confirmation from time-dependent ribbon models (Zirnstein et al. 2015b(Zirnstein et al. , 2020c. Analyses of the GDF signal, which is the focus of our study, usually concentrate on either the unique temporal evolution of the GDF as a function of energy and direction in the sky (e.g., McComas et al. 2012McComas et al. , 2014McComas et al. , 2017aMcComas et al. , 2020 and references therein), or narrow their attention to the spectra index (e.g., Dayeh et al. 2011Dayeh et al. , 2012Dayeh et al. , 2014Desai et al. 2015Desai et al. , 2016Desai et al. , 2019 and its thermodynamic properties or temporal evolution. While the GDF spectral evolution is an important and topical subject, the fundamental source of the GDF in the IHS is of equal importance, as it informs us of the bulk plasma flows and pressures within the IHS region. For example, Livadiotis et al. (2011 used IBEX-Hi ENA observations to derive the IHS plasma thermodynamic properties, under the assumption that the spectra can be represented by a kappa distribution. These authors derived the spatial distribution of the plasma temperature (∼1 × 10 6 K) and kappa index (∼1.7), as well as determining that only the low-latitude ENA flux may originate from a single kappa distribution due to breaks in the ENA spectra at high latitudes (see Dayeh et al. 2012).
There is, however, an important effect that has not yet been fully taken into account in studies of the plasma source of the GDF ENA signal, which is that the IHS plasma is moving away, or in some cases, transverse, to the IBEX line of sight (LOS) at speeds ranging from ∼40 to 200 km s −1 , based on Voyager observations (Krimigis et al. 2011;Richardson & Decker 2014), and likely faster speeds at higher latitudes during solar minimum (McComas et al. 2000;Pogorelov et al. 2013). Because these flow speeds are comparable to the speeds at which ENAs strike the IBEX detector, the Compton-Getting (CG) effect (Compton & Getting 1935; see also Gleeson & Axford 1968;Forman 1970;Ipavich 1974) will greatly influence the intensity and spectral slope of the ENA and parent proton flux derived in the spacecraft reference frame Zirnstein et al. 2013). Studies that separated the ribbon signal from the GDF (e.g., Schwadron et al. 2014) did include an estimate for the IHS plasma flow speed to derive the partial IHS plasma pressure but only using a single value for the radial IHS flow speed in the entire sky. Because the IHS plasma flow speed may be a significant fraction of the ENA speeds observed by IBEX-Hi (a plasma flow speed of ∼250 km s −1 is ∼80% the speed of 0.5 keV ENAs observed at the low end of IBEX-Hiʼs energy passband 2, with central energy ∼0.7 keV), and the IHS plasma flow changes in the sky, we utilize simulations to help determine how a more realistic transformation may affect inferences of the IHS plasma properties.
It is the goal of this study to present a method for transforming IBEX ENA observations (the methods apply to both IBEX-Lo and IBEX-Hi, but we specifically analyze IBEX-Hi observations in this study) to the IHS plasma reference frame and derive proton fluxes from IBEX data. This requires estimations of the LOS-averaged IHS plasma and neutral properties, which we extract from a 3D simulation of the SW-VLISM interaction. After deriving the proton fluxes in the IHS plasma frame, we calculate the proton spectral indices between each energy and analyze their properties and discuss the implications of our results on our understanding of the IHS and prior analyses of GDF ENA fluxes.

IBEX Observations
We analyze IBEX observations of heliospheric ENAs collected from 2009 through 2016, before the heliospheric response to the SW dynamic pressure increase in late 2014 , and specifically utilize measurements in the spacecraft reference frame-this allows us to directly transform from the spacecraft frame to the IHS plasma frame without introducing any small but finite uncertainties from the CG corrections made to produce data sets in the solar inertial frame (e.g., McComas et al. 2010). The data are survival-probability corrected, meaning that the ENA fluxes observed at 1 au are corrected for the losses of ENAs by ionization as they travel from their source near or beyond 100 au to 1 au, providing a better representation of the ENA source spectrum, albeit with additional uncertainties from this correction. Finally, we restrict our analysis to data collected in IBEXʼs ram frame, when its motion has a component toward its look direction, because the ram data have a higher signal-tonoise ratio compared to the anti-ram data. For more information on the ram, survival-probability-corrected IBEX-Hi data set, see McComas et al. (2020) and references therein.
A subset of the IBEX-Hi observations is shown in Figure 1. We present sky maps taken in IBEXʼs ram frame in 2009, 2013, and 2016 for IBEX-Hi energy passbands 2-6 (Funsten et al. 2009a), which have central energies 0.71, 1.11, 1.74, 2.73, and 4.29 keV. Note that we removed the ribbondominated regions of the sky from our analysis, because it is widely believed that the ribbon does not originate from the IHS and so its spectra cannot be included in our analysis of the IHS plasma. The ribbon is masked out by finding all pixels in the sky that are within 20°of a circle with center (218°.33 40°.38) in ecliptic J2000 coordinates and radius 74°.81 . We note that while the ribbon's center and width do change with ENA energy (e.g., Funsten et al. 2009bFunsten et al. , 2013Fuselier et al. 2009;Dayeh et al. 2019), it is important to keep pixels consistently over all energies. Therefore, we use the ribbon center and radius averaged over all energies and time from Dayeh et al. (2019) with a width of 40°to mask out the ribbon consistently over all energies. A width of 40°is chosen such that it masks out most of the ribbon observed in energy passband 6 (∼4.3 keV), which appears wider than that at lower energies (e.g., Schwadron et al. 2014). Moreover, any pixels whose uncertainty is greater than the flux or have gaps due to, e.g., high background are also removed (representing ∼4% of the data set).
where θ and f are the angles in spherical coordinates, v is the speed, the prime variables are in the plasma frame, and the tilde variables are in the spacecraft frame. The velocity vectors in these two reference frames are related as ¢ =v v u rel  , where u rel = u p − u sc is the relative velocity between the plasma frame (u p ) and spacecraft frame (u sc ). The change in flux is  1.11, 1.74, 2.73, and 4.29 keV) in 2009, 2013, and 2016 in the spacecraft ram frame. Ribbon fluxes and any data points with uncertainty greater than the flux are removed (gray pixels). Directions toward the nose, Voyager 1 ("V1"), Voyager 2 ("V2"), port ("P"), and starboard ("S") directions are also shown.
given by the ratio where α is the angle between the ENA velocity in the spacecraft frame, v , and u rel . The frames are connected as where u p is the bulk plasma velocity in the solar inertial frame extracted from our 3D simulation of the heliosheath, and u sc is the IBEX spacecraft velocity (where u sc = 30 km s −1 ) in the solar inertial frame (as we perform our calculations in ecliptic J2000, u sc,z = 0). The IBEX-Hi boresight follows a great circle perpendicular to the Sun-pointing spacecraft axis. For this study, we analyze only the ram data in which the spacecraft is moving toward the incoming atoms. Consequently, the spacecraft velocity is opposite to the projection of ENAs' velocity on the ecliptic plane. An illustration of the frames of reference is shown in Figure 2. Next, the rate of charge exchange between H + and neutral H atoms is energy dependent, and we must convert the ENA fluxes in the spacecraft frame, J ṽ (˜), into proton fluxes in the plasma frame, , using the following equations as derived in Appendix A (see also Zirnstein et al. 2013 is the particle speed in the solar inertial frame, ¢ ¢ f v p ( ) is the proton distribution in the plasma frame, σ ex is the energy-dependent, charge-exchange cross section given by Lindsay & Stebbings (2005), and v rel is the relative speed of interaction between the ENA (or parent proton) at speed v′ in the plasma frame and the background neutral H population moving with bulk speed u H . Recent studies of the H-H + charge-exchange cross section show there is considerable uncertainty in using the cross section from Lindsay & Stebbings (2005) at low interaction energies (0.1 keV; Swaczyna et al. 2019;Bzowski & Heerikhuisen 2020). The cross section from Lindsay & Stebbings (2005), however, is sufficient at the energies we analyze. All variables on the right-hand side of Equation (4) are solved in the frame-transformed direction Ω. Note that the differential time Δt′ = Δt = Δl/v is invariant to frame transformations (for nonrelativistic particles) and therefore can be written in terms of distance Δl and speed v in the solar inertial frame. The relative speed v rel is given by (e.g., Zirnstein & McComas 2015; note, however, the typo in their Equation (6), fixed here) Figure 2. Illustration of the reference frames used in this study. Our analysis transforms ENA fluxes (ENA velocity v , black arrow) measured in IBEXʼs ram reference frame (spacecraft velocity u SC , gray arrow) to the IHS plasma frame ( ¢ v , red arrow) using estimates of plasma flow velocities from a 3D MHD simulation of the heliosphere. The angle α between the ENA velocity in the spacecraft frame, v , and the IHS plasma flow velocity relative to the spacecraft frame, u rel , is also shown (not to scale).
where ω is the relative speed between the ENA/proton and the background neutral H population (transformed to the same reference frame, which we define here in the solar inertial frame) scaled by the background H thermal speed = v k T m 2 th,H 2 B H H . We note that while the proton speed v ? u H at the energies we consider, in general, Equation (4) should be solved in the frame of the background neutral H population as we do in this study. This is especially important when computing low-energy ENA fluxes (0.1 keV; e.g., Zirnstein et al. 2013).
Equation (4) can be written as an average over the LOS, yielding i.e., the covariance of the model proton flux with respect to G along an LOS through the IHS, normalized by the total ENA flux it produces, is similar to the same quantity in the real heliosphere. This does not require the model to reproduce the observed ENA flux, but rather that the proton flux distributed as a function of distance l along the LOS, normalized by the total ENA flux it produces, is similar to reality. Then, Equation (9) can be simplified as is the model ENA flux in the spacecraft frame, and á ¢ ¢ ñ J v p,m ( ) is the average model proton flux in the plasma frame. Note that Equation (10) is equivalent to the LOS-weighted mean of the term Notice that Equations (4)-(11) include the term v rel /v that accounts for the fact that the rate of charge exchange is proportional to the speed of the proton speed relative to the neutral H frame, v rel . This ratio is close to 1 at IBEX-Hi energies but deviates from 1 at low energies due to the relative speeds of the background neutral H population with respect to the parent proton. The same care must be taken when evaluating the charge-exchange cross section at relative speed v rel instead of v. We also emphasize that while the distance L is not known for each direction in the sky, L is not required to determine the spectral index of the proton flux because L is defined to be the same for all ENA energies along an LOS (it is not the ENA source length, but the total integrated distance from the TS to the heliopause). This also applies to the neutral H density: while our simulation may slightly underestimate the neutral H density in the IHS (see discussion in Section 2.3.4), this does not significantly impact the results of our study.
While we ignore the covariance terms in Equations (10) and (11), we have calculated the ratio ˜(˜) from our model (see Appendix B for the model of protons in the IHS) to see how large the covariance term is compared to the model ENA flux. We determine from our model that this term is negative and small (<15%) in the entire sky for all IBEX-Hi passbands (maximum of ∼15% at ∼0.7 keV and smaller for higher ENA energies). If the actual covariance ratio is similar to that derived in our model, then these terms can be ignored in Equations (10) and (11). Even if the actual covariance term deviates from the model predictions by a significant amount, e.g., 50%, then Equations (10) and (11) are still accurate to <10%.
An equivalent form of Equation (11) can be written for the observed ENA flux in the plasma frame (i.e., without converting ENA flux to proton flux, assuming that á ¢ ñ v v where the terms in parentheses ... ( ) are defined as weighted means in the same way as in Equation (11). Note also that Equation (12) can be written as a function of the weighted mean of the ratio ¢ v v 2 2 , but the ENA flux is the probability weight rather than the proton flux. Equation (12) also yields the energy transformation for each IBEX-Hi energy passband , which we use to transform the observed energies to the plasma frame. While the majority of the results from our analysis are derived from Equation (11), examples of ENA flux in the plasma frame derived from Equation (12) are shown in Figure 6.

Simulation of the Heliosphere
In order to transform ENA fluxes into the IHS plasma frame using Equations (11) and (12), information on the plasma flow properties must be assumed in some fashion. For this reason, we utilize a 3D MHD/kinetic simulation of the heliosphere to calculate the LOS-average plasma and neutral H flow properties for each IBEX direction in the sky. The simulation we use solves a coupled system of ideal MHD equations for the SW and interstellar plasmas and Boltzmann's equation for neutral H atoms (for more details, see Pogorelov et al. 2008;). The SW properties at 1 au are derived from a combination of measurements from the OMNI database at low latitudes and Ulysses measurements at high latitudes, because we aim to mimic SW conditions near the end of solar cycle 23 and the beginning of solar cycle 24 (∼2009-2011). We also test the results of our study using different SW boundary conditions, which is explained in more detail in Appendix D.

Solar Wind Boundary Conditions from Ulysses and OMNI
We choose a period of time to extract SW measurements and apply as boundary conditions for our simulation that approximately accounts for the delay between the SW propagation past 1 au and ENA measurement by IBEX. Because it takes ∼1-1.5 yr for SW with speed ∼450 km s −1 to travel from the Sun to the TS (ranging from ∼90 to 140 au from the Sun, based on Voyager observations and our prior simulations), a similar time for ENAs to travel back from its IHS source, and ∼1.5-4.7 yr for heliosheath plasma with flow speed ∼100 km s −1 to move ∼30-100 au (the typical distances over which >1 keV ENAs are created), we estimate that ∼2004-2009 is a reasonable time period to extract SW observations to be used as our simulation's inner boundary conditions.
At high latitudes, we utilize measurements from Ulysses taken during the third fast polar scan in 2007(McComas et al. 2008) and calculate the average SW properties scaled to 1 au, specifically plasma speed (743 km s −1 ), density (2.23 cm −3 ; including both protons and alphas), temperature (2.98 ×10 5 K), and radial magnetic field assuming a Parker spiral (34.7 μG) averaged at latitudes >|±37°|. The latitude boundary of |±37°| approximately separates the low-latitude streamer belt from the fast SW emitted from the polar coronal holes during Ulysses' third polar scan (Zirnstein et al. 2020b), which we use as the boundary separating our fast and slow SW boundary conditions.
At low latitudes (<|±37°|), we extract measurements from the OMNI database taken during 2004-2009, as described above. We note that this time period is only an approximation; while a longer time period would better correspond to the 2009-2016 IBEX data set that we analyze, it would only be accurate in a fully time-dependent simulation, which is beyond the scope of our study (see further discussion below). We determine the average SW plasma speed (449 km s −1 ), density (6.53 cm −3 ; including both protons and alphas), temperature (1.02 × 10 5 K), and radial magnetic field component assuming a Parker spiral (37.4 μG) over this time period and apply them at low latitudes. The low-and high-latitude SW values are scaled to 10 au under adiabatic expansion, the inner boundary of our simulation, and the simulation is iteratively run with the coupled neutral code until a quasi-steady state is achieved. Details about the coupling between the MHD and neutral codes can be found in Section 2.3.5.
We note that IBEX has performed a full solar cycle of imaging the heliosphere from 2009 through 2019 . Our choice to average SW boundary conditions from the period 2004-2009 approximately applies to the first half of IBEX observations (2009)(2010)(2011)(2012)(2013) when the ENAs were created in the outer heliosphere during solar-minimum-like conditions. Technically, a fully time-dependent simulation over a longer period of time that overlaps the entire IBEX epoch would be a better representation of the SW conditions that affected the IBEX measurements. However, truly timedependent simulation and time-dependent frame transformation of the IBEX ENA fluxes are currently beyond our capability, requiring us to account for the latitudinally asymmetric evolution of the high-latitude SW (e.g., Sokół et al. 2020). Rather than increasing the time period for averaging the SW properties, instead, we apply our frame transformation analysis of IBEX ENA observations with the SW properties as described above to IBEX observations from 2009-2016 for the majority of this study, a time before the large increase in SW dynamic pressure. We then test our results by analyzing IBEX observations over shorter time periods (2009-2013 versus 2014-2016) in Section 3.5, and by utilizing different heliosphere simulations with different SW boundary conditions in Appendix D.

VLISM Boundary Conditions from IBEX
The VLISM outer boundary conditions of our simulation are the same as those used in some of our previous works, primarily derived from IBEX-Lo and IBEX-Hi observations. We use values for the interstellar plasma speed (25.4 km s −1 ; McComas et al. 2015), density (0.0856 cm −3 ; Zirnstein et al. 2016b), temperature (7500 K; McComas et al. 2015), and total magnetic field magnitude (2.93 μG) and direction (227°.28, 34°.62) in ecliptic J2000 (Zirnstein et al. 2016b) implemented at an outer boundary distance of 1000 au from the Sun. The VLISM neutral atom properties are described in Section 2.3.4.

Heuristic Presence of Interstellar He +
Protons constitute most of the VLISM plasma (nearly 60% of the dynamic pressure and number density), He + a smaller amount (∼40% of the dynamic pressure or ∼10% of the number density), and higher mass particles have trace contributions. While our simulation does not strictly include interstellar He + as a separate population, our single-fluid MHD plasma does act as the total effective plasma density and pressure. In this sense, while maintaining the same pressure on the heliosphere, we note that the presence of interstellar He + can be heuristically included by scaling the MHD plasma density during charge-exchange events between protons and neutral H atoms, such that the MHD plasma acts as only protons during charge exchange (because only ∼60% of the dynamic pressure of the plasma is actually protons). We estimate the presence of interstellar neutral He by scaling the MHD plasma density, n MHD , at position r in the simulation by a constant fraction. The local proton density outside the heliopause in the OHS, r n p,OHS ( ), is given as We also include the presence of SW alpha particles (because our Ulysses-and OMNI-derived plasma densities contain both protons and alphas) by performing a similar type of scaling to the MHD plasma density inside the heliopause. The proton density inside the heliopause during charge exchange, r n p,SW ( ), is calculated as which is calculated using a relative density ratio for alpha particles to protons of 4%, estimated from Ulysses and Wind observations (McComas et al. 2008;Kasper et al. 2012).

Interstellar Neutral H Properties
The distribution of interstellar neutral H atoms is solved kinetically in our simulation (e.g., , where they are initialized as a Maxwell-Boltzmann distribution at the upstream outer boundary with the same bulk velocity and temperature as the interstellar plasma. In order to determine the interstellar neutral H density, we first note that the filtration of interstellar neutral H through the heliosphere is affected by the scaling of the MHD plasma density using Equations (13) and (14) during charge exchange. We tested different values for the interstellar neutral H density at the simulation's outer boundary to derive a density such that the simulation approximately reproduces the neutral H density at the upwind TS inferred from spacecraft measurements extrapolated from near Earth (∼0.09 cm −3 ; Bzowski et al. 2009). We find that in order to reproduce this value, we must take a density for the interstellar neutral H at the outer boundary (1000 au) of 0.11 cm −3 .
A recent analysis of New Horizons' SWAP observations of interstellar H + PUIs by Swaczyna et al. (2020) has revealed that the interstellar H density near the upwind TS is actually ∼40% higher than previously thought (∼0.13 cm −3 instead of ∼0.09 cm −3 ). This change is significant, as described by Swaczyna et al. (2020), because a higher neutral H density would (1) increase the number of PUIs produced by charge exchange in the SW and (2) increase the rate of production of ENAs in the IHS. Implementing a higher interstellar H density in our simulation is currently beyond our capabilities because this will require extensive testing to derive the proper VLISM neutral H density due to the complex filtration of H through the upwind heliosphere boundary regions. A high neutral H density would proportionally increase the charge-exchange rate at all energies; therefore, this change does significantly affect our proton spectral index results. However, it may change the charge-exchange momentum and energy source terms significantly enough to make a noticeable difference in the plasma flow properties. While the global heliosheath flow is not likely much different than in the currently used model, as our model reproduces a number of heliospheric observations (e.g., Zirnstein et al. 2016b), other boundary conditions of the model (e.g., interstellar plasma density) would require a respective adjustment with the new value of the interstellar H density. Therefore, this will be the subject of a future analysis.

Simulated Plasma Distribution Function during Charge Exchange
The coupling of the MHD and kinetic modules by charge exchange necessarily requires an assumption for the proton distribution function throughout the SW-VLISM interaction. In most regions of the heliosphere (supersonic SW, OHS, VLISM), the bulk of the plasma can be represented by a distribution close to, but not necessarily exactly, a Maxwell-Boltzmann distribution; thus, we assume the plasma distribution is Maxwell-Boltzmann in these regions. However, in the IHS, the plasma consists of a superposition of relatively cold SW ions and nonthermal PUIs that propagate with the bulk plasma flow (McComas et al. 2017b) and are preferentially accelerated at the TS (Zank et al. 1996Kumar et al. 2018), forming a suprathermal PUI distribution in the IHS from which ENAs are created. It is not clear whether the proton distribution is strictly a kappa distribution, but several studies have suggested that the distribution is likely similar to a kappa function/power law at IBEX-Hi energies in many directions of the sky (e.g., Livadiotis et al. 2011;Fuselier et al. 2014;Zirnstein et al. 2018b). Therefore, we assume the proton distribution in the IHS is a kappa distribution with kappa index = 2.0 during charge exchange. We note that we have tested different values for kappa index in the simulation (e.g., kappa index = 1.63; Decker et al. 2005) and its effects on our analysis in Appendix D, but the results of our study are not significantly affected by the choice of the proton kappa index (within a reasonable range of ∼1.6-3.0), because this does not significantly change the IHS plasma flow (Heerikhuisen et al. 2015). Zirnstein & McComas (2015) showed that, under the presence of charge exchange, the spectral index of a proton distribution fit with a kappa function over ∼0.1-10 keV does not change significantly with distance along a streamline through the IHS (see their Figure 8) if the initial proton distribution downstream of the TS is described by a kappa distribution with kappa index between ∼2 and 3, which appears to be a reasonable assumption as our results will show.
As the supersonic SW flows to the outer heliosphere, the accumulation of PUIs by photoionization and charge exchange creates a filled shell-like distribution of PUIs superposed on the Maxwell-Boltzmann-like core SW ion distribution (Elliott et al. 2016;McComas et al. 2017b). The number density of PUIs increases steadily until reaching ∼15%-30% of the total proton density upstream of the TS (likely to be even more with a larger interstellar neutral H density; Swaczyna et al. 2020), and therefore a more accurate description of the plasma distribution during charge exchange in our simulation would be to include the PUI-filled shell. While this can be included in a future study, it should not significantly affect our results because we are most concerned with the bulk plasma flow properties in the IHS, and the MHD simulation does account for the total energy exchange between protons and neutral H atoms throughout the SW-VLISM interaction.

Extracting LOS-weighted Quantities from the Simulation
Equations (11) and (12) show the LOS-weighted quantities in parentheses ...
( ) used to transform IBEX ENA fluxes in the spacecraft frame to proton fluxes in the plasma frame. However, we also present other LOS-weighted variables derived from our simulation in this study, such as the plasma flow velocity components shown in Figure 4 and the frametransformed energy, angle, and velocity shown in Figure 5. All variables presented in this study are calculated as an LOSweighted average by stepping along each IBEX LOS from the TS (l TS ) to the heliopause (l HP ) in each direction of the sky, q f W = , (˜)  , and weighting the local desired variable/quantity by the local proton flux produced by our simulation as done in Equation (11). In a more general form, for any quantity A i , its LOS-weighted value is calculated as )is the local proton flux in the plasma frame predicted by our model. A model of the proton distribution along each LOS must be used in order to perform the frame transformation of IBEX data. A description of the model proton distribution is described in Appendix B.
It is also important to note that IBEX ENA fluxes are functions of the instrument energy passbands, with log-Gaussian-like response functions (Funsten et al. 2009a). Therefore, the variables W A v , i (˜), such as the CG/chargeexchange conversion variable shown in Equations (11) and (12), are also integrated over the energy passband response functions. In the general form, Equation (15) is integrated over multiple energies per passband as here W E j (˜) is the energy response for energy passband ("PB" in Equation (16)) j (where j = 2, 3, 4, 5, or 6) at energy Ẽ in the spacecraft frame. The energy ranges for each passband span approximately twice the FWHM of the response function, which is effectively zero at their ends (see Funsten et al. 2009a). The integration of Equation (16) is effectively performed by calculating the LOS-weighted variables in each direction of the sky multiple times for a range of energies per IBEX energy passband (properly transformed to the plasma frame), weighted by the corresponding calibration response function.

Simulated Heliosphere and Effects of the Frame Transformation
Here we describe the results of our simulated heliosphere, beginning first with a comparison to Voyager observations of the heliosphere boundaries. The heliosphere simulation with SW and VLISM boundary conditions described in Sections 2.3.1-2.3.4 yields distances of ∼90 and 88 au from the Sun to the TS in the Voyager 1 and 2 directions, respectively, and ∼145 and 138 au from the Sun to the heliopause (producing IHS thicknesses of ∼55 and 50 au, respectively), which overestimate the distances to the heliopause from those observed at the times of the Voyagers' crossings (Gurnett et al. 2013;Stone et al. 2019). No simulations, in fact, can reproduce the thickness of the IHS observed by the Voyagers, but comparisons are better if timevarying SW conditions (not just latitudinal-varying conditions) are included and the heliosphere boundaries can move over time (Borovikov & Pogorelov 2014;Pogorelov et al. 2017;Washimi et al. 2017;Izmodenov & Alexashov 2020), or the IHS plasma flow properties are modified beyond ideal MHD (e.g., Malama et al. 2006;Izmodenov et al. 2014;Pogorelov et al. 2016;Heerikhuisen et al. 2019;Opher et al. 2020). However, our boundary distances compare just as well, if not better, than most heliosphere steady-state simulations.
The dynamical motion of the heliosphere boundaries during ∼2004-2012 (e.g., Pogorelov et al. 2017;Washimi et al. 2017) suggests that the average IHS thickness in the Voyager directions may actually be larger than the distances inferred from the instantaneous Voyager crossings. Note, however, that the boundary distances also depend on the VLISM conditions imposed in the simulation, which is why utilizing other observations to constrain the VLISM conditions (e.g., IBEX ribbon; Zirnstein et al. 2016b) is important. Because our simulation is in a steady state and does not account for the dynamic motion of the heliosphere boundaries, it is not surprising that it overestimates the thickness of the IHS and overestimates the distance to the heliopause by ∼20 au compared to the Voyagers' observations, at least during their times of crossing. Because we cannot realistically expect to reproduce both the distances to the TS and heliopause in a steady state, we can at least aim to reproduce one of them. As the goals of our analysis are to utilize plasma and neutral flow properties within the IHS to transform IBEX ENA fluxes into proton fluxes in the IHS plasma frame, we find it is better to approximately reproduce the distances to the TS rather than the heliopause. While our simulated TS locations do not reproduce the asymmetry observed by the Voyagers (94 versus 84 au; Stone et al. 2005Stone et al. , 2008, the average of our simulated TS distances is consistent with the average TS distance from the Voyagers (∼89 au). Figure 3 displays cross sections of our simulated heliosphere in the solar equatorial and meridional planes, showing the total plasma speed, radial plasma speed component, and transverse plasma speed component (i.e., transverse to the radial direction). The range of plasma flow speeds in the IHS (between ∼250 km s −1 near the polar directions and close to ∼0 km s −1 near the upwind heliopause region) greatly affects the frame transformation in different directions in the sky. Also, in many parts of the sky, there are both radial and transverse plasma flow components with respect to an observer near the Sun. While radial flows produce the largest changes in ENA energy and flux during the frame transformation, transverse flows also contribute a significant amount to the transformation, especially at low ENA energies (<1 keV). This is explained in more detail below. Examples of the plasma flow velocity components extracted from our simulation, u p , using Equations (15) and (16), are shown in Figure 4 as projections onto a sky map, including the total plasma speed, the radial plasma speed, and the transverse speed. We note first that the LOS-weighted properties are calculated as a function of the ENA energy passband. However, within the energy ranges in the plasma frame (∼1-6 keV), the energy-dependent term s v v ex rel,p rel,p ( ) only changes by a factor of ∼1.5 (e.g., Zirnstein & McComas 2015). Added to the fact that LOS averaging is a weighted integral along flow streamlines from the TS through the IHS and that the flow properties are integrated over wide energy ranges for each energy passband, the LOS-averaged plasma quantities are not significantly different as a function of ENA energy.
The LOS-weighted plasma flow is slower near the nose of the heliosphere where the simulation's stagnation region is approximately located-a product of the steady-state assumption and not necessarily in a dynamic state. Note that our stagnation region appears to be ∼5°-10°north of the stagnation region derived from IBEX data by McComas and Schwadron et al. (2014), but our simulated stagnation region is offset from the nose toward the southern-port direction (left and down in Figure 4) in a direction similar to that deduced by McComas and Schwadron et al. (2014). The fast SW at high latitudes is clearly visible in all flow components. As one looks away from the nose toward the flanks, the transverse flow increases (at longitudes approximately between −30 and 0°and −150 and −180°) due to the deflection of the plasma away from the heliopause, and the radial flow increases (at longitudes approximately between 0 and +30°and 180 and +150°) as the plasma flows down the heliotail. The total plasma flow speed increases (asymmetrically) on the port and starboard sides of the heliosphere, indicating that the IHS cross-sectional area perpendicular to the plasma flow direction is decreasing and the flow speed increases to maintain flow continuity. Closer to the central heliotail direction (∼60°-90°longitude), the plasma flow speed decreases because the flow is mostly radial, and charge exchange with interstellar neutrals continues to remove energy and momentum from the plasma, reducing both the plasma pressure and flow speed.
Transforming IBEX ENA fluxes from the spacecraft to the IHS plasma frame yields changes in the ENA energy (and thus flux) and direction of emission in the plasma frame. This is illustrated in Figure 5, which shows the change in ENA energy, change in ENA emission angle, and change in velocity (transverse component) when transformed to the plasma frame (also calculated using Equations (15) and (16)). The change in ENA energy is clearly related to the IHS plasma flow velocity relative to the observer. At high latitudes, where the IHS plasma flow is faster, the source plasma is moving quickly with respect to the observer, and the ratio of ¢ E Ẽ is much higher compared to lower latitudes. The change in ENA energy is smallest near the nose direction, where the IHS flow slows down and is diverted around the heliopause, and near the central heliotail where the plasma is slowing down due to momentum loss by charge exchange. The change in ENA energy is less pronounced at higher energies (∼4.3 keV) than lower energies (∼0.7 keV), which is important to note as this alters the spectral slope of the proton spectrum (see Section 3.2). Lastly, there is a moderate increase in ENA energy near the port and starboard sides that form a ring-like structure around the flanks of the heliosphere. As we discussed earlier, the IHS plasma speed increases on the flanks of the heliosphere as it propagates toward the heliotail (see Figure 4). Figure 5 also shows the change in emission angle of ENAs when transformed to the plasma frame, i.e., the angle between the ENA velocity in the plasma frame and the inward radial direction (see Figure 2). Similar to the previous panels, the angular direction of the emitted ENA is significantly different in the plasma frame at high latitudes due to the faster IHS plasma speed. The change in ENA angular direction exceeds ∼18°at high latitudes, or ∼3 IBEX pixels. Also, similar to the  (15) and (16)), but visible differences between the maps at different energies are small. Here we define = E E sc˜, the ENA passband central energy in the spacecraft frame. energy transformation, the ENA emission angle is significantly offset from the radial direction near the port and starboard sides of the heliosphere due to the significant transverse-flow component of the IHS plasma. There is no significant change in emission angle near the nose and heliotail directions because the IHS plasma flow is either stagnating or moving radially away from the observer (change in emission angle only depends on the transverse-flow component).
Finally, Figure 5 shows the change in ENA velocity direction emitted in the plasma frame. Specifically, this shows the transverse component of the ENA velocity in the plasma frame divided by the particle speed in the spacecraft frame. While IBEX observes ENAs coming from the radial direction in the spacecraft frame (see Figure 2), the ENAs were actually emitted with large transverse velocity components in the plasma frame due to the transverse flow of the IHS plasma. The vector arrows illustrate the direction of the ENA particle emitted in the plasma frame and the length/color of the arrow illustrates the relative magnitude. The majority of the vector arrows point toward the nose of the heliosphere because the IHS plasma flows away from the point of highest pressure, i.e., near the nose or stagnation region. This does not mean that ENAs observed by IBEX are coming from directions different from what it observes in the spacecraft frame (IBEX-Hi and IBEX-Lo are "single-pixel" cameras with relatively narrow fields of view). Rather, it implies that IBEX observes some ENAs emitted from the nonradial part of the parent plasma distribution.

Proton Spectra in the Heliosheath Plasma Frame
Utilizing the plasma and neutral properties extracted from the MHD simulation, we transform the observed ENA fluxes to LOS-averaged proton fluxes in the plasma frame. Figure 6 shows examples of frame-transformed spectra from the Voyager 1 and southern polar directions at various times throughout the IBEX epoch. ENA fluxes in the spacecraft frame are plotted at the central energies of each of IBEX-Hiʼs energy passbands.
Transforming the ENA flux to the plasma frame shifts the observed energies higher and increases the fluxes (as the IHS plasma is moving away from IBEX). Converting the ENA flux to proton flux (Equation (11)) does not change the energy of the particle, but its flux is increased due to the energydependent rate of charge exchange. The ENA and proton spectra change over time, although we note that the change in time is not reproduced in the heliosphere simulation. In each panel of Figure 6, we also show the probability of the proton spectra being consistent with single power laws, which is described in more detail in Section 3.4. The probability is commonly referred to as the "p value." As described in Appendix C.5, the p value is determined by calculating the chisquare (χ 2 ) goodness-of-fit statistic for a line fit to the data in log-log space for each direction in the sky and calculating the integral of the χ 2 probability distribution function.
Next, we present our spectral indices of the derived proton spectra in the plasma frame. The spectral index is calculated between each energy passband in the plasma frame, ¢ E j , where the subscript j corresponds to each IBEX-Hi energy passband, corrected for the change in energy from the spacecraft to plasma frame (see labels in Figure 6). The spectral index in the plasma frame, g ¢ j , is calculated as ( where we propagate the uncertainties of the ENA fluxes directly through the frame transformation and charge-exchange conversion. As frame transformations are energy dependent, the spectral indices of ENAs observed in the spacecraft frame are different than those of their parent protons in the plasma frame. As demonstrated in Figure 7, the spectral indices for protons in the plasma frame are larger than those for ENAs in the spacecraft frame on average by ∼25% for ∼0.7-1.1 keV and ∼3% for ∼2.7-4.3 keV (see Appendix C.3 for equations of the mean and standard deviation). The general increase in the spectral index is due to the motion of the IHS plasma away from IBEX  (11) and (12) assuming L = 30 au. For results presented later in this study, we calculate the spectral index of the proton flux between each energy passband, i.e., g¢ 1 , g¢ 2 , g¢ 3 , and g¢ 4 , and over the entire energy range as a single-power-law fit. The normalization of the proton flux (i.e., the values for n H and L) does not significantly affect the spectral index results. The probabilities of the proton spectra being consistent with single power laws are shown in the top right as the p value. with speeds comparable to the observed speed at the detector (see Figures 3 and 4). Note, however, that the distribution of spectral indices over the entire sky is quite wide (standard deviation ∼1). Figure 8 shows histograms of the proton spectral indices in the plasma frame between each energy passband accumulated over different regions of the sky: (1) all sky, (2) front hemisphere, (3) back hemisphere, (4) high latitudes, and (5) low latitudes. Proton spectra in the back hemisphere are, on average, slightly steeper than those in the front hemisphere. At high latitudes, the proton spectra are harder at the highest energies, presumably due to the fast SW flows, while the opposite is true at low latitudes where there appears to be a significant rollover above ∼2 keV. As the spectra are already Figure 8. Histograms of proton spectral index in the plasma frame between each energy passband, g¢ i . We show results over the entire sky from 2009 through 2016, as well as separately for the front hemisphere (255°. 7-90°< f p < 255°. 7 + 90°), back hemisphere (75°. 7-90°< f p < 75°. 7 + 90°), high latitudes (|90°-θ p | > 45°), and low latitudes (|90°-θ p | 45°). Vertical lines below the histograms show the locations of the mean indices. "corrected" for the relative plasma flow velocity, these spatially dependent features are likely indications of (1) the TSprocessed PUI distributions from slow versus fast SW plasma parcels and (2) modifications to the distributions as they propagate through the IHS in different directions of the sky.
It must be made clear again that at various times throughout the IBEX epoch that the ENA and proton spectra appear to exhibit spectral breaks (both "knees" and "ankles") in different directions of the sky, as shown in Figures 6 and 8. These spectral breaks, however, are less prevalent than previously thought earlier in the mission, which we discuss later in Section 3.6, and may not be statistically significant, which we discuss more in Sections 3.4.

Importance of Transverse Plasma Flows
Before we continue our analysis of the proton spectra, we briefly demonstrate the importance of accounting for the transverse-flow component of the IHS plasma. The frame transformation from the spacecraft frame to the IHS plasma frame (Equation (2)) shows that the plasma flow component parallel to the observer's LOS (i.e., radial component) introduces a first-order transformation proportional to the ratio of the relative frame speed and the measured ENA speed, u v rel˜. The transverse-flow component (i.e., transverse to the radial direction) only contributes to a second-order transformation proportional to u v rel 2 (˜) . However, depending on the ENA energy, the contribution of the IHS transverse flow to the frame transformation can be significant. Figure 9 illustrates this for ENAs observed by IBEX-Hi at ∼0.7 and ∼4.3 keV. At the lowest energy passband for IBEX-Hi, the change in ENA energy (and flux) due to the transverse-flow component of the IHS plasma is at most ∼30% due to the fast SW at mid to high latitudes on the noseward side of the heliosphere. Near the flanks of the heliosphere where the IHS plasma is diverted around the heliopause and begins to the flow down the heliotail, the transverse component accounts for <20% of the frame transformation. At the highest energy passband for IBEX-Hi, the transverse-flow component accounts for <15% of the frame transformation.
Although we do not analyze IBEX-Lo observations in this study, we can show how the IHS plasma transverse-flow component becomes more important at lower ENA energies (note that we include the effects of IBEX-Loʼs energy response assuming a log-Gaussian response with dE/E = 0.7; McComas et al. 2012). Figure 9 shows results for IBEX-Lo energy passband 4 (central energy ∼0.11 keV). At these low ENA energies, the transverse component of the IHS plasma flow has a significant impact on the change in ENA energy and flux between the two reference frames. The transverse component accounts for up to ∼60% of the change in ENA energy and flux at mid to high latitudes on the noseward side of the heliosphere and up to ∼40% at the flanks, due to the high transverse component of the IHS plasma as it flows toward the heliotail.

Mean Proton Spectral Index-"Single" Power Law?
We continue our analysis of the proton spectra by determining whether or not the proton spectra derived from IBEX-Hi are statistically consistent with single power laws, which may possibly be evidence for kappa distributions of protons in the IHS (e.g., Livadiotis et al. 2011;. A similar type of analysis was performed for IBEX-Lo and IBEX-Hi observations in the Voyager 1 direction by Fuselier et al. (2014) but without transforming the data to the IHS plasma reference frame. We calculate the mean proton spectral index per pixel in the sky by fitting a line to each spectrum in log-log space including statistical uncertainties (see Appendix C.2) and demonstrate the goodness of the fits to the data by the chi-square statistic and p value (see Appendix C.5). Figure 10 shows histograms of the weighted mean proton spectral index per pixel in the sky, only including pixels where the p value is above some threshold value. Typically, a p value <0.05 is used to reject the null hypothesis of the test (in our case, the null hypothesis is that the spectra are consistent with a single power law per pixel). However, cases where the null hypothesis cannot be rejected (p value >0.05) are not necessarily strong evidence that the spectra are single power laws. Therefore, to illustrate higher probabilities of singlepower-law spectra in the sky, we show several cases of histograms where the p-value threshold is increased and how it changes the weighted mean and standard deviation of the accepted population.
First, Figure 10 demonstrates that each population with different p-value thresholds yields nearly identical estimations of the mean and standard deviation of the parent population (mean and standard deviation of the histograms; see Appendix C.4). Second, as is evident from the right panel of Figure 10, the percentage of the data set that is excluded from a specific p-value threshold is approximately consistent with that expected from statistical fluctuations. This is evidence that the majority of the IBEX-Hi data analyzed in this study are consistent with the hypothesis that the proton spectra are described by single power laws, where the mean index over the sky is ∼2.1 and the standard deviation ∼0.4. It is also interesting to note that there are very few occurrences of spectral indices below 1.5 (∼5% of the entire data set), which is a lower limit for kappa distributions (e.g., . Figure 11 expands on the significance of our results by showing the observed chi-square distribution of the fitted results compared to the expected distribution where deviations from a power law are due to statistical fluctuations. The agreement between the two (see also the cumulative distribution function on the right side of Figure 11) further supports the fact that the majority of the IBEX-Hi data set appears to be consistent with the single-power-law hypothesis.
Regardless of the statistical probabilities of individual spectra being single power laws, it is clear that the mean spectral index in the sky is approximately ∼2.1 with a standard deviation of ∼0.4. In other words, the ±1σ slope range of the LOS-weighted proton fluxes in the IHS is á ¢ ñ µ ] . If we interpret the proton spectra as being kappa distributions, then the mean kappa index is most likely close to ∼2.1, and the ±1σ range k á ¢ñ ∼ [1.7, 2.5]. Under this interpretation, the majority of the proton spectra lie in the far-equilibrium regime (<2.5; . Livadiotis et al. (2011) found a kappa index mode of ∼1.7 from IBEX observations during 2009, except that by including the relative motion of the IHS bulk plasma with respect to the observer, we find that the spectral index we obtain is likely larger by, on average, ∼20%. However, the power-law slope over a specific energy range of a kappa distribution only reaches the value of the kappa index at energies much larger than the mean speed of protons in the plasma frame.
Moreover, the energy range of IBEX-Hi observations in the spacecraft frame (∼0.5-6 keV) is lower than the energies transformed into the plasma frame (see, e.g., Figure 6). Therefore, the proton power-law spectral indices derived from IBEX observations are not necessarily equivalent to the kappa index of a proton kappa distribution. To demonstrate this, we have calculated the spectral slopes of a kappa distribution with different kappa indices, temperatures, and plasma flow velocities with respect to the spacecraft. For large plasma velocity away from the spacecraft, small plasma temperature, and small kappa index, the observed spectral index will be closer to the actual kappa index. Typical values for the IHS plasma temperature based on simulations is ∼1 to 4 × 10 6 K, and the plasma flow velocity away from the IBEX spacecraft is between ∼50 and 250 km s −1 . In order to reproduce the mean proton spectral index from the results of our analysis (2.1) Figure 10. (Left) Histograms of the mean proton spectral index for data from 2009 through 2016. The significance of the spectral index as being consistent with a single power law is determined by only including pixels whose p value is greater than 0.05, 0.2, and 0.4 (maximum p value is 0.5). We also show the weighted mean index and weighted sample standard deviation of each histogram (see calculation in Appendix C.4). Percentages of the available data set that are included in each histogram are also shown. (Right) Population percentages for each p-value threshold between 0 and 0.5. The percentage of the data with p values above a certain threshold (dashed curve), added to twice the p-value threshold (solid curve), is shown as the total (dotted curve). This ideally equals 100% if deviations of measurements from a single-power-law fit are purely due to random, statistical fluctuations.
within the range of plasma temperatures and speeds described above, the actual kappa index of the proton distribution lies between ∼2.15 and 2.7. This implies that if the proton distribution in the IHS can be represented by a single kappa distribution in each direction of the sky, the mean kappa index is likely ∼3% to 30% larger.
The spatial dependence of the statistical significance of the single-power-law hypothesis is demonstrated in Figure 12, which presents sky maps of the spatial distribution of the proton spectral index presented in Figure 10, shown for different thresholds of the p value. For a given 6°× 6°pixel in the map, we weight-averaged any spectral indices over the IBEX data set (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)) that satisfied the p-value condition, as in Figure 10. There is no significant spatial or latitudinal dependence on spectra that can or cannot be represented by a single power law, in contradiction to the results of Livadiotis et al. (2011. The most likely reason for this, which we explain in more detail in Section 3.6, is a source of nonheliospheric background. Finally, similar to the ENA spectral indices observed in the spacecraft frame (e.g., McComas et al. 2020), the port and starboard heliotail lobes exhibit the steepest spectra, which appear as depletions in ENA flux at higher energies from the flanks compared to the nose and tail directions (Figure 12; see also Figure 5 in Zirnstein et al. 2016a). The origin of the port and starboard ENA lobes was explained in detail by Zirnstein et al. (2016a), owing to a combination of (1) the relative PUI to total proton density at the TS, (2) total proton density, (3) PUI cooling length as a function of angle from the nose to the tail, and (4) the plasma flow structure (see Figures 7-9 in Zirnstein et al. 2016a). The north and south heliotail lobes (i.e., where the heliotail ENA lobes "split" at energies >2 keV) exhibit harder spectra due to the relatively high mean energy of the fast SW plasma propagating down the heliotail ).

Time Dependence of the Proton Spectra
For the majority of this analysis, we applied our frame transformation and ENA-to-proton conversion to all IBEX-Hi observations from 2009-2016, to provide the best statistics, where we assumed the average SW structure was like solar minimum. However, it is clear from IBEX ENA observations (e.g., Zirnstein et al. 2017aZirnstein et al. , 2020aMcComas et al. 2020) that the evolution of the solar cycle and disappearance of the fast SW at high latitudes during the recent solar maximum had a significant impact on the ENA spectra. This suggests that, at high latitudes, our analysis overestimates to some degree the ENA/proton spectral index during ∼2014-2016, a few years after the solar cycle went through the recent solar maximum in ∼2012-2016 (see, for example, Figure 8 in Zirnstein et al. 2020a)-however, we note again that this is not strictly true because the northern and southern hemispheres of the heliosphere evolved differently owing to the asymmetric evolution of the fast SW.
We explore this time dependence further in Figure 13, which shows all-sky histograms of the mean proton spectral indices in the plasma frame for 2009-2013 versus 2014-2016. It is clear that earlier in the mission, when the ENA fluxes were created during largely solar-minimum-like conditions (more accurately represented by our heliosphere simulation), the spectral indices are smaller by ∼10% compared to the results shown in Figure 10. In the latter period (2014-2016), the spectral indices are ∼15% larger. A fully time-dependent simulation of the heliosphere and frame transformation analysis would likely be necessary to better understand the time-dependent proton spectral indices, which will be the topic of a future study.

Broken Spectra and the "Ion Gun Effect"
As alluded to in Section 3.3, our analysis reveals that there is no significant latitudinal dependence on whether or not proton spectra have a high probability of being single power laws, i.e., like a kappa distribution. Earlier analyses of proton spectra derived from IBEX-Hi observations suggested that the highlatitude spectra early in the mission (∼2009-2011) cannot be represented by a single power law due to significant breaks in the spectra . The high-latitude broken spectra were analyzed in detail by Dayeh et al. (2012), who explained with reasonable evidence that the broken spectra may be a result of a superposition of ENA populations produced from both fast and slow SW parcels moving through the IHS (see, for example, Figure 5 in Dayeh et al. 2012). A significant spectral break with sufficiently small statistical uncertainty would produce a small p value and a low probability of originating from a single-power-law distribution. However, during our analysis, we tested our methods by comparing our results to those presented in earlier studies, e.g., Dayeh et al. (2012). In doing so, we noticed significant differences in the ENA spectral indices. Figure 14 shows histograms of the proton spectral indices between each energy passband over the entire sky (not the mean spectral index over all passbands), only including observations from 2009 to 2011. This is similar to the data sets analyzed by Livadiotis et al. (2011  Interestingly, there is a clear difference in spectral indices between energy passbands 2 and 3 (gray histograms), a moderate difference between energy passbands 4 and 5 (green histograms), and less difference between other energies. While not shown here, these differences are even more pronounced at high latitudes, like that shown by Dayeh et al. (2012). By testing all other IBEX all-sky map data releases (5 yr-Data Release #7, and 7 yr-Data Release #10), we found that the spectral indices were similar for releases #7, #10, and #16, but vastly different in release #4. The differences were not related to CG or survival-probability corrections. In fact, the major change in data processing that occurred between Data Releases #4 and #7 was the identification and removal of the "ion gun" background. As explained by McComas et al. (2014) in Data Release #7, a background associated with ambient neutral atoms (unrelated to heliospheric ENAs; e.g., molecules outgassed from the spacecraft) that were ionized largely by SW electrons within the IBEX-Hi positive collimator region was identified as a significant background source. This background mostly affected IBEX observations during 2009-2012 when the negative collimator voltage was decreased in order to prevent SW ions from being pulled into the instrument. In 2013, the negative collimator voltage was increased enough to prevent any significant number of SW ions and electrons from entering the instrument, the ion gun background was also significantly reduced for observations post-2012, and the ion gun background was removed in all observations for all data releases starting with Data Release #7 (McComas et al. 2014).
The removal of the ion gun background, which was most prevalent during operations in 2009-2012 and data analyses prior to Data Release #7, effectively changed the spectral index distribution and removed evidence suggesting that there is a significant break in the ENA (and thus proton) spectra, as we show in Figure 14. This background removal also implies that the explanation provided by Dayeh et al. (2012), which still may be true in some situations, is not as significant an effect as previously thought. It also implies that the conclusion by Livadiotis et al. (2011) that the majority of the high-latitude ENA spectra cannot be represented by a kappa distribution, due to the spectral breaks, is also likely not true. We note, however, that there are still some spectral breaks in the ENA measurements, but there does not appear to be any latitudinal correlation (Figure 13), and the statistical uncertainty of the spectral breaks is large. A time-dependent analysis of the proton spectra in the plasma frame may reveal any potential latitudinal correlation with single-power-law proton spectra; however, seeing that there are no significant spectral breaks in 2009-2011 during solar minimum (Figure 14), it seems unlikely to be a major effect.

Summary and Conclusion
IBEX is capable of probing the IHS plasma via the measurement of ENA emissions. However, due to the relative motion of the source plasma in the IHS with respect to the spacecraft, the parent proton distributions (which are typically assumed to be isotropic in the plasma frame at these energies) inferred from ENA observations will not be accurate if the frame transformation (i.e., CG effect) between the reference frames is not properly taken into account. To demonstrate this, we simulate the SW-VLISM interaction with a 3D MHD/ kinetic simulation, extract the LOS-averaged IHS plasma and neutral H flow properties from the simulation, transform the ENA fluxes observed by IBEX to proton fluxes in the IHS plasma frame, and analyze the parent proton spectral properties.
The frame transformation affects both the energy of the ENA and its direction of propagation with respect to the moving plasma frame. These effects are minimal near the high-pressure stagnation region of the upwind heliosphere where the plasma slows down considerably before encountering and diverting around the heliopause, and also near the center of the heliotail where the IHS plasma slows down due to continuous loss of energy and momentum by charge exchange with neutrals from the IHS and interstellar medium. However, at mid to high latitudes and near the flanks of the heliosphere, the IHS plasma has both high radial and transverse-flow speeds and thus the frame transformation greatly affects the energy and angle in the sky in which the ENA is emitted in the plasma frame. Moreover, the transformation is greatest at high latitudes when the Sun emits fast SW speeds (∼750 km s −1 ). Using our heliosphere simulation, we estimate that the parent proton energies are higher by a factor >2 compared to ENAs observed by IBEX-Hi energy passband 2 (∼0.7 keV) at high latitudes and <1.5 for IBEX-Hi energy passband 6 (∼4.3 keV). Because the IHS plasma generally flows away from the high-pressure region of the upwind heliosheath, the majority of ENAs actually originate from IHS protons whose velocities before charge exchange had noseward components in the plasma frame.
Typically, only the radial (or LOS-parallel) component of the source plasma motion is considered when calculating the CG effect. However, the IHS plasma has a significant transverseflow component in many directions of the sky. While the transverse component is mainly responsible for the change in ENA emission angle between the spacecraft and plasma reference frames, it also changes the energy and flux of the ENAs as a second-order effect. We demonstrated that the  Comparison of proton spectral indices over the whole sky from 2009-2011, extracted from two different IBEX data releases. The first (Release #4) did not account for the removal of background due to the ion gun effect. All releases following #4, including the latest data release (Release #16), remove this background. transverse component of the IHS plasma flow accounts for up to ∼30% of the change in ENA energy and flux for IBEX-Hi energy passband 2 (∼0.7 keV in the spacecraft frame) measurements when transforming to the plasma frame, and up to ∼15% for passband 6 (∼4.3 keV in the spacecraft frame). However, IBEX-Lo measurements are much more sensitive. For IBEX-Lo energy passband 4 (central energy ∼0.1 keV), the transverse-flow component contributes up to ∼60% in certain directions of the sky. Therefore, we suggest that analyses that aim to transform the lowest energy IBEX-Hi ENA fluxes and any IBEX-Lo ENA fluxes to proton fluxes in the IHS plasma reference frame (e.g., Fuselier et al. 2014;Galli et al. 2017) should account for the transverse-flow component, which can be estimated from a 3D simulation of the heliosphere, but also to use the term s v v v ex rel rel ( ) instead of s v ex ( ) during the conversion. Moreover, it is worth noting that the frame transformation implies that ENAs observed by IBEX-Lo at ∼0.1 keV actually originate from protons with energies ∼0.5 keV in the plasma frame.
By calculating the spectral index of the ENA fluxes between each IBEX-Hi energy passband, we compared the spectra of ENAs in the spacecraft frame to those of protons in the plasma frame. We found that, on average, the inferred proton spectral index is ∼20% higher than the observed ENA spectrum, with small variations depending on the direction in the sky and the energy range. Moreover, by calculating the mean proton spectral indices by performing a linear fit over the IBEX-Hi energy range in log-log space, we find that the majority of observations are statistically consistent with the hypothesis of being single power laws, perhaps originating from a kappa distribution of protons in the IHS, where the mean proton spectral index in the sky over 2009-2016 is ∼2.1 with a standard deviation of ∼0.4. This implies that the proton distribution in the energy range IBEX-Hi observes lies in the range á ¢ ñ µ ¢ ¢ -- ] . There does not appear to be any significant spatial dependence on this result, although this may depend more on time variability, which we currently do not account for in our methodology.
We find that earlier analyses of IBEX-Hi observations that identified significant breaks in a large number of ENA spectra in ∼2009-2011, especially at high latitudes, were actually a product of the ion gun background (McComas et al. 2014). Analyses that utilized observations in Data Release #4 or earlier did not account for this background, which mostly affected fluxes detected in energy passbands 2 (∼0.7 keV) and 5 (∼2.7 keV). Comparing our derived proton spectral indices in 2009-2011 using Data Release #4 and comparing observations in the recent Data Release #16 (McComas et al. 2020) reveal that spectra where fluxes between passbands 2 and 3 were significantly steeper than other energies were likely a product of this ion gun background. The background was analyzed and has been removed for every data release since Data Release #7 (McComas et al. 2014) for all years. This analysis highlights some of the complexities in making heliospheric ENA observations, the care that the IBEX team has demonstrated in removing every known background in the ENA measurements, and the importance for the heliospheric community to utilize the most up-to-date IBEX data releases when analyzing or modeling ENA observations. Finally, based on the results presented in this study, it is clear that any future analyses of the IHS plasma properties utilizing heliospheric ENA measurements should take into account the CG effect due to the relative motion of the IHS plasma with respect to the observing spacecraft, including both the radial and transverse-flow components. This is particularly important for observations of ENAs made at energies below ∼10 keV in the spacecraft frame, applicable to IBEX and Interstellar Mapping and Acceleration Probe (IMAP) measurements with the IMAP-Lo and IMAP-Hi instruments ).
This work was funded by the IBEX mission as a part of the NASA Explorer Program (80NSSC20K0719) and the IMAP mission as a part of NASA's Solar Terrestrial Probes (STP) mission line (80GSFC19C0027). E.J.Z. also acknowledges support from NASA grants 80NSSC18K1212 and 80NSSC20K0783. M.A.D. acknowledges support from NASA grant 80NSSC20K0783 and thanks George Livadiotis for helpful discussions related to this manuscript. Ulysses SW data and OMNI SW data were extracted from the SPDF (https:// spdf.gsfc.nasa.gov/). The work reported in this paper was partly performed at the TIGRESS high-performance computer center at Princeton University, which is jointly supported by the Princeton Institute for Computational Science and Engineering and the Princeton University Office of Information Technology's Research Computing department.
where r is the position vector, and we use κ = 2 (for our main results). The distribution ¢ f p,t represents protons that are transmitted across the TS, including both core SW ions and PUIs that may experience preferential acceleration at the TS and form a suprathermal tail. This enables us to reproduce the mean power-law-like spectra observed by IBEX. The evolution of the distribution function as a function of distance l along a streamline through the IHS is given as where the density ratio introduces compression effects along the streamline, the exponential expression describes the loss of energetic protons by charge exchange, σ ex is evaluated at v rel , and v rel is the relative speed of interaction of a proton at speed ¢ v in the plasma frame with a Maxwellian neutral H distribution r The power-law slope of the charge-exchange cross section, ν, is evaluated as , Δu/v th,κ = 1 or Δu/v th,κ ? 1, ν + 2 > 0, and κ − ν − 1 > 0, which are satisfied for our purposes.
The pressure source term H P is given as the sum of two different sources, H P = H P,1 + H P,2 , where H P,1 is the gain in injected PUI energy by charge exchange between the proton kappa distribution and the neutral H distribution (DeStefano & Heerikhuisen 2020, Equation (41)

C.2. Mean Index per Pixel
When we calculate the mean spectral index over the entire IBEX-Hi energy range per pixel, we assume that the measurements represent a single power law and any deviations from their mean is due to measurement uncertainties. This assumption is tested by using the p value. To calculate the mean index per pixel, we fit a line of the form A + Bx to the proton flux in log-log space via linear regression and minimize the chi-square statistic (see Equation (55)). The spectral index is determined by the best-fit value for −B (as B is the spectral slope, −B is the spectral index). We use this method for results presented in Figure 6 (to calculate the p value), the individual data points shown in Figures 10 and 11, and the individual data points before time-averaging in Figure 12.
We note that the method of fitting a line to log-transformed data, under the assumption that the data are well represented by a power law, is a common practice to derive the power-law slope. However, it is not always the case that this method is appropriate. Depending on the distribution of the observed sample and their uncertainties, it may be more appropriate to fit a power law directly to the nontransformed data via nonlinear, least-squares regression. To test the sensitivity of our results to these different methods, we repeated our analysis of Figures 10 and 11 by fitting a power law of the form Ax − B to the nontransformed data using the Levenberg-Marquardt method (Markwardt 2009). We find that while the best-fit spectral indices appear to be slightly smaller when using the power-law fitting method, the results are not significantly different. For the power-law fitting method, the mean spectral index for all three histograms in Figure 10 (left panel) is 2.10, and the standard deviation is between 0.43 and 0.44. This is a change of ∼2%-3% for the mean spectral index and ∼5%-8% for the standard deviation. The chi-square distribution in Figure 11 is similar in both methods, where the majority of the observations (90%) still appear to be consistent with a single power law.

C.3. Mean of the Spectral Index Histogram between Energy Passbands
For results presented in Figures 7, 8, 13, and 14 that show histograms of spectral indices between each IBEX-Hi energy passband, we make the simplifying assumption that the measurement uncertainties are small compared to the parent standard deviation and calculate the unweighted mean and sample standard deviations of the individual histograms (shown as colored text in the figure panels) using Equation (50) and the unweighted version of Equation (49) (when s s ¢ ¢ j  ).

C.4. Mean of the Spectral Index Histogram for a Single Power law
For results presented in Figure 10, which shows histograms of the mean spectral indices culled with different p-value criteria, we utilize the derivation of Equation (48) to calculate the weighted mean of the histograms and the corresponding standard deviation. As the parent standard deviation is not known, we use an iterative procedure to determine both the mean and standard deviation. First, we estimate the mean of the data in each histogram in Figure 10 using Equation (50). Then, the parent standard deviation is estimated as where we subtract the measurement uncertainties from the total deviations of the measurements about their mean. The derivation of Equations (47) and (48) can also be expressed as the minimization of the χ 2 estimator, given as where m ¢ x represents a range of possible values to find the minimum of χ 2 . In the next step, we use the estimate of the parent standard deviation (s ¢ est ) from Equation (52) to obtain a new estimate of the mean m ¢ est from the minimum χ 2 in Equation (53) over a range of values for m ¢ x . After finding the value for m ¢ x that minimizes χ 2 , we substitute m m ¢ = ¢ x est in Equation (52) and repeat the process until consecutive values for m ¢ est and s ¢ est converge to within an accuracy of 0.005. We note that only a few iterations are required before convergence.
The results of the final iteration are shown in Figure 15 as the minimum of the reduced χ 2 (ν = degrees of freedom). Reduced χ 2 curves for all three histograms in Figure 10 are shown in Figure 15, but the results are nearly identical, and they overlap. Uncertainties of the mean are calculated using the propagation of error of Equation (48) ,min , differencing the results from s ¢ min , and taking the maximum difference. Note that the uncertainties of the mean and standard deviation are much smaller than the parent standard deviation.

C.5. Calculating the p Value of a Single-power-law Fit
To determine the goodness of fit of a single power law to each spectrum (per individual pixel at a certain time), or in other words the probability that the spectrum can be represented by a single power law, we utilize the χ 2 statistic The probability given by Equation (57) determines the goodness of fit of a single power law (i.e., the mean) to each spectrum, or in other words, spectra with a high p value cannot reject the null hypothesis (that the spectral indices have zero variance about their mean). Because the fit is either bad when the p value → 0 (χ 2 → ∞ ) or overfit when the p value → 1 (χ 2 → 0), we decide to find the minimum of the integrals on either extreme of χ 2 ([0, χ 2 ] versus c ¥ , 2 [ )) such that the highest probability occurs when the p value ≅ 0.5 or χ 2 ≅ ν.

Appendix D Effect of Different Heliosphere Simulations on Spectral Index Results
In this section, we test the dependence of our spectral index analysis results presented in Figures 10 and 11 on the heliosphere simulation used to calculate the CG frame transformation and charge-exchange conversion of the IBEX data. Specifically, we report the mean and standard deviation of proton spectral index for all data whose p value >0.2, similar to the results shown in Figure 10, for different heliosphere simulations. Scenarios #1-5 utilize simulations performed in this study with slight variations. Scenario #6 utilizes a simulation from Shrestha et al. (2020). Unless otherwise specified, all other modeling assumptions and boundary conditions are the same as those described for the main results of this study.
(1) This study: (mean, standard deviation) = (2.15, 0.41) (2) Same as (1), except proton distribution κ = 2.5 in the IHS: (2.13, 0.40) (3) Same as (1), except proton distribution κ = 3.0 in the IHS: (2.16, 0.41) (4) Same as (1) We find that there is no significant dependence of our results on the proton distribution kappa index (#2-3 show ∼1% change). There is also no significant dependence on using solar-maximum-like SW conditions for IBEXʼs observations in 2014-2016 (#4 shows ∼3% change), when ENA spectra are steepening as a response to the recent solar maximum in solar cycle 24 (Zirnstein et al. 2020a). There is a minor change in the results if the SW is assumed to be uniformly slow at all latitudes for the entire IBEX data set in 2009-2016 (#5-6 show ∼7% change). However, scenarios #5 and #6 are not realistic, as we know there is fast SW at mid to high latitudes corresponding to IBEX observations in 2009-2013 (e.g., Sokół et al. 2015Sokół et al. , 2020Zirnstein et al. 2017aZirnstein et al. , 2020a. We only show scenarios #5 and #6 as examples of the possible change in our results in unrealistically different SW conditions. We note that in all of these scenarios, the majority of the proton spectral indices derived from the IBEX data set are still statistically consistent with single power laws, with results similar to those shown in Figures 10 and 11. Our results appear robust based on the fact that different kappa indices of the proton distribution function do not significantly change our results, and the introduction of slower SW at higher latitudes in 2014-2016 as a simple approach to mimic the behavior of the SW in solar cycle 24 only slightly decreases the mean proton spectral index in the sky. In reality, however, the heliosphere plasma properties are likely more complicated than we can account for, a prime example of which is the large SW dynamic pressure increase observed at 1 au in late 2014 that then propagated to the outer heliosphere and induced a time-, energy-, and directional-dependent evolution in the heliospheric ENA signal after 2016 Zirnstein et al. 2018a). Because of this complex interaction, we did not include 2017-2019 IBEX data in our analysis.
While scenario #4 is possibly more realistic than #1, reality is more complicated based on the fact that the ENA spectra from the downwind side of the heliosphere do not appear to significantly change over time as it does in the forward hemisphere, even in 2014-2016, due to delayed reactions to changes in the SW within the extended heliotail . Moreover, observations of the evolution of the Sun's polar coronal holes (e.g., Karna et al. 2014; Figure 15. Calculation of the mean proton spectral index in the sky for Figure 10 by minimizing the reduced χ 2 statistic. The reduced χ 2 for indices with p value >0.05 (gray), >0.2 (blue), and >0.4 (green) are nearly identical and their curves overlap. We only show the text of the analysis results for spectra with p value >0.2 as an example. The propagated uncertainties of the mean and standard deviation are also shown. See the text in Appendix C.4 for more details.