Relating Energetic Ion Spectra to Energetic Neutral Atoms

Heliospheric energetic neutral atoms (ENAs) originate from energetic ions that are neutralized by charge exchange with neutral atoms in the heliosheath and very local interstellar medium (VLISM). Since neutral atoms are unaffected by electromagnetic fields, they propagate ballistically with the same speeds as parent particles. Consequently, measurements of ENA distributions allow one to remotely image the energetic ion distributions in the heliosheath and VLISM. The origin of the energetic ions that spawn ENAs is still debated, particularly at energies higher than ∼keV. In this work, we summarize five possible sources of energetic ions in the heliosheath that cover the ENA energy from a few keV to hundreds of keV. Three sources of the energetic ions are related to pickup ions (PUIs): those PUIs transmitted across the heliospheric termination shock (HTS), those reflected once or multiple times at the HTS, i.e., reflected PUIs, and those PUIs multiply reflected and further accelerated by the HTS. Two other kinds of ions that can be considered are ions transmitted from the suprathermal tail of the PUI distribution and other particles accelerated at the HTS. By way of illustration, we use these energetic particle distributions, taking account of their evolution in the heliosheath, to calculate the ENA intensities and to analyze the characteristics of ENA spectra observed at 1 au.


Introduction
Heliospheric energetic neutral atoms (ENAs) originate from energetic ions that are neutralized by charge-exchange processes with inflowing interstellar neutral atoms (generally hydrogen but also helium and other heavy atoms) in the heliosheath and very local interstellar medium (VLISM). Neutral atoms are unaffected by electromagnetic fields, and hence propagate ballistically with speeds the same as their parent particles. Measurements of ENA distributions allow one to remotely image the energetic ion distributions in the heliosheath and VLISM, and infer the global heliospheric structure .
There are two well-known relevant populations of ions in the heliosphere, solar wind particles and pickup ions (PUIs). The solar wind particles are the dominant population and follow a Maxwellian distribution. PUIs are ionized LISM neutrals mainly resulting from the charge exchange with solar wind ions and (or) photoionization of neutral atoms by solar UV radiation. PUIs contribute to the power-law tails observed in the solar wind plasma distribution. A simple way to add a power-law tail in the distribution function, thereby modeling the solar wind particles and PUI populations together as a single distribution, is to assume a kappa distribution (Heerikhuisen et al. 2008). The low energy part of the kappa distribution follows a Maxwellian distribution, and at higher energies the kappa distribution has a power-law form. The limitations of the kappa distribution include that the low-and higher-energy relative abundance is fixed by the parameter k, and the observed flat-topped PUI population is absent in the kappa distribution (Heerikhuisen et al. 2008). Moreover, taking a (multiple) kappa distribution(s) to model the particle distribution is a convenient but phenomenological method, since the formation of kappa distributions is not clear.
More efforts are needed to clarify the origin of energetic ion distributions in the heliosheath. Zank et al. (1996) proposed a mechanism for the acceleration of PUIs by repeated reflection from the electrostatic cross-shock potential of a quasiperpendicular shock. The acceleration mechanism, multiply reflected ion acceleration, is efficient and leads to reflected PUIs that could obtain considerable energy (when the thickness of the shock is small). The model predicted that singly and multiply reflected particles can provide the primary dissipation mechanism at the quasi-perpendicular heliospheric termination shock (HTS) whereas thermal solar wind protons experience comparatively little heating across the HTS. The predicted dissipation and preferential heating of PUIs were confirmed by the Voyager 2 observations at the HTS (Richardson 2008;Richardson et al. 2008). Based on the Zank et al. (1996) model, Zank et al. (2010) developed a theoretical model that incorporates the microphysical processes occurring at the quasi-perpendicular HTS mediated by PUIs to construct the shape and relative contributions of three distinct downstream ion populations (hereafter, the Zank10 model): transmitted solar wind ions, transmitted PUIs, and PUIs that are first reflected and then transmitted into the downstream region of the HTS (reflected PUIs). Desai et al. (2012Desai et al. ( , 2014 compared the ∼0.5-6 keV ENA spectra measured by IBEX-Hi along the lines of sight of Voyager 1 and 2 to the Zank10 model. They found that the observed ENA spectra cannot be produced by the core solar wind ions only, but agree to within ∼50% of the Zank10 model, and concluded that transmitted PUIs in the ∼0.5-5 keV energy range and reflected PUI above ∼5 keV energy contribute to the observed ENA flux. At energies higher than the PUI cutoff, suprathermal power-law tails on solar wind distributions are observed both in the quiet and disturbed solar wind near the Earth (Gloeckler 2003). The tails are found to persist to at least 40 au from the Sun (Kollmann et al. 2019), and might continue to be present in the outer region of the heliosphere and be transmitted across the HTS. These suprathermal particles and the reflected PUIs are likely further accelerated at the HTS via the diffusive shock acceleration mechanism. By way of illustration, we use these possible energetic particle distributions to calculate the ENA intensity and discuss the characteristics of ENA spectra.

High-energy Particle Distributions in the Heliosheath
Since the kappa distribution is widely used to model the particle distribution when calculating the ENA flux, a short discussion of the propertiesof the kappa distribution is presented. We then revisit the Zank10 model by invoking diffusive shock acceleration for the reflected particles. The presence and acceleration of the suprathermal PUI tail are also discussed.

The Kappa Distribution
The kappa distribution, proposed to fit both the thermal as well as the suprathermal part of the observed velocity distribution of magnetospheric electrons (Olbert 1968;Vasyliunas 1968), is widely applied in space physics including to the proton distribution in the heliosheath in order to fit the ENA observations (Heerikhuisen et al. 2008;Prested et al. 2008). The isotropic, three-dimensional form of the kappa distribution is written as where Γ(x) is the Gamma function, and θ is the most probable particle speed at which v 2 f (v) becomes maximum. By considering the second moment of the distribution function, the mean energy per particle E m is related to the most probable speed θ through the associated characteristic energy The kinetic temperature T k is defined through E m = 3/2k B T k (Maksimovic et al. 1997;Oka et al. 2013).
Approximately near the core, where  q v k, the distribution function is (Livadiotis & McComas 2013) We can further write it in a Maxwellian form as where q q = + c k k 1 . This is a standard Maxwellian distribution with thermal speed q + k k 1 ( ) and modified density ) . Since Γ(k + α) → Γ(k)k α for k → ∞ , the factor in the front of the density n converges asymptotically to 1. This is another perspective that illustrates the equivalence between a kappa distribution in the limit of k → ∞ and a Maxwellian distribution. Obviously, for  q v k, the distribution function has a power low form ∝v −2 k−2 . The kparameter is a free parameter that serves as a measure of the departure from its Maxwellian characteristics and determines the high-energy power-law tail of the distribution (Lazar et al. 2016). As k decreases, the number of suprathermal particles increases compared to the Maxwellian distribution. Figure 1 illustrates different examples of kappa distributions. All are normalized at v = 0, such that f (0) = 1.
The kappa distribution has a critical limitation for the velocity moments of l order which diverges for k < (l + 1)/2 (Scherer et al. 2017;Lazar & Fichtner 2021). By requiring the existence of the kinetic temperature (l = 2; see, e.g., Equation (2)), the kappa distribution must fulfill the condition k > 3/2. Scherer et al. (2017) introduced a generalization of the standard kappa distribution by adding an exponential cutoff to the power-law tail which yields nondivergent velocity moments, where k ä (0, ∞ ). The combination of multiple regularized kappa distributions can fit the measured ion spectra up to 344 MeV (Scherer et al. 2022). The regularized kappa distribution is given by where U is Kummer-U or Tricomi function and ξ is the cutoff parameter. However, applying (regularized) kappa distributions to model the particle spectrum without the thermal component is questionable, and the derived density and temperature are problematic.
k, the kappa distribution function approaches a Maxwellian distribution, and for  q v k, it has a nonthermal power-law tail. For large values of k, the distribution approaches a Maxwellian distribution. It is obvious that the most probable speed for a kappa distribution is θ.

PUIs and Suprathermal Tails
PUIs are created in the inner heliosphere from the inflowing interstellar neutral gas, mostly via charge-exchange interactions with solar wind ions and/or photoionization by solar UV radiation. The newborn ions from the neutral atoms immediately gyrate around the interplanetary magnetic field with an initial velocity inherited from the bulk velocity of the neutral gas flow. The motional electric field of the solar wind accelerates the newborn ions to the speed of the solar wind. The scattering in pitch angle caused by heliospheric magnetic field fluctuations and waves excited by an instability associated with the ring distribution results in the PUIs being isotropically over the outermost shell in velocity space (Isenberg 1986;Lee & Ip 1987;Williams & Zank 1994;Isenberg & Lee 1996; see Zank 1999Zank , 2015 for extensive reviews). Adiabatic cooling causes the shell to shrink toward the center, while the outermost shell is replaced by newly generated PUIs. A filled shell in the PUI velocity space is thus formed in the outer heliosphere. Since adiabatic cooling is less effective in the outer heliosphere, the cutoff speed of PUIs is not necessarily twice the solar wind speed at the HTS (320 km s −1 ; Kumar et al. 2018).
The filled-shell distribution for PUI protons at the outer heliosphere in the solar wind frame is Vasyliunas & Siscoe (1976, Equation 11(b) where n p is the PUI proton density, c = v − u is proton speed in the solar wind frame, u is the solar wind speed, v is the particle speed in the inertial frame, and Θ is the Heaviside step function. Gloeckler & Geiss (2001) found that the PUI proton spectra measured by Ulysses have a sharp cutoff at v/u ≅ 2 which is the expected cutoff speed of freshly born PUIs in the spacecraft reference frame.
where c is the PUI speed at distance r, and c 0 is the particle speed at the heliocentric distance r 0 . Thus the generalized PUI distribution function is given by For adiabatic cooling, α c = 3/2, and this equation reverts to the standard Vasyliunas & Siscoe model. This extended equation considers nonadiabatic cooling with stronger cooling when α c < 3/2 and additional heating of the particle distribution when α c > 3/2 (McComas et al. 2021. It is worth noting that if heating could significantly change the adiabatic index (particle spectrum), it should also increase the maximum energy of particles. The Ulysses observations of H + PUIs are consistent with adiabatic cooling with a cooling index 1.5 (Gloeckler & Geiss 2001). However, based on SWAP observations between ∼22 and 47 au, McComas et al. (2021) found acooling index with mean 2.1 and standard error 0.45. These measurements strongly suggest that PUIs are experiencing additional heating as they are advected in the outer heliosphere toward the HTS. By modeling the cooling index and the radial distance from the Sun as a power-law relationship, the extrapolated cooling index of the PUI distribution in the HTS is 2.9 ± 0.2. However, the experimental determination of the cooling index depends on the knowledge of ionization rates and their spatial variation, the spatial distribution of interstellar neutral atoms, and the survival probability of PUI (Chen et al. 2014;Swaczyna et al. 2020;McComas et al. 2021). Moreover, the empirical power-law relationship is unlikely to work in the outer heliosphere, e.g., the heating may result from the dissipation of turbulence energy, but turbulence is strongly spatial dependent and may behave differently in the outer heliosphere ). It has been found that a transmitted PUI distribution with a −3/2 spectrum will produce too large a high-energy ENA flux and cannot match the observations (Desai et al. 2014;Zirnstein et al. 2018). A larger cooling index will significantly enhance the deviation (see discussion in Section 3.2). PUIs may be accelerated to higher energies to form a suprathermal population. Within 10 au, the suprathermal component has been observed to be f ∝ v − ν with ν ≈ 5 above the PUI cutoff speed  Hill & Hamilton (2010) analyzed Cassini CHEMS observations of H + , He + , and He 2+ suprathermal tails and showed the presence of v −5 power-law tails in the solar wind frame both during quiet periods as well as during active periods in the solar wind between ∼1 and 9 au from the Sun. Kollmann et al. (2019) used the PEPSSI instrument on the New Horizons spacecraft and measurements from the CHEMS observations on Cassini to present a statistical study of suprathermal ions in the keV to hundred keV energy range. The trajectory of New Horizons is similar to the Voyager 2 trajectory ). Kollmann et al. (2019) found that suprathermal He + ions with energies above the PUI cutoff up to ∼100 keV at 5-40 au have a power-law distribution f ∝ v − ν with 4  ν  6 and the mean of ν is not exactly equal to 5. The measurements do not support a significant radial dependence of the average suprathermal exponent. The distribution function in the solar wind frame can be written as where c mt is the maximum particle speed. Gloeckler & Fisk (2010), Baliukin et al. (2020) incorporated these suprathermal ions in the heliosheath to model the ENA flux. Note that in , they included the suprathermal tail in the simulation, but the distribution function between the suprathermal tail and the standard PUI distribution is continuous at the PUI cutoff speed. With that model, they ruled out the presence of the suprathermal tail of PUI at the HTS since the position of the PUI overshoot will be significantly shifted away from the HTS and thus result in an inconsistency between the simulated and measured magnetic field profile along the shock normal (Kumar et al. 2018). However, the continuity of distribution function between the PUI distribution and the power-law suprathermal tail is not observed by spacecraft (see e.g., Hill et al. 2009;Zhang et al. 2019). To match observations for a suprathermal tail intensity that is much smaller than that of PUIs at the cutoff energy ( f tail (c = u) = f p (c = u)), we require n tail /n p = α c /(ν − 3).

Revisiting the Zank10 Model
The charge separation between ions and electrons across quasi-perpendicular shock results in the formation of a crossshock potential. As discussed by Zank et al. (1996Zank et al. ( , 2010, Lee et al. (1996), ions with mass m i and normal velocity v x in the shock frame such that , 9 x i spec ( ) are reflected by the cross-shock potential f, where x is the direction parallel to the shock normal.
The number density of the reflected PUIs in terms of threshold speed V spec and PUI density n p for the Vasyliunas & Siscoe model is given by Shrestha et al. (2021). Following that method of derivation (Shrestha et al. 2021), the number density of reflected PUIs for a generalized PUI filled-shell distribution is, where α = 3 − α c . Since the number of reflected particles is comparatively small, the transmitted PUI distribution can be approximated by the filled-shell distribution with number density = n n n p t p p r . Once specularly reflected back upstream, a PUI is accelerated by the motional electric field. Some reflected PUIs can experience multiple reflections at the shock front with the result that those ions acquire considerable energy. Yao et al. (2021Yao et al. ( , 2021b produced quasi-perpendicular magnetized collisionless shocks in the laboratory. They found that protons from ambient gas could be energized to reach considerable energies of hundreds of keV, and that shock surfing acceleration can be regarded as the sole mechanism for accelerating these ions. The multiple reflection spectra for an initial PUI power-law distribution can be approximated as (Zank et al. 1996;Lipatov et al. 1998;Rice et al. 2001) where η is the power-law index, v 0 is the reference speed, and c mp is the maximum speed. In this work, we take v 0 = u 1 /2 = 160 km s −1 . The maximum particle energy is determined by the balance of the particle's Lorentz force and the electric field force from the electrostatic cross-shock potential (Zank et al. 1996. The estimated maximum kinetic energy for a reflected PUI is determined by the ratio between the thickness of the shock ramp (L ramp ) and the gyroradius of a PUI in the shock upstream (r g1 ), The maximum reflected PUI speed in the solar wind frame c mp is v mp − u 1 . The thickness of the ramp of the termination shock TS-3 observed by Voyager 2 is of the order of 6000 km. The observed TS-3 crossing has been studied extensively in the literature because of its very clear and characteristic, almost classical, quasi-perpendicular structure with a clearly defined foot, ramp, and magnetic field overshoot. This particular crossing therefore provides relatively clear estimates of the various shock structure scalings. By contrast, the other crossings that were observed appear to be during times of, possibly shock reformation or shock rippling, all suggestions that have been made in various works. (e.g., Yang et al. 2015;Lembège & Yang 2016 discussed both the (Zank et al. 1996 analysis and the nonstationarity of the termination shock in detail.) For a PUI upstream of the HTS, the PUI gyroradius is ∼63,000 km. This results in an estimated maximum energy for reflected PUIs of ∼15 keV, but this is, of course, uncertain. The maximum energy can be increased by the fine structure of the shock and the turbulence spectrum (Gedalin 1997;Lipatov et al. 1998;Burrows et al. 2010).
To estimate the speed of particles transmitted across the HTS (v 2 ), we must account for the deceleration by the cross-shock potential, thus . 13 We use a Taylor expansion to express Equation (13) in terms of the upstream flow speed u 1 , . Substituting this expression into Equation (14), we get Thus, the speed of a transmitted PUI particle in the solar wind frame is

Diffusive Shock Accelerated Particles
Part of the multiply reflected PUI population could be further accelerated by the diffusive shock acceleration process provided that the particle speed is above a threshold or injection energy. The injection problem in the theories of shock acceleration is still under debate (Giacalone & Jokipii 1999;Zank et al. 2006Zank et al. , 2001Giacalone 2012;Giacalone et al. 2022;Perri et al. 2022). Parker et al. (2014) argued that the injection energy must be at least a few times the upstream thermal energy so that the particle can make an initial crossing of the shock boundary. The upstream and downstream core solar wind temperatures at the HTS are ∼2 × 10 4 K and ∼2 × 10 5 K, respectively (Richardson & Wang 2010). Hence the thermal speed of solar wind protons is only a few tens km s −1 . Following Neergaard Parker & Zank (2012), Parker et al. (2014), we take the threshold speed c inj as a free parameter.
The downstream accelerated particle distribution function f (c) is given by (Parker et al. 2014;Zank 2014), where u 1(2) is the upstream (downstream) flow velocity, q = 3r/(r − 1) is the power-law index for particles with energies exceeding the maximum energy of seed particles, ) is a pre-existing high-energy particle distribution function far upstream of the shock, p ¢ ¢ Q c c 4 2 ( ) ( ) is the injected particle population at the shock, usually modeled as a delta function. In this work, we only consider an upstream background population.
As discussed above, the pre-existing high-energy particles could be the suprathermal tail of PUI protons and multiply reflected PUIs. These particles could have considerable energy, and thus have the possibility of being further accelerated via diffusive shock acceleration. Both the multiply reflected particles and the PUI high-energy tail have power-law form distribution functions f ( − ∞ , c) ∝ c − γ in the speed range (c l , c m ). The spectrum for the accelerated particles given by Equation (19) is, under the assumption that γ ≠ q. For the special case, γ = q, the spectrum follows Thus, for a source spectrum that is softer than c − q (γ > q), the accelerated spectrum will resemble a single power law with index q. But if the source has a spectrum harder than c − q (γ < q), the accelerated particle spectrum will have a broken power-law form with the break occurring at the maximum speed of the source c m . For c < c m , the accelerated spectrum has approximately the same slope as the source. However, for c > c m , the spectral index for accelerated particles is q. Explicitly, the downstream distribution function is given by The downstream distribution function of transmitted particles is determined by Equation (18), and the maximum speed is rc inj . In total, there are five possible types of high-energy protons (>∼1 keV) in the heliosheath: transmitted PUIs ( f p t ), multiply reflected PUI that are transmitted ( f p r t , ), multiply reflected PUIs that are accelerated via diffusive shock acceleration ( f p r a , ), transmitted high-energy PUI tail ( f tail t ), and a shock accelerated high-energy tail ( f tail a ). These distribution functions just behind the HTS have the following forms,

Distribution Functions in the Heliosheath
The distribution function f in the heliosheath is obtained by solving the Parker transport equation in the form provided that diffusion in spatial and velocity space in the IHS are negligible (Fahr & Lay 2000;Fahr & Scherer 2004). Here u is the bulk solar wind velocity, β = (n H σ cx,H + n He σ cx,He )v rel is the loss rate due to charge exchange between protons and neutral hydrogen and helium. σ cx,H(He) is the charge-exchange cross section between the proton and neutral hydrogen (helium). The analytical formulae for the cross sections are discussed in Appendix B. For simplification, we approximate the relative speed v rel by the proton speed v. The first term on the right hand of Equation (30) describes the de/acceleration of protons due to the divergence of the large-scale flow. For the power-law distribution function f ∝ c − γ( x) , this term can be written as (Zirnstein & McComas 2015;Schwadron & Bzowski 2018 ( ) · is the timescale for adiabatic heating/ cooling. The acceleration timescale is of the same order as the charge-exchange timescale in the heliosheath τ ex = 1/β (Zirnstein & McComas 2015). Using observations from Voyager 2 from the years 2008 to 2018, Schwadron & Bzowski (Schwadron 2018) found that the acceleration ratedominates the charge-exchange rate for more than 70% of the time period studied. Thus the adiabatic heating/cooling term should not be ignored.
Equation (30) can be solved by the method of characteristics (see e.g., Baliukin et al. 2020). The characteristic particle trajectory and the change in particle velocity are Using the above two equations, the evolution of particle speed is related to the plasma flow speed through = -dc dl c du dl u 3 .

( )
If du/dl < 0 (du/dl > 0), then adiabatic heating/cooling occurs in the heliosheath along the streamline. The particle speed along the streamline is calculated as Equation (30) can be written as df/dt = − βf, and the solution is given by where the exponential term describes the extinction of particles along the path due to charge exchange. In this work, we assume that the solar wind speed in the heliosheath linearly decreases from the HTS to the heliopause (Czechowski et al. 2005) according to where L s characterizes the gradient of the plasma flow speed. We set L s = 50 au. The thickness of the heliosheath L is ∼35 au along the Voyager 2 direction but can vary at different locations (see e.g., Zank 2015). Equation (36) is a very idealized description and may be a good approximation only near the nose of the heliosphere. Also, the flow streamline near the Voyager 2 direction is not radial meaning that plasma originating from multiple locations on the HTS could contribute to the production of ENAs. Figure 2. The ENA flux in the inertial frame resulting from different assumptions about the proton distribution functions. Here "MR" means multiply reflected and "DSA" means diffusive shock acceleration. Here "Tail" refers to the tail of the PUI distribution function, not the heliotail direction.

ENA Flux
The expression for calculating the ENA flux in the inertial frame is given by (see and v = v(c l , l u ), respectively. We set the total proton density in the upstream to be 0.0013 cm −3 , with 25% being PUIs . Figure 2 shows the ENA flux at 1 au in the inertial frame for different proton distributions in the heliosheath. The measurements of ENA from IBEX-Hi and INCA are over the years 2009-2012 along the Voyager 2 direction (McComas et al. 2020). We use the ENA flux data from IBEX-Hi Data Release 16, which is survival probability and Compton-Getting corrected. We also include the HSTOF data between 1996 and 1997 for ENAs that originate within ±45°from the apex direction (Czechowski et al. 2008). Figure 2(a) illustrates the ENA flux from the transmitted PUI population. We set the lower limit of the PUI velocity as c l = 0.1u 1 and V spec = 0.25u 1 for illustration. As noted before by Desai et al. (2014), Zirnstein et al. (2018), the ENA flux from transmitted PUIs far exceeds the observed flux when the filled-shell distribution is used for downstream transmitted PUIs. For the cooling index α c = 2.5, which implies the existence of extra heating for PUIs along the path to the outer heliosphere, there is insufficient ENA flux at low energies while the excess problem becomes more severe for the high-energy flux. It seems that the IBEX-Hi observations do not support the extrapolated large cooling index. The contradiction needs to be further investigated. The dashed-dotted line corresponds to the model for which there is no adiabatic heating in the heliosheath. For this case, we set the plasma flow speed in the heliosheath to a constant 50 km s −1 . The difference between models with and without adiabatic heating is significant. The ENA flux from models with adiabatic heating can extend to energies a few times higher than the PUI cutoff energy. And the ENA spectra at the PUI cutoff energy do not have a sharp cutoff structure. These properties seem not to have been addressed in the previous investigations. Figure 2(b) shows the ENA flux from multiply reflected PUIs that are then directly transmitted downstream of the HTS. For the model without adiabatic heating, the ENA flux approaches a power law in the energy range [rv 0 , rc inj ]. The upper limit is determined by the injection velocity. Increasing c inj could extend the power-law spectrum to higher energy. If c inj is set to the lower limit of the speed of multiply reflected particles (v 0 ), there are no transmitted protons since all experience diffusive shock acceleration. With adiabatic heating, the maximum ENA particle speed near the heliopause is + u L rc u u L inj 2 1 3 ( ) ( ( )) . The flux gradually increases (decreases) above the lower (upper) limit of the proton energy range immediately downstream of the HTS.
The multiply reflected PUIs form a power-law distribution with spectral index η in the energy range [v 0 , c mp ]. Protons with energy smaller than the injection speed c inj are transmitted downstream across the HTS. For protons with energies larger than the injection speed c inj , we assume they are further accelerated at the HTS via diffusive shock acceleration. The resulting ENA spectra are shown in Figure 2(c). The accelerated proton spectral index not only depends on the properties of the preaccelerated protons (η, and c mp ), but also on the spectral index q( = 3r/(r − 1)) determined by the shock compression ratio r. In this work, we set r = 2.5, q = 5, and c mp = 5u 1 . The solid line and dashed line show the ENA flux has a spectral break around c mp due to the multiply reflected PUI spectrum having a power-law spectrum with η = 4 < q. For η > q, the spectrum resembles a single power law with spectral index q. Setting c inj to ∞, ensuring that there are no accelerated protons, implies that all multiply reflected PUIs are transmitted.
The suprathermal PUI tail also follows a power-law spectrum with spectral index ν in the energy range [u 1 , c mt ] upstream of the HTS. In this work, we set c mt = 10u 1 , and the density of the tail to be 1/400 of the PUI density which is in the range of measured ratio that is between a few tens to thousand in the inner heliosphere (see e.g., Zhang et al. 2019). The ENA flux resulting from the transmitted tail is shown in Figure 2(d). Similar to the multiply reflected PUIs, the suprathermal PUI tail with speeds above c inj is further energized by diffusive shock Figure 3. The velocity between different reference frames. Here, ṽ is the ENA velocity in the spacecraft reference frame, u s is the spacecraft velocity, u p is the plasma flow velocity, ¢ v is the ENA or parent particle velocity in the plasma frame.
acceleration. This is shown in Figure 2(e). The properties are similar to the accelerated multiply reflected PUIs.

Conclusion
Interpretation of ENA observations requires a detailed understanding of the energetic proton distributions and their transport process in the heliosheath. This in turn provides considerable insight into the kinetic structure, properties, and physical processes of the HTS and in the heliosheath. In this work, we investigate various energetic proton populations in the heliosheath and illustrate the resulting ENA fluxes in the range of ∼keV to 100 keV. Five assumptions for the proton distribution function are considered: (1) PUIs transmitted downstream of the shock without experiencing reflection or diffusive shock acceleration, (2) multiply reflected PUIs that are transmitted downstream and without being diffusively shock accelerated, (3) multiply reflected PUIs with speeds above the injection speed and thus further energized via diffusive shock acceleration, (4) the PUI high-energy tail that is transmitted downstream of the HTS, (5) and the PUI highenergy tail with particle speeds above the injection speed that are energized via diffusive shock acceleration. All these particles have a power-law form just downstream of the HTS.
Due to adiabatic heating, protons change their energy along the streamline from the HTS to the heliopause. Because the adiabatic acceleration timescale and the charge-exchange timescale are of the same order, adiabatic acceleration generally cannot be ignored in modeling the proton distributions in the heliosheath and hence the computed ENA fluxes. The resulting particle distribution function has a relatively simple form related to that immediately downstream of HTS and the plasma flow speed in the heliosheath. We take into account the influence of adiabatic heating by adopting a simple linear relationship between plasma flow speed and the distance from the HTS. This oversimplifies the adiabatic heating rate, but, it is reasonable to address this effect. In this work, we do not present a parametric study but illustrate the various possibilities. Detailed application and comparison to observations are presented in Kornbleuth et al. (2022).
We illustrate the ENA flux resulting from five different proton populations. The computation is restricted to the ENA flux from the Voyager 2 direction. For models with adiabatic heating, the ENA flux produced by the transmitted PUIs is higher than that observed at energies above the PUI cutoff energy. This feature could provide a constraint on the adiabatic heating rate in the heliosheath. However, a detailed model for the plasma flow speed along the streamline is needed. The extrapolation of the cooling index for PUIs at the HTS is almost twice that of the adiabatic cooling index. This results in a very hard PUI spectrum, and thus it produces a too-low (high) ENA intensity at low (high) energy to match the observations. The relative contribution of high-energy transmitted protons (multiply reflected PUIs and a high-energy tail) and protons further energized via diffusive shock acceleration at the HTS is determined by the threshold speed for injection. For a lower threshold speed, there are fewer transmitted particles. With a high threshold speed, more particles are transmitted downstream rather than accelerated. The high-energy ENAs could originate from the combinations of these energetic protons. Although all these particles have a power-law form at immediately downstream of the HTS, adiabatic heating can modify the shape of spectra, especially around the lower and upper limit of the proton energy at the HTS. Shock accelerated particles could provide the ENA flux at high energies and may be imaged remotely by further ENA experiments. The results obtained in this work will be useful in guiding relevant models investigating the ENA fluxes, especially for the upcoming IMAP mission. We derive a formula to calculate the ENA differential flux. In the inertial frame, the omnidirectional proton distribution function is f (v). The number density of particles with speed around v is For one target particle with a cross section σ, ΔN A particles interact with it in the time interval Δt, giving the rate as Consider a volume element at the position l, with volume ΔV = ΔlΔA and the ENA number density n a . Then, the number of ENAs in the volume is n a ΔV. The total number of particles that interact with these ENAs per unit time Δt is given by where v rel is the relative speed between the ion and the neutral atom. For nonrelativistic particles, we have D = The observed ENA flux can thus be written as where β = n A σv rel is the ionization rate, and J is the differential flux of the ions. Since the distribution function is invariant in a different reference frame, we can substitute f (v) in the inertial frame by ¢ ¢ f v ( ) in the plasma reference frame. Here, ¢ =v v u is the particle velocity in the plasma frame, where u is the plasma flow velocity. In general, the proton velocity in the heliosheath is a function of position and initial velocity at the termination shock, and may be written as If there is only adiabatic heating/cooling, ξ = 1/3. To get the ENA flux for a specified energy, the protons which produce ENAs should have the same velocity (energy) at every point along the path. Assuming that the proton velocity is in the range of ¢ ¢ v v , l m [ ] immediately downstream of the HTS, the ENA particle velocity is in the range of The ENA flux is given by where the integration lower limit and upper limit, l d and l u , are determined by the ENA velocity v and satisfy v = v(v m , l d ) and v = v(v l , l u ).
The ENA flux J ENA in the inertial frame can be written as The differential function of ENA in a spacecraft/observer frame f ENAo is also identical to the distribution function in the inertial frame f ENA , Finally, the differential ENA flux in the spacecraft/observer frame is given as (see also Zirnstein et al. 2013Zirnstein et al. , 2021 The primed variables are in the plasma frame, and the tilde variables are in the spacecraft frame. ṽ is the particle in the spacecraft frame, u s is the velocity of the spacecraft. The relative velocity between the plasma frame and spacecraft frame is u rel = u − u s . These speeds are shown in Figure 3  where α is the angle between the ENA velocity and u rel , and f is the angle between the plasma flow speed and u rel .

Appendix B Charge-exchange Cross Sections
The charge-exchange cross section for proton and hydrogen given by Lindsay & Stebbings (2005) is widely adopted to calculate the ENA flux. They presented a critical review of the published experimental measurements and provided the parameterized cross sections by fitting the experimental data using empirical equations. The parameterized cross section between a proton and a hydrogen atom is given as ( where = E mv 1 2 rel 2 is the collision energy. v rel = v − u H is the relative velocity of an ENA or a parent proton in the inertial frame with a background hydrogen neutral atom moving with the bulk speed u H . The formula is valid for an energy range from 0.005 to 250 keV with estimated absolute accuracy of ±10%. The cross section data complied by Barnett et al. (1990; Ba90; data can be accessed through the websites 6 ) are consistent with the formula from Lindsay & Stebbings (2005) for collision speeds greater than 300 km s −1 . However, the differences increase with a decreasing energy from ∼15% at 200 kms −1 (200 eV) to ∼50% at 10 km s −1 (0.5 eV; Bzowski & Heerikhuisen 2020). Bzowski & Heerikhuisen (2020) used the recommended charge-exchange cross section compiled by Barnett et al. (1990), and developed a new expression valid for 10 −4 < E < 1 keV (BH20), We show these two formulae in Figure 4 in the energy range 10 −4 < E < 1 keV together with the data recommended by Barnett et al. (1990). Based on the recommended cross section for a proton and neutral helium by Barnett et al. (1990)