Strength of the Termination Shock Inferred from the Globally Distributed Energetic Neutral Atom Flux from IBEX

Bishwas L. Shrestha; Eric J. Zirnstein; Jacob Heerikhuisen; Gary P. Zank

doi:10.3847/1538-4365/abf659

1. Introduction

Supersonic solar wind (SW) plasma emitted from the Sun travels outward, interacts with the local interstellar medium (LISM), and forms the heliosphere. The Sun and our entire solar system are moving through the interstellar medium at a relative speed of 25.4 km s⁻¹ (McComas et al. 2015). On the journey of the SW plasma from the Sun toward the LISM, the flow transitions from a supersonic to a subsonic flow at the heliospheric termination shock (HTS). Charge-exchange interaction between plasma ions with neutral atoms coming from the LISM creates a hot population of protons called pickup ions (PUIs) that are picked up by the motional electric field and comove with the thermal core SW plasma (Vasyliunas & Siscoe 1976; Zank 1999). Beyond the HTS, the SW flow deflects away from the interstellar plasma flow. The two opposing flows are separated by a tangential discontinuity called the heliopause (HP). The region between the HTS and HP is called the inner heliosheath (IHS), characterized by a higher temperature and density and a lower flow speed compared to the supersonic SW. The region beyond the HP is the LISM that consists of higher plasma and neutral hydrogen density compared to the outer heliosphere. Extensive reviews of the interaction of the SW with the LISM and the structure of the heliosphere can be found in Zank (1999, 2015).

Voyager 2 (V2) crossed the HTS in 2007 (Richardson et al. 2008; Stone et al. 2008) at a heliocentric distance of 84 au and provided the first in situ measurements of the plasma properties at the termination shock. Unfortunately, the Voyager 1 plasma instrument was not working when it crossed the HTS in 2004 (Stone et al. 2005). Measurements by V2 show that the SW plasma bulk flow speed decreased from ∼400 to ∼125 km s⁻¹ in two steps; it first decreased from ∼400 km s⁻¹ at about 1 au before the shock to ∼320 km s⁻¹ just upstream of the shock (in the form of an extended ∼0.3–1 au foreshock; Decker et al. 2008; Florinski et al. 2009) and then to ∼125 km s⁻¹ just downstream of the shock (Richardson et al. 2008). In contrast, the temperature of the thermalized SW increased from about 20,000 to 180,000 K (Richardson 2008; Richardson et al. 2008). The rather small increase in SW temperature indicates that the thermalized SW downstream of the HTS contains only about 20% of the upstream SW flow energy, and the extra energy is believed to be carried by the PUIs created in the supersonic SW (Zank et al. 1996, 2010; Richardson et al. 2008). Some of these PUIs are heated preferentially at the HTS due to reflection from electrostatic cross-shock potential (Lee et al. 1996; Zank et al. 1996; Kumar et al. 2018) and referred to as reflected PUIs. Though V2 was not able to measure the PUI properties downstream of the HTS, the idea of preferential heating of the PUIs at the shock is supported by measurement of PUIs at an interplanetary shock by the SWAP instrument on New Horizons (Zirnstein et al. 2018b) at about 34 au from the Sun and recent simulations of this event by Lembége et al. (2020). On the other hand, some PUIs are only heated adiabatically at the HTS (similar to the thermal core SW) and referred to as transmitted PUIs. Thermal core SW and PUIs can charge-exchange with neutral atoms of interstellar origin in the IHS, downstream of the HTS, and generate energetic neutral atoms (ENAs) with keV energies. Some of these ENAs can be detected by NASA's Interstellar Boundary EXplorer (IBEX) spacecraft orbiting the Earth, if their motion is directed toward the inner heliosphere.

IBEX was launched in 2008 (McComas et al. 2009) to explore the global interaction between the SW and interstellar medium (see McComas et al. 2004). Since then, it has provided remote measurements of the SW–LISM interaction by measuring ENAs created in the outer heliosphere. These ENA measurements are presented in the form of an all-sky map every 6 months. It has two ENA sensors, IBEX-Lo and IBEX-Hi, each with an angular resolution of 6°. IBEX-Hi (Funsten et al. 2009) measures ENAs in six energy bins from ∼0.5 to 6 keV (the data from the lowest energy bin, with a central energy of 0.45 keV, are neglected due to excessive noise present in the data). IBEX-Hi measurements show a distinct feature in the all-sky maps: a narrow, circular structure that shows an enhancement in the ENA fluxes over the surrounding globally distributed fluxes (GDF). This feature is called the ribbon and was one of the unexpected discoveries of the IBEX mission in 2009. The ribbon fluxes are believed to be created through the so-called secondary ENA mechanism (McComas et al. 2009; Heerikhuisen et al. 2010; Zirnstein et al. 2015; McComas et al. 2017; Zirnstein et al. 2018a) outside the HP, which causes the ribbon to be aligned perpendicular to the draped interstellar magnetic field. The GDF, on the other hand, is primarily created from the charge exchange of PUIs in the IHS. IBEX has measured ENAs over a full solar cycle (solar cycle 24) covering 11 yr from 2009 to 2019 (McComas et al. 2020). This unprecedented measurement by IBEX will be followed by the Interstellar Mapping and Acceleration Probe (IMAP; McComas et al. 2018) mission with improved energy, angular, and temporal resolution.

The core SW and hot suprathermal PUIs are not in thermal equilibrium in the IHS and a multicomponent plasma description is required to model the total proton populations in the IHS. Heerikhuisen et al. (2008) and Prested et al. (2008) simulated the heliosheath ENA fluxes assuming a κ distribution (Livadiotis & McComas 2013) for heliosheath protons. The κ distribution function consists of both a Maxwellian-like core and power law–like tail, and it reduces to a Maxwellian for large values of κ. One drawback of using a κ distribution for the IHS protons is that the assumed value of κ fixes the ratio between the core and tail number densities. Thus, one cannot independently change the characteristics of the core without making a similar change to the tail of the distribution. Zank et al. (2010, hereafter Z10) developed a model for the heliosheath plasma with a superposition of three Maxwellian distributions corresponding to three different populations of protons in the IHS (thermal SW protons and transmitted and shock-reflected PUIs) to simulate the heliosheath ENA fluxes. Zirnstein et al. (2014) extended this model by considering the energy-dependent extinction of protons by charge exchange as they advect with the bulk plasma flow through the IHS. Their model also incorporates the fourth population of PUIs (injected PUIs) created in the IHS by charge exchange of any three of the IHS proton populations with the interstellar neutrals. Heerikhuisen et al. (2019) ran simulations where they included the self-consistent feedback on the plasma of charge exchange between neutral hydrogen and the four proton populations. They found that different model assumptions resulted in different rates for the cooling of the IHS plasma. Shrestha et al. (2020) used the four population formulations for IHS protons to quantitatively compare the simulated IHS ENA fluxes with observed ENA flux from IBEX-Hi and demonstrated that the injected PUI populations are crucial to reproduce the IBEX-Hi data quantitatively. Kornbleuth et al. (2018, 2020) also used the three population descriptions for protons in the IHS to simulate the global ENA fluxes using a croissant-like heliosphere (Opher et al. 2015, 2020) as a background. Baliukin et al. (2020) processed a PUI filled-shell distribution across the HTS assuming conservation of magnetic moment, yielding a single, compressed filled shell of PUIs downstream of the HTS to simulate the IHS ENA fluxes in the framework of the Moscow group model of the heliosphere (Izmodenov & Alexashov 2015, 2020).

In this paper, we estimate the HTS compression ratio at multiple directions in the sky using a quantitative comparison of the ENA flux observed by IBEX-Hi with simulated ENA fluxes from the IHS. We first present a method to estimate the HTS compression ratio and then validate this method by a comparison of the compression ratio estimation in the V2 direction with V2 observations. Finally, we use this method to estimate the shock compression ratio in other directions of the sky (close to the ecliptic plane).

2. Simulation Method

The simulation of the hydrogen ENA fluxes from the IHS consists of two steps. First, we simulate the SW–LISM interaction (Section 2.1), which provides a background heliosphere for the ENA flux calculations. Second, we postprocess the simulation result of the background heliosphere (Section 2.2) to calculate the IHS ENA fluxes at 1 au for comparison to IBEX data.

2.1. Simulation of the Background Heliosphere

We simulate the SW–LISM interaction using a 3D steady-state MHD plasma/kinetic neutral code based on MS-FLUKSS (Pogorelov et al. 2011), which couples an MHD code for plasma (Pogorelov et al. 2006, 2009) to a Monte Carlo particle code for neutral hydrogen (Heerikhuisen et al. 2006, 2008; Heerikhuisen & Pogorelov 2010). The flow of plasma at all points of the SW–LISM domain is described by a single fluid solved using ideal MHD equations. The evolution of neutral hydrogen is described by a (6+1)D Boltzmann equation that tracks particles on ballistic trajectories. The neutral code is coupled to the MHD code by charge-exchange source terms in momentum and energy, describing the charge exchange between a neutral hydrogen atom distribution and protons. In the case of ionization by photons, a mass source term is included based on an ionization rate of 8 × 10⁻⁸ s⁻¹ at 1 au (Bzowski 2008), which scales as r⁻² with distance from the Sun.

The PUI distribution in the IHS is described by a multicomponent PUI model with three different populations of protons in the IHS (Z10; Zirnstein et al. 2013). This model partitions the total thermal energy of protons obtained from MHD into three separate Maxwellian distributions: core SW, transmitted PUIs, and reflected PUIs. See Heerikhuisen et al. (2019) for the details of the number density and total thermal energy partitioning in the IHS. The MHD plasma, kinetic neutral, and PUI modules are repeated iteratively to obtain a steady-state heliosphere, the result of which is collected in a spherical grid to calculate the hydrogen ENA flux at 1 au.

The SW conditions at 1 au are the same as the uniform slow SW case of Shrestha et al. (2020): plasma density 5.74 cm⁻³, bulk plasma flow speed 450 km s⁻¹, total plasma temperature 51,000 K, and radial component of magnetic field 37.5 μG. We choose the uniform slow SW conditions at the simulation inner boundary because the analysis presented in this paper is limited close to the ecliptic plane where the SW speed is always slow. The SW dynamic pressure ( $\tfrac{1}{2}{m}_{{\rm{p}}}{n}_{{\rm{p}}}{v}_{{\rm{p}}}^{2}$ ) used here also roughly corresponds to the average SW dynamic pressure in the ecliptic plane that IBEX had observed over 2009–2011 (see the shaded region in Figure 1). Note that it takes around 2–4 yr for the SW condition at 1 au to be observed by IBEX as ENAs in the upwind part of the heliosphere, whereas it takes more than 6 yr in the tail direction. Then the SW values at the simulation inner boundary (10 au) are obtained by following the 1 au values advected by simple adiabatic expansion to 10 au. The LISM boundary condition at 1000 au is the same as the 3 μG case of Zirnstein et al. (2016): plasma density 0.09 cm⁻³ and neutral hydrogen density 0.154 cm⁻³, both at a temperature of 7500 K. The flow speed for both species is 25.4 km s⁻¹ (McComas et al. 2015) coming from (255 fdg 7, 5 fdg 1) in ecliptic J2000 coordinates, and the magnetic field with strength 3 μG is directed at (226 fdg 99, 34 fdg 81) from ecliptic J2000 coordinates.

**Figure 1.** Variation of the SW dynamic pressure with time in the ecliptic plane at 1 au from the OMNI data set. The dark and light gray regions roughly correspond to the average SW condition that IBEX had observed over 2009–2011 in the upwind and downwind regions of the heliosphere, respectively.
Download figure:
Standard image High-resolution image

2.2. IHS ENA Flux Calculation

We simulate the IHS ENA flux in the framework of a multi-Maxwellian description for proton distribution at all points in the IHS (Z10; Zirnstein et al. 2014, 2017). The total proton distribution is constructed as a sum of three Maxwellian distributions: (i) core SW, (ii) transmitted PUI, and (iii) reflected PUIs. We also consider the energy-dependent extinction of these protons by charge exchange, as they are advected with the bulk plasma flow through the IHS. Details of the IHS ENA flux calculation are described in Sections 2.2.1, 2.2.2, and 2.2.3 below.

2.2.1. PUI Properties Downstream of the HTS

In the first step, we partition the total number density immediately downstream of the HTS obtained from MHD into three separate Maxwellian distributions as

$\begin{eqnarray}&&{n}_{\mathrm{MHD}}={n}_{{\rm{p}}}^{\mathrm{core}}+{n}_{{\rm{p}}}^{\mathrm{tr}}+{n}_{{\rm{p}}}^{\mathrm{ref}},\end{eqnarray} \tag{ 1 }$

where ${n}_{{\rm{p}}}^{\mathrm{core}}$ , ${n}_{{\rm{p}}}^{\mathrm{tr}}$ , ${n}_{{\rm{p}}}^{\mathrm{ref}}$ , and n_MHD are the densities of the core SW, transmitted PUIs, reflected PUIs, and total protons immediately downstream of the HTS, respectively. The number density of PUIs ( ${n}_{\mathrm{PUI}}={n}_{{\rm{p}}}^{\mathrm{tr}}+{n}_{{\rm{p}}}^{\mathrm{ref}}$ ) just upstream of the HTS at a direction (θ, ϕ) in the sky is obtained from the simulation by integrating the production of PUIs along a radial direction by charge exchange, photoionization, and electron impact ionization of neutral hydrogen (see Equation (6) in Heerikhuisen et al. 2019). The PUI density fraction at the HTS is defined as a ratio of the total number density of PUIs just downstream of the HTS to the number density obtained from the MHD model (i.e., α = n_PUI/n_MHD).

If we assume that PUIs follow a filled-shell distribution (Vasyliunas & Siscoe 1976) in the SW upstream of the HTS, then the density fraction of reflected PUIs (β), which is defined as the ratio of reflected PUI density to total PUI density (i.e., $\beta ={n}_{{\rm{p}}}^{\mathrm{ref}}/{n}_{\mathrm{PUI}}$ ), can be written as

$\begin{eqnarray}\beta =\left\{\begin{array}{ll}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}-{\left(1-\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}\right)}^{3/2}\left(\tfrac{1}{\sqrt{1-\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}}}-1\right)\right], & \mathrm{if}\ {V}_{\mathrm{spec}}\leqslant {\text{}}{u}_{1}\\ \left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}+{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}-1\right)}^{3/2}\left(\displaystyle \frac{1}{\sqrt{\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}-1}}-1\right)\right], & \mathrm{if}\ {V}_{\mathrm{spec}}\gt {\text{}}{u}_{1},\end{array}\right.\end{eqnarray} \tag{ 2 }$

where u₁ is the upstream plasma flow speed (also the PUI cutoff speed), and V_spec is the specular reflection velocity for protons in the shock rest frame (i.e., protons with v_x < V_spec in the shock rest frame experience specular reflection). See the Appendix for a detailed derivation of Equation (2). The value of the specular reflection velocity can be estimated using the following expression (Zank et al. 1996):

$\begin{eqnarray}&&{V}_{\mathrm{spec}}^{2}=\displaystyle \frac{2\eta (r-1)}{{M}_{{\rm{A}}1}^{2}}{u}_{1}^{2},\end{eqnarray} \tag{ 3 }$

where η is empirically found to be about 2, r is the shock compression ratio, and M_A1 is the upstream Alfvénic Mach number (≈8, V2 observation; Richardson et al. 2008). The value of β depends on the shock compression ratio (r), and for a shock compression ratio of 2.5 and upstream flow speed (u₁) of 320 km s⁻¹, the reflected PUI density fraction becomes about 0.04. This estimate is lower than that obtained in Z10, which is ∼0.1 (note that there is a typo in Equation (11) in Z10). Finally, the densities of the different populations of protons in terms of α and β can be expressed as

$\begin{eqnarray}&&{n}_{{\rm{p}}}^{{\rm{i}}}={{\rm{\Pi }}}_{{\rm{p}}}^{{\rm{i}}}{n}_{\mathrm{MHD}},\end{eqnarray} \tag{ 4 }$

where ${{\rm{\Pi }}}_{{\rm{p}}}^{\mathrm{core}}=(1-\alpha )$ , ${{\rm{\Pi }}}_{{\rm{p}}}^{\mathrm{tr}}=\alpha \,(1-\beta )$ , and ${{\rm{\Pi }}}_{{\rm{p}}}^{\mathrm{ref}}=\alpha \,\beta$ are density fractions of the core SW, transmitted PUIs, and reflected PUIs, respectively.

We assume that the temperature of the core SW just downstream of the HTS in all directions is similar to the V2 measurements, i.e., ${T}_{{\rm{p}}}^{\mathrm{core}}={\rm{181,000}}\,{\rm{K}}$ , though the exact value has little impact on our results. The reason for this is that due to their low temperature, core SW ions do not produce ENAs within the IBEX-Hi energy range and therefore do not affect our results. The temperature of the transmitted PUIs immediately downstream of the HTS is defined in terms of the shock compression ratio as (Z10)

$\begin{eqnarray}&&{T}_{{\rm{p}}}^{\mathrm{tr}}={r}^{2}{T}_{{\rm{p}},\ 1}^{\mathrm{PUI}},\end{eqnarray} \tag{ 5 }$

where ${T}_{{\rm{p}},1}^{\mathrm{PUI}}$ is the PUI temperature upstream of the HTS. The value of ${T}_{{\rm{p}},1}^{\mathrm{PUI}}$ can be obtained by taking the scalar pressure moment of the PUI distribution (assumed to be a filled shell)

$\begin{eqnarray}&&{T}_{{\rm{p}},1}^{\mathrm{PUI}}=\displaystyle \frac{1}{{n}_{\mathrm{PUI}}}\displaystyle \frac{{m}_{{\rm{p}}}}{3{k}_{{\rm{B}}}}{\int }_{0}^{{u}_{1}}f({\bf{c}}){c}^{2}{d}^{3}c,\end{eqnarray} \tag{ 6 }$

where m_p is the proton mass, n_PUI is the PUI number density, c is the particle speed in the SW frame, and u₁ is the cutoff speed of the filled-shell distribution in the SW frame. For a bulk flow speed (u₁) of 320 km s⁻¹, the upstream PUI temperature becomes 1.77 × 10⁶ K. To obtain the upstream PUI temperature in other directions, we first rescale the upstream flow speed there (from MHD) to the observed value at the V2 direction (see Equation (19)) and then use the rescaled flow speed in Equation (6). The filled-shell distribution for PUIs in the SW frame in the limit that the distance from the Sun is much larger than the ionization cavity (Vasyliunas & Siscoe Vasyliunas & Siscoe 1976; Z10) is given by

$\begin{eqnarray}&&f({\boldsymbol{c}})=\displaystyle \frac{3{n}_{\mathrm{PUI}}}{8\pi {u}_{1}^{3/2}}{c}^{-3/2}{\rm{\Theta }}({u}_{1}-c),\end{eqnarray} \tag{ 7 }$

where Θ(x) is the Heaviside step function. An estimate for the reflected PUI temperature downstream of the HTS is provided in Z10 as

$\begin{eqnarray}&&{T}_{{\rm{p}}}^{\mathrm{ref}}=\displaystyle \frac{{m}_{{\rm{i}}}}{3{k}_{{\rm{B}}}}\left(1+\displaystyle \frac{1}{4}{\left({r}_{{\rm{g}}1}/{L}_{\mathrm{ramp}}\right)}^{2}\right){u}_{1}^{2},\end{eqnarray} \tag{ 8 }$

where r_{g 1} is the gyroradius of an upstream PUI, L_ramp is the shock ramp thickness, and u₁ is the upstream plasma flow speed. Using r_{g 1} ∼ 55,000 km and L_ramp ∼ 6000 km as observed by V2 (Burlaga et al. 2008) yields ${T}_{{\rm{p}}}^{\mathrm{ref}}\sim 9.31\times {10}^{7}\,{\rm{K}}$ . To calculate the reflected PUI temperature in other directions, the upstream flow speed there is also rescaled with the observed upstream flow speed at the V2 direction, and we assume that the ratio r_{g 1}/L_ramp does not change.

2.2.2. PUI Distributions in the IHS

As plasma flows away from the HTS along a streamline through the IHS, protons are depleted by charge exchange with hydrogen atoms coming from interstellar space and replaced by a newly injected population of ionized hydrogen atoms. Hence, we update the density of protons in the IHS by taking into account the extinction effect as (Zirnstein et al. 2017)

$\begin{eqnarray}&&\begin{array}{l}{n}_{{\rm{p}}}^{{\rm{i}}}({\boldsymbol{r}})={{\rm{\Pi }}}_{{\rm{p}}}^{{\rm{i}}}\,{n}_{\mathrm{MHD}}({\boldsymbol{r}})\\ \quad \times \exp \left(-{\int }_{{r}_{\mathrm{HTS}}}^{r}{n}_{{\rm{H}}}(r^{\prime} ){\sigma }_{\mathrm{ex}}({v}_{\mathrm{rel}}^{{\rm{i}}}(r^{\prime} )){v}_{\mathrm{rel}}^{{\rm{i}}}(r^{\prime} )\displaystyle \frac{{ds}}{{u}_{{\rm{p}}}(r^{\prime} )}\right),\end{array}\end{eqnarray} \tag{ 9 }$

where the integration is carried out over a flow streamline of segment ${ds}={u}_{{\rm{p}}}(r^{\prime} ){dt}$ , where ${u}_{{\rm{p}}}(r^{\prime} )$ is the bulk plasma flow speed in the IHS, and ${n}_{{\rm{H}}}(r^{\prime} )$ is the neutral hydrogen density along the same streamline. The energy-dependent charge-exchange cross section ${\sigma }_{\mathrm{ex}}({v}_{\mathrm{rel}}^{{\rm{i}}}(r^{\prime} ))$ for each species is obtained from Lindsay & Stebbings (2005), and the relative interaction velocity ${v}_{\mathrm{rel}}^{{\rm{i}}}(r)$ of each proton distribution (Maxwellian) with neutral hydrogen distribution (Maxwellian) is approximated by Pauls et al. (1995),

$\begin{eqnarray}&&\begin{array}{l}{v}_{\mathrm{rel}}^{{\rm{i}}}(r^{\prime} )\\ =\ \sqrt{\tfrac{4}{\pi }\left({\left({v}_{\mathrm{th},{\rm{p}}}^{{\rm{i}}}(r^{\prime} )\right)}^{2}+{\left({v}_{\mathrm{th},{\rm{H}}}(r^{\prime} )\right)}^{2}\right)+| {{\boldsymbol{u}}}_{{\rm{p}}}(r^{\prime} )-{{\boldsymbol{u}}}_{{\rm{H}}}(r^{\prime} ){| }^{2}},\end{array}\end{eqnarray} \tag{ 10 }$

where ${v}_{\mathrm{th},{\rm{p}}}^{{\rm{i}}}(r^{\prime} )$ is the thermal velocity for each proton distribution, and ${v}_{\mathrm{th},{\rm{H}}}(r^{\prime} )$ is the thermal velocity for the local neutral hydrogen distribution along a streamline. Note that it is important for the global simulations to describe neutrals kinetically. However, for the postprocessing, it is sufficient to model only the lower-energy hydrogen atoms, since the ENAs produced in the SW represent on the order of 1% of the total neutral hydrogen density in the heliosheath. Since the temperatures of the LISM and outer heliosheath (OHS) are much smaller than the temperature of the protons in the IHS and the IHS flow speed, combining LISM/OHS neutrals into one Maxwellian is sufficient in Equation (10). The charge-exchange interaction removes high-energy protons and replaces them with low-energy (∼0.05 keV) injected PUIs. Since we do not include any physical mechanism that might heat the injected PUIs to the IBEX-Hi energy range (>0.7 keV), we neglect this population when calculating the IHS ENA fluxes.

To calculate the temperature of the core SW, transmitted PUIs, and reflected PUIs at different locations in IHS, we use their temperature ratios at the HTS with the total proton temperature. These ratios are assumed to remain constant along a streamline throughout the IHS (i.e., ${{\rm{\Gamma }}}_{{\rm{p}}}^{{\rm{i}}}={T}_{{\rm{p}}}^{{\rm{i}}}/{T}_{{\rm{p}}}={T}_{{\rm{p}}}^{{\rm{i}}}({\boldsymbol{r}})/{T}_{{\rm{p}}}({\boldsymbol{r}})$ ). Hence,

$\begin{eqnarray}&&{T}_{{\rm{p}}}^{{\rm{i}}}({\boldsymbol{r}})={{\rm{\Gamma }}}_{{\rm{p}}}^{{\rm{i}}}\times {T}_{{\rm{p}}}({\boldsymbol{r}}),\end{eqnarray} \tag{ 11 }$

where T_p is the total proton temperature at the HTS, and ${{\rm{\Gamma }}}_{{\rm{p}}}^{\mathrm{core}}={T}_{{\rm{p}}}^{\mathrm{core}}/{T}_{{\rm{p}}}$ , ${{\rm{\Gamma }}}_{{\rm{p}}}^{\mathrm{tr}}={T}_{{\rm{p}}}^{\mathrm{tr}}/{T}_{{\rm{p}}}$ , and ${{\rm{\Gamma }}}_{{\rm{p}}}^{\mathrm{ref}}={T}_{{\rm{p}}}^{\mathrm{ref}}/{T}_{{\rm{p}}}$ are the temperature fractions of the core SW, transmitted PUIs, and reflected PUIs, respectively. These fractions are calculated by dividing the temperature of the core SW (V2 observations), transmitted PUIs (Equation (5)), and reflected PUIs (Equation (8)) by the total proton temperature (T_p) immediately downstream of the HTS. Note that these fractions are not constant, as in previous simulations (Zirnstein et al. 2014, 2017; Shrestha et al. 2020), but rather a function of the total proton temperature just downstream of the HTS at any shock location.

2.2.3. Line-of-sight ENA Flux at 1 au

We compute the hydrogen ENA flux at 1 au for a line of sight (LOS) by calculating the charge-exchange rates of all three populations of protons at all points in the IHS. The differential hydrogen flux for an individual population in the inertial frame, assuming a Maxwellian interstellar neutral hydrogen distribution, is given by Zirnstein et al. (2013),

$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{J}^{{\rm{i}}}({\boldsymbol{r}},{\boldsymbol{v}}) & = & \displaystyle \frac{1}{{m}_{{\rm{H}}}}{f}_{{\rm{p}}}^{{\rm{i}}}({\boldsymbol{r}},{{\boldsymbol{v}}}_{{\rm{p}}}){v}^{2}P({\boldsymbol{r}},{\boldsymbol{v}}){n}_{{\rm{H}}}({\boldsymbol{r}}){v}_{\mathrm{rel}}(| {\boldsymbol{v}}-{{\boldsymbol{u}}}_{H}({\boldsymbol{r}})| )\\ & & \times {\sigma }_{\mathrm{ex}}({v}_{\mathrm{rel}}({\boldsymbol{r}})){\rm{\Delta }}t,\end{array}\end{eqnarray} \tag{ 12 }$

where m_H is the mass of the hydrogen atom, ${f}_{{\rm{p}}}^{{\rm{i}}}({\boldsymbol{r}},{{\boldsymbol{v}}}_{{\rm{p}}})$ is the distribution function for different proton species i in the plasma frame, v _p is the parent proton velocity in the plasma frame (i.e., v _p = v − u _p, where u _p is the bulk plasma flow velocity), v is the ENA velocity in the inertial frame, P( r , v ) is the ENA survival probability, n_H( r ) is the background neutral hydrogen density, v_rel is the relative velocity between the parent proton and background neutral hydrogen distribution (Ripken & Fahr 1983; Heerikhuisen et al. 2006), u _H( r ) is the neutral hydrogen bulk velocity, σ_ex is the energy-dependent charge-exchange cross section (Lindsay & Stebbings 2005), and Δt is the integration time step. We integrate Equation (12) from the HTS to the HP along a radial vector (r) for each population and then sum their results to get the total IHS ENA flux along that LOS. The survival probability of hydrogen ENAs is given by P( r , v ) = 1 − ∫β( r , v )dt, where β( r , v ) is the ionization rate of a single hydrogen ENA due to charge exchange and photoionization (we ignore loss due to electron impact ionization). We only consider losses to the hydrogen ENAs beyond 100 au (average radial distance to the HTS from the Sun) because the IBEX measurements are corrected for losses between 1 and 100 au. We also take into account the energy response function of the IBEX-Hi detectors (Funsten et al. 2009) to calculate the differential ENA fluxes in the IBEX-Hi energy range. The energy response data of the IBEX-Hi detector is available on the IBEX public release website at https://ibex.princeton.edu/DataRelease1/CalibrationData.

3. IBEX Data

The data used for the comparison of the simulated ENA flux are from IBEX-Hi over the energy range 0.71–4.29 keV (central energy) observed during the first 3 yr of the mission, from 2009 to 2011 (McComas et al. 2020). The choice of the IBEX data over this period is to ensure that it is as close to the V2 HTS crossing time as possible, while having a few years of IBEX data to average and improve statistics. This data set is available to the public as Data Release 16 at the IBEX website: https://ibex.princeton.edu/DataRelease16. In this data set, the IBEX team has optimized the data processing compared to previous releases, and it covers a full solar cycle (11 yr from 2009 to 2019) of observation. Each 6 month data set used is survival probability corrected between 1 and 100 au (McComas et al. 2014) and Compton–Getting corrected to the solar inertial frame. We average the six data sets (2009A, 2009B, 2010A, 2010B, 2011A, and 2011B) for each 6° × 6° pixel (total of 1800) with a weight proportional to the detector exposure time for that pixel. Note that we use the original data set (not the one with the ribbon subtracted) and avoid the ribbon region in the analysis to make sure that we compare the simulated ENA fluxes directly with IBEX data.

4. Compression Ratio Estimation

We estimate the HTS compression ratio for a direction in the sky using a quantitative comparison of the observed ENA fluxes from IBEX-Hi with the simulated IHS ENA fluxes for a range of compression ratios. To reiterate, our model of the IHS ENA fluxes depends on the HTS compression ratio according to Equations (2)–(5). First, we simulate the IHS ENA fluxes at an LOS of the sky for a range of compression ratios (1.3 ≤ r ≤ 3.5). Then, the simulated fluxes corresponding to each compression ratio are compared with the 9 pixels of IBEX-Hi data centered on that LOS that are weight-averaged using the measurement uncertainties. The comparison is performed by calculating the fractional difference in the spectral slope between the simulated ENA fluxes and IBEX data as

$\begin{eqnarray}&&| {\rm{\Delta }}{\gamma }_{\mathrm{fr}}| =\left|\displaystyle \frac{({\gamma }_{\mathrm{data}}-{\gamma }_{\mathrm{sim}})}{{\gamma }_{\mathrm{data}}}\right|,\end{eqnarray} \tag{ 13 }$

where γ_data and ${\gamma }_{\mathrm{sim}}$ are the spectral slopes of the ENA fluxes from the data and simulation, respectively. The ENA spectral slopes are calculated by fitting a power law J(E) ∝ E^−γ over the IBEX-Hi energy range (0.71–4.29 keV), where J is the ENA flux from data (or simulation) at energy E, and γ is the spectral slope over the fitted energy range. The power-law fit involves a χ² minimization of a linear function fit in logarithmic space that requires a variance of the flux logarithm "ln J " given by ${\sigma }_{\mathrm{ln}\,J}^{2}={({\sigma }_{J}/J)}^{2}$ . We also consider an estimate of the systematic uncertainties as 20% of the ENA flux for IBEX-Hi (Fuselier et al. 2012, 2014) and add them in quadrature to the statistical uncertainties for ENA count rates while calculating the uncertainty (σ_J) in the ENA flux data. However, we cannot assign any uncertainty to the simulated ENA flux, so we compute their spectral slope by assuming equal uncertainty for logarithmic flux ( ${\sigma }_{\mathrm{ln}J}=1$ ) for each energy channel. Finally, the best compression ratio in a direction is given by the shock compression ratio corresponding to the minimum value of ∣Δγ_fr∣.

5. Results and Discussions

5.1. Compression Ratio in the V2 Direction

We simulate the IHS ENA flux in the V2 direction using a 3D steady-state background heliosphere. The background heliosphere is simulated using a uniform slow SW at all latitudes that corresponds to case III of Heerikhuisen et al. (2019), where the electron temperature in the IHS is assumed to be equal to the MHD temperature, i.e., T_e = T_MHD = T_p. The bulk plasma flow speed upstream of the HTS is taken from the V2 observations (320 km s⁻¹) that give an upstream PUI temperature ( ${T}_{{\rm{p}},1}^{\mathrm{PUI}}$ ) of 1.77 × 10⁶ K (Equation (6)). The reflected PUI temperature immediately downstream of the HTS is obtained from the Z10 estimate (Equation (8)), and its value is around 9 × 10⁷ K. The comparison of the simulated IHS ENA fluxes for each shock compression ratio with the time exposure–averaged IBEX-Hi data for the first 3 yr of the mission is shown in Figure 2. The left panel of Figure 2 shows the variation of the fractional difference in spectral slope between IBEX-Hi data and simulated ENA flux for a range of HTS compression ratios (1.3 ≤ r ≤ 3.5). The compression ratio corresponding to the minimum value of ∣Δγ_fr∣ represents the best HTS compression ratio in the V2 direction that matches the spectral slope from the IBEX-Hi data. We obtain a compression ratio of 2.4 for the V2 direction that is very close to the large-scale compression ratio measured by V2 (∼2.5) in its crossing of the HTS in 2008 (Richardson et al. 2008).

**Figure 2.** (Left panel) Variation of the fractional difference in spectral slope between IBEX-Hi data and simulated IHS ENA fluxes for a range of shock compression ratios (1.3 ≤ r ≤ 3.5). (Right panel) ENA spectrum in the V2 direction from both data and simulation using the best shock compression ratio of 2.4. The power-law fit to the ENA spectrum from both data and simulation is also shown by dotted red and blue lines, respectively. The contribution to the simulated fluxes by individual proton species is also shown on the same plot: core SW (dashed black), transmitted PUIs (dashed cyan), and reflected PUIs (dashed green). Simulated ENA fluxes are also shown beyond the IBEX-Hi energy range (shaded region on the right) to see the effect of reflected PUIs on the total ENA flux.
Download figure:
Standard image High-resolution image

The ENA spectra in the V2 direction from both IBEX-Hi data and simulation using the best-fitted compression ratio of 2.4 are plotted in the right panel of Figure 2. A power-law fit to the ENA spectrum from both data and simulation is also shown on the same plot (dotted red and blue lines). As the simulated ENA fluxes are lower by a factor of ∼3 (see also Figure 3) compared to the data in the V2 direction, the comparison of the spectral indices between them serves a better purpose to quantify their differences. A large difference in flux between simulation and data was also reported by many authors (e.g., Zirnstein et al. 2017; Kornbleuth et al. 2020; Shrestha et al. 2020). All of these models include a Maxwellian description for PUIs downstream of the HTS, unlike in Baliukin et al. (2020). The investigation of the exact reason for this discrepancy is a subject of future study. The underestimation can be partially explained by a higher interstellar neutral hydrogen density, as recently reported by Swaczyna et al. (2020). For the purpose of our current study, we assume that the underestimation of the observed ENA fluxes is uniform over energy and would not significantly change our comparison of spectral slopes. Note that the data, and especially the simulation, are not exactly a pure power law in any direction of the sky (see also Figure 9). The contribution to the simulated ENA flux from each proton population is also shown on the same plot to see their effect on the total ENA flux. The contribution to the simulated ENA flux by the core SW (black dashed line) is so low that it is not visible in the scale chosen. It is clear from this figure that most of the ENA flux in the IBEX-Hi energy range is coming from the transmitted PUIs. Also, the reflected PUIs only contribute a significant number of ENAs in the highest energy channel of IBEX-Hi (central energy of 4.29 keV). We have also shown the simulated ENA flux for energies higher than the IBEX-Hi energy range (up to 10 keV) in the right part of the plot (shaded in gray). The ENA flux for the higher energy shows that the contribution of transmitted PUIs to total ENA flux drops sharply after 5 keV, at which point the reflected PUIs begin to dominate. Future observations by IMAP, which will measure ENAs in this energy range and beyond, will be essential to verify this result and to allow for the presence of reflected PUIs to influence comparisons with data.

**Figure 3.** The ENA spectrum in the V2 direction from IBEX-Hi along with simulated ENA spectra for different shock compression ratios. Note that the simulated fluxes are multiplied by a normalization factor equal to the ratio of the observed ENA flux at 1.74 keV to the simulated ENA flux corresponding to the best compression ratio at the same energy.
Download figure:
Standard image High-resolution image

5.1.1. Effect of Compression Ratio on the ENA Spectrum

Simulated ENA spectra for three different shock compression ratios (1.3, 2.4, and 3.5) along with the ENA spectrum from IBEX-Hi data in the V2 direction are shown in Figure 3. The simulated ENA fluxes are multiplied by a normalization factor equal to the ENA flux at 1.74 keV from observations to the simulated ENA flux at the same energy using a compression ratio of 2.4 for a better comparison. Evidently, the ENA fluxes increase with increasing shock compression ratio for higher energy channels. However, for lower energy channels, the ENA fluxes first increase with the compression ratio and then decrease after a certain value of shock compression. This behavior is related to the evolution of a Maxwellian distribution for the transmitted PUIs with an increase in the compression ratio. An increase in the compression ratio increases the temperature of the transmitted PUIs by a factor of r² (see Equation (5)) that results in PUIs gaining energy and populating the higher-energy side of the distribution. This leads to an increase in the number of ENAs on the higher-energy side of IBEX-Hi. At the same time, an increase in temperature flattens the Maxwellian distribution, leaving a lower number of particles on the lower-energy side and hence fewer ENAs around 1 keV.

5.1.2. Differential ENA Flux along the IHS and the Footpoint Locations of the Flow Streamline

The differential ENA fluxes created along each radial step outward from the HTS (the integral of which gives the total simulated ENA flux) in the V2 direction for all five IBEX-Hi energies are shown in Figure 4. An interesting feature in this figure is that the differential ENA flux does not decrease monotonically from the HTS through the IHS but rather reaches a maximum at a certain distance from the HTS and then drops thereafter. Four significant quantities affecting the number of ENAs created at any point in the IHS are (i) bulk plasma flow speed, (ii) density of neutral hydrogen, (iii) densities of individual proton species, and (iv) temperatures of individual proton species. If the first two quantities are constant at all points in the IHS, one would expect that the ENA production decreases almost exponentially from the HTS toward the HP. But this is not true in the heliospheric plasma flow, where the radial plasma flow speed decreases gradually outward from the HTS while the density of neutral hydrogen increases. That increases the charge-exchange rate between protons and neutral hydrogen and hence also increases ENA production. However, the densities of different proton species are depleted by charge exchange as plasma flows away from the HTS, which results in a lower ENA production rate. Finally, the increased tangential flow close to the HP further decreases the observed ENA flux. The net effect is that ENA production first slowly increases and then swiftly decreases monotonically with radial distance through the IHS. The proton temperature can also affect the number of ENAs created in the IHS by increasing (or decreasing) the number of available ions for charge exchanging into ENAs. But we do not see any significant correlation between the temperature and the ENA production rate as a function of distance through the IHS, as represented in Figure 4. The reason for this is likely due to a slower variation of the PUI temperature with distance compared to the flow speed. The cumulative differential fluxes from the HTS in the V2 direction for each IBEX-Hi energy channel are also shown in the same figure.

**Figure 4.** Simulated differential ENA flux along the V2 direction along each radial step outward from the HTS through the IHS that contributes to the total LOS ENA flux at 1 au for each IBEX-Hi energy channel. The cumulative differential flux along the IHS is also plotted with dashed lines. The cumulative ENA fluxes are normalized by the total ENA flux for each IBEX-Hi energy channel.
Download figure:
Standard image High-resolution image

The footpoint locations at the HTS for the plasma flow streamlines from each radial step outward of the HTS along the V2 LOS are plotted in a Mollweide projection and shown in Figure 5. These footpoint locations are obtained by tracing the plasma flow vectors backward along each radial step outward from the HTS until it reaches the HTS at any location. The footpoint locations are color-coded by a fractional ENA flux created from each radial step. The fractional ENA flux along an LOS for each energy channel is calculated by dividing the ENAs produced at a certain radial distance by the maximum number of ENAs produced along that LOS. The average value of fractional ENA fluxes for five energy channels is used to color-code the footpoint locations. The footpoint locations are also binned into 6° × 6° pixels and plotted as semitransparent rectangular bins in the same figure. The HTS bins are also color-coded by the average value of the fractional fluxes for all coordinates lying on that bin. It is clear from the figure that the HTS locations contributing to the plasma distribution along the V2 LOS are very close to that direction, and the footpoint origins of streamlines farther along the V2 LOS connect closer to the nose of the heliosphere. This is a direct result of the antinoseward flow of the IHS plasma away from the IHS pressure maximum.

**Figure 5.** Footpoint locations of plasma flow streamlines at the HTS for plasma flowing along each radial step in the V2 direction plotted in a Mollweide projection. The footpoint locations are color-coded by a fractional ENA flux produced from the plasma flowing from that location to the V2 LOS. The fractional flux used here is the average of the fractional ENA fluxes for five energy channels of IBEX-Hi. The footpoint locations are also binned into a 6° × 6° HTS bin and color-coded by the average fractional flux for that bin.
Download figure:
Standard image High-resolution image

5.1.3. Effect of Reflected PUI Temperature

All of the analyses presented above are done by using Equation (8) for the temperature of the reflected PUIs (Z10 estimate) and assuming that the electron temperature in the IHS is equal to the MHD temperature (T_e = T_MHD). Here we investigate how the shock compression ratio changes for different ways to estimate the reflected PUI temperature. Note that we still estimate the core SW and transmitted PUI temperature immediately downstream of the HTS from V2 observations and Equation (5), respectively. The two other estimates for the reflected PUI temperature we investigate here are as follows.

1.
We estimate the reflected PUI temperature immediately downstream of the HTS using MHD pressure balance while still assuming that the electron temperature in the IHS is equal to the MHD temperature (T_e = T_MHD). The reflected PUI temperature immediately downstream of the HTS in this case is given by
$\begin{eqnarray}&&{T}_{{\rm{p}}}^{\mathrm{ref}}=\displaystyle \frac{{T}_{\mathrm{MHD}}-(1-\alpha ){T}_{{\rm{p}}}^{\mathrm{core}}-\alpha (1-\beta ){T}_{{\rm{p}}}^{\mathrm{tr}}}{\alpha \beta }.\end{eqnarray} \tag{ 14 }$
The total proton temperature immediately downstream of the HTS is ∼2 × 10⁶ K (=T_MHD), and the reflected PUI temperature using Equation (14) is similar to the previous case (∼9 × 10⁷ K). The temperature of the core SW, transmitted PUIs, and reflected PUIs at any point in the IHS is given by
$\begin{eqnarray}\begin{array}{rcl}{T}_{{\rm{p}}}^{\mathrm{core}}({\boldsymbol{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{core}}}{{T}_{\mathrm{MHD}}}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\boldsymbol{r}}),\\ {T}_{{\rm{p}}}^{\mathrm{tr}}({\boldsymbol{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{tr}}}{{T}_{\mathrm{MHD}}}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\boldsymbol{r}}),\\ {T}_{{\rm{p}}}^{\mathrm{ref}}({\boldsymbol{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{ref}}}{{T}_{\mathrm{MHD}}}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\boldsymbol{r}}).\end{array}\end{eqnarray} \tag{ 15 }$
The plots in the left panels of Figure 6 correspond to this case. The top panel shows the variation of the fractional difference in spectral slope between the data and simulation for a range of compression ratios. The best-fit shock compression ratio obtained in this case is the same as before (i.e., r = 2.4). The bottom panel shows the ENA spectra from both data and simulation using the best shock compression ratio. As the reflected PUI temperature here is similar to that used in Section 5.1, we do not see much difference in ENA fluxes from the reflected PUIs compared to the right panel of Figure 2.
2.
Next, we consider an electron temperature in the IHS that is different from the ion temperature (see also case IV of Heerikhuisen et al. 2019). More precisely, we assume that the electron temperature immediately downstream of the HTS is 0.1 times the MHD temperature and increases linearly up to 800 au along the plasma flow streamline, after which it becomes equal to the proton (and therefore MHD) temperature. The MHD temperature fraction for electrons immediately downstream of the HTS is obtained by assuming their temperature to be equal to the core SW temperature ( ${T}_{e}={T}_{{\rm{p}}}^{\mathrm{core}}\,\mathrm{and}\,{\gamma }_{{\rm{e}}}={T}_{{\rm{p}}}^{\mathrm{core}}/{T}_{\mathrm{MHD}}\sim 0.1$ ). Hence, the total proton temperature at any point downstream of the HTS is given by
$\begin{eqnarray}&&\,{T}_{{\rm{p}}}({\boldsymbol{r}})\,=\,{T}_{\mathrm{MHD}}({\boldsymbol{r}})\left[2-\left(\displaystyle \frac{\min (l,L)(1-{\gamma }_{{\rm{e}}})}{L}+{\gamma }_{{\rm{e}}}\right)\right],\,\end{eqnarray} \tag{ 16 }$
where l is the distance along a streamline from the HTS, L is the distance beyond which we assume electron and proton temperatures are equal, and γ_e is the fraction of the MHD temperature for electrons at the HTS. As mentioned before, we use L = 800 au and γ_e = 0.1. The total proton temperature immediately downstream of the HTS from Equation (16) is almost double the MHD temperature, T_p = T_MHD(2 − γ_e). The reflected PUI temperature immediately downstream of the HTS using MHD pressure balance is given by
$\begin{eqnarray}&&\,{T}_{{\rm{p}}}^{\mathrm{ref}}\,=\,\displaystyle \frac{{T}_{\mathrm{MHD}}(2-{\gamma }_{{\rm{e}}})-(1-\alpha ){T}_{{\rm{p}}}^{\mathrm{core}}-\alpha (1-\beta ){T}_{{\rm{p}}}^{\mathrm{tr}}}{\alpha \beta }.\end{eqnarray} \tag{ 17 }$
This expression gives a higher temperature for the reflected PUIs compared to case 1 (Equation (14)) because of the higher total proton temperature just downstream of the HTS. The temperature of the core SW, transmitted PUIs, and reflected PUIs at any point in the IHS is given by
$\begin{eqnarray}\begin{array}{rcl}{T}_{{\rm{p}}}^{\mathrm{core}}({\bf{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{core}}}{{T}_{\mathrm{MHD}}(2-{\gamma }_{{\rm{e}}})}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\bf{r}})\\ & & \times \left[2-\left(\displaystyle \frac{\min (l,L)(1-{\gamma }_{{\rm{e}}})}{L}+{\gamma }_{{\rm{e}}}\right)\right],\\ {T}_{{\rm{p}}}^{\mathrm{tr}}({\bf{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{tr}}}{{T}_{\mathrm{MHD}}(2-{\gamma }_{{\rm{e}}})}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\bf{r}})\\ & & \times \left[2-\left(\displaystyle \frac{\min (l,L)(1-{\gamma }_{{\rm{e}}})}{L}+{\gamma }_{{\rm{e}}}\right)\right],\end{array}\end{eqnarray}$

$\begin{eqnarray}\begin{array}{rcl}{T}_{{\rm{p}}}^{\mathrm{ref}}({\bf{r}}) & = & {\left(\displaystyle \frac{{T}_{{\rm{p}}}^{\mathrm{ref}}}{{T}_{\mathrm{MHD}}(2-{\gamma }_{{\rm{e}}})}\right)}_{\mathrm{HTS}}{T}_{\mathrm{MHD}}({\bf{r}})\\ & & \times \left[2-\left(\displaystyle \frac{\min (l,L)(1-{\gamma }_{{\rm{e}}})}{L}+{\gamma }_{{\rm{e}}}\right)\right].\end{array}\end{eqnarray} \tag{ 18 }$
The total proton temperature in this case (Equation (16)) changes along a streamline by a factor $\left[2-\left(\tfrac{\min (l,L)(1-{\gamma }_{{\rm{e}}})}{L}+{\gamma }_{{\rm{e}}}\right)\right]$ of the MHD temperature. As the total proton temperature changes along a streamline, the individual ion temperatures change to preserve the same ratio. The plots in the right panels of Figure 6 correspond to this case; the top and bottom panels follow the same explanations as in case 1. The higher temperature for the reflected PUIs, represented by a Maxwellian distribution, results in more particles having energy higher than the IBEX-Hi range. This reduces the number of PUIs in the IBEX-Hi energy range, hence ensuring a smaller number of ENAs. As there are fewer ENAs from the reflected PUIs at 4.29 keV compared to case 1, more ENAs from the transmitted PUIs are required to match the spectral slope. This is achieved by a higher value of the shock compression ratio (r = 2.6) and hence a higher transmitted PUI temperature in this case.

**Figure 6.** (Top panels) Variation of the fractional difference in spectral slopes between IBEX-Hi data and simulated IHS ENA fluxes for a range of shock compression ratios (1.3 ≤ r ≤ 3.5) in the V2 direction. (Bottom panels) ENA spectra in the V2 direction from both data and simulation using the best shock compression ratio corresponding to the minimum value of ∣Δγ_fr∣. The two columns correspond to two different estimates for the reflected PUI temperature downstream of the HTS: (left panels) MHD T_e and ${T}_{{\rm{p}}}^{\mathrm{ref}}$ from the MHD pressure balance and (right panels) variable T_e and ${T}_{{\rm{p}}}^{\mathrm{ref}}$ from the MHD pressure balance.
Download figure:
Standard image High-resolution image

**Figure 6.** (Top panels) Variation of the fractional difference in spectral slopes between IBEX-Hi data and simulated IHS ENA fluxes for a range of shock compression ratios (1.3 ≤ r ≤ 3.5) in the V2 direction. (Bottom panels) ENA spectra in the V2 direction from both data and simulation using the best shock compression ratio corresponding to the minimum value of ∣Δγ_fr∣. The two columns correspond to two different estimates for the reflected PUI temperature downstream of the HTS: (left panels) MHD T_e and ${T}_{{\rm{p}}}^{\mathrm{ref}}$ from the MHD pressure balance and (right panels) variable T_e and ${T}_{{\rm{p}}}^{\mathrm{ref}}$ from the MHD pressure balance.
Download figure:
Standard image High-resolution image

5.2. Compression Ratio in Other Directions

We now estimate the shock compression ratio in several directions of the sky using the method developed in Section 4 and validated from V2 observations (Section 5.1). As we are using IBEX-Hi data between 2009 and 2011 that correspond to the IHS plasma coming from the Sun during solar minimum, and the simulated background heliosphere corresponds to solar maximum, we limit our analysis close to the ecliptic plane (within ±30°). Note that there is a time delay of 2–6 yr (depending on the direction in the sky) between the SW measurement at 1 au and the ENA measurement at 1 au (see Zirnstein et al. 2017). To avoid the influence of the ribbon fluxes on this analysis, we choose directions sufficiently far away from the ribbon location. The electron temperature downstream of the HTS is assumed to be equal to the MHD temperature (T_e = T_MHD), and the temperature of the reflected PUIs just downstream of the HTS is estimated using Equation (8) (Z10 estimate). The upstream bulk flow speeds at different locations of the HTS are taken from the MHD kinetic simulation and then rescaled to the upstream bulk plasma flow speed observed along the V2 direction; i.e., the upstream flow speed at (θ, ϕ) in the sky is given by

$\begin{eqnarray}&&{u}_{1}(\theta ,\phi )={u}_{1,\mathrm{MHD}}(\theta ,\phi )\times {\left(\displaystyle \frac{{u}_{1,\mathrm{data}}}{{u}_{1,\mathrm{MHD}}}\right)}_{{\rm{v}}2}.\end{eqnarray} \tag{ 19 }$

The variation of ∣Δγ_fr∣ with the shock compression ratio for six different directions in the sky is shown in Figure 8. These six directions are represented by the violet Greek letters , ζ, ι, λ, ξ, and τ in Figure 7. The best compression ratio corresponding to the minimum value of ∣Δγ_fr∣ for each direction is tabulated in Table 1. The compression ratios in these directions are close to the compression ratio at V2, with a maximum difference of 0.3. A slight variation in the compression ratio in different directions is likely due to (1) different values of the upstream bulk plasma flow speed and (2) flow changes in the IHS in different directions. The upstream bulk flow speed is maximum in the nose and decreases toward the flanks of the heliosphere. The minimum value of the compression ratio, r = 2.1, is observed on the port side of the heliosphere at (27°, 15°). The low compression in this direction is likely due to the steeper IBEX spectrum on the flanks of the heliosphere, which requires a smaller compression ratio (lower PUI temperature) to compensate.

**Figure 7.** Footpoint locations of plasma flow streamlines at the HTS for plasma flowing along each radial step outward for six different LOSs in a Mollweide projection, similar to Figure 5.
Download figure:
Standard image High-resolution image

epsilon — **Figure 8.** Variation of ∣Δγ_fr∣ with shock compression ratio at six different directions in the sky. These directions are marked by the violet Greek letters , ζ, ι, λ, ξ, and τ in Figure 7.
Download figure:
Standard image High-resolution image

Table 1. HTS Compression Ratio

Direction	(Long., Lat.)	Compression	Flux	Distance to
		Ratio (r)	Factor	HTS (au)
V2	(291°, −33°)	2.4	3.3	75.8

Nose 1 ()	(255°, 21°)	2.2	4.7	74.0
Nose 2 (ζ)	(255°, 15°)	2.3	3.8	77.6
Port 1 (ι)	(333°, −15°)	2.2	4.6	84.3
Port 2 (λ)	(27°, 15°)	2.1	3.7	107.6
Starboard (ξ)	(117°, −9°)	2.2	3.5	113.5
Downwind (τ)	(69°, 3°)	2.3	3.0	117.2

Note. The flux factor for a given direction is calculated by dividing the ENA flux from IBEX-Hi at energy 1.74 keV by the total simulated ENA flux for the same energy. The radial distance to the HTS (from the Sun) in the front of the heliosphere (V1 and V2 directions) is less than 80 au. The smaller distance from the simulation compared to the observations is because we try to match the observed HP crossing distance in the V1 and V2 directions when simulating the background heliosphere.

Download table as: ASCII Typeset image

The ENA spectra in all six directions from both data and simulation corresponding to the best shock compression ratio are shown in Figure 9. Each plot in Figure 9 also contains the contribution to the total simulated ENA fluxes from each of the three populations of protons. The flux factors between the simulated ENA flux and data for all six directions are shown in Table 1. The simulated differential ENA fluxes along each radial step from the HTS to the HP at six different directions in the sky for each IBEX-Hi energy channels are shown in Figure 10. We can see that in the downwind tail (τ) direction, most of the ENA flux comes from within a few hundred au of the HTS, even though the tail is long. The footpoint locations of the plasma flow streamlines for each radial step outward from the HTS for all six LOSs are shown in Figure 7.

**Figure 10.** Differential ENA flux along each radial step from the HTS to the HP at six different directions in the sky for each IBEX-Hi energy channel. The cumulative normalized differential flux along each radial step in the same direction is also shown in the plot (dashed lines).
Download figure:
Standard image High-resolution image

6. Summary and Conclusions

We have presented a method to estimate the HTS compression ratio at several directions in the sky by a quantitative comparison of the ENA flux observed by IBEX-Hi with the simulated ENA flux for a range of compression ratios. This paper serves as a first step to fulfill one of the fundamental objectives of the IBEX mission: estimation of the global strength of the HTS. The simulations are based on a multi-Maxwellian description of protons in the IHS with a 3D steady-state background heliosphere simulated using a uniform slow SW at the inner boundary. The data used for the comparison are from IBEX-Hi over the 3 yr period between 2009 and 2011. This data set is chosen to ensure that (1) it is as close to the V2 HTS crossing time as possible, while (2) having a few years of IBEX data to average and improve the statistics. As IBEX data for this period correspond to the supersonic SW plasma during the solar minimum condition and our background heliosphere corresponds to solar maximum, we have limited our analysis to the lower-latitude region. Though the compression ratio presented here is only in directions close to the ecliptic, this method is equally applicable to other locations in the sky. The compression ratio estimation for higher latitudes requires the simulation of the background heliosphere with a latitude-dependent SW profile, i.e., a slow SW at low-to-mid-latitudes and a fast SW at high latitudes. The magnetic field structure over the polar region has a different geometry compared to the low-latitude region due to the nature of the Parker spiral. The magnetic field structure over the polar region will likely affect the accuracy of the shock compression estimated using this method because the PUI properties downstream of the HTS represented by Equations (2)–(5) may no longer be valid as the shock obliquity angle deviates away from 90°. A detailed analysis of this is beyond the scope of this paper and will be pursued in future.

The compression ratio estimated by this method is in good agreement with the large-scale compression ratio observed by V2. The shock compression ratios in other directions are close to that in the V2 direction. Unfortunately, we do not have a measurement of the plasma properties at the HTS in other directions of the sky against which to compare these results. A possible future data set for comparison will be provided by the New Horizons spacecraft if it continues to operate for a further 40–50 au and perhaps by a future interstellar probe mission. It is clear from V2 observations that the HTS evolves on short timescales (Richardson et al. 2008; Stone et al. 2008) yielding different local compression ratios, as the V2 spacecraft experienced multiple crossings before finally entering the IHS. The large-scale HTS compression ratio, however, should be more stable over time but may change with large-scale changes of the solar cycle. Of course, the timing of any crossing, solar maximum or minimum, will be important for any spacecraft crossing of the HTS, as with the possibility of interplanetary shocks interacting with the HTS (e.g., Story & Zank 1997; Donohue & Zank 1993). In that case, we have to use the ENA flux data for a later time from IBEX or the future IMAP mission for better comparison.

We acknowledge support from NASA IBEX grant 80NSSC20K0719, subaward SUB 0000395, and NASA grant NNX16AG83G. We also acknowledge the SW data at 1 au compiled by the OMNI database at the Space Physics Data Facility of NASA/GSFC (https://spdf.gsfc.nasa.gov/pub/data/omni/low_res_omni/). G.P.Z. and B.S. acknowledge the partial support of NSF EPSCoR RII-Track-1 Cooperative Agreement OIA-1655280, a NASA IMAP subaward under NASA contract 80GSFC19C0027, and NASA grant 18-DRIVE18_2-0029, Our Heliospheric Shield, 80NSSC20K0603.

Appendix: Reflected PUI Density from a Filled-shell Distribution

The PUI distribution in the supersonic SW far from the Sun is a filled shell, as represented by Equation (7) in the SW flow frame (Vasyliunas & Siscoe 1976; Z10). A sketch of the velocity phase-space portrait of this type of distribution is shown in Figure 11, where the green horizontal plane represents a specular reflection velocity (V_spec) in the shock rest frame (dotted black line) below which all particles are reflected; i.e., all particles with the x-component of velocity in the shock rest frame satisfying v_x < V_spec experience specular reflection. The sketch in the left panel is for the case of V_spec ≤ u₁, and the sketch in the right panel is for V_spec > u₁.

**Figure 11.** Sketch of the velocity phase-space portrait of a PUI filled-shell distribution in the SW frame just ahead of the HTS. The gray labels represent corresponding velocities in the shock rest frame (dotted black line). The sketch in the left panel is for the case of V_spec ≤ u₁, and the sketch in the right panel is for V_spec > u₁.
Download figure:
Standard image High-resolution image

The number density of reflected PUIs can be calculated by integrating the PUI filled-shell distribution in the SW frame over the shaded volume in Figure 11 as

$\begin{eqnarray}&&{n}_{{\rm{p}}}^{\mathrm{ref}}={\int }_{\mathrm{shadedvolume}}f(c)\,{d}^{3}c,\end{eqnarray} \tag{ A1 }$

i.e.,

$\begin{eqnarray}{n}_{{\rm{p}}}^{\mathrm{ref}}=\left\{\begin{array}{ll}\displaystyle \frac{3{n}_{\mathrm{PUI}}}{8\pi {u}_{1}^{3/2}}{\displaystyle \int }_{\phi =0}^{2\pi }{\displaystyle \int }_{\theta ={\theta }_{0}}^{\pi }{\displaystyle \int }_{c=\tfrac{| {v}_{{\rm{x}},{\rm{p}}}| }{\cos (\pi -\theta )}}^{{u}_{1}}{c}^{-3/2}\,{c}^{2}\,{dc}\,\sin \theta d\theta \,d\phi , & \mathrm{if}\ {V}_{\mathrm{spec}}\leqslant {u}_{1}\\ {n}_{\mathrm{PUI}}-\displaystyle \frac{3{n}_{\mathrm{PUI}}}{8\pi {u}_{1}^{3/2}}{\displaystyle \int }_{\phi =0}^{2\pi }{\displaystyle \int }_{\theta =0}^{{\theta }_{0}}{\displaystyle \int }_{c=\tfrac{{v}_{{\rm{x}},{\rm{p}}}}{\cos \theta }}^{{u}_{1}}{c}^{-3/2}\,{c}^{2}\,{dc}\,\sin \theta d\theta \,d\phi , & \mathrm{if}\ {V}_{\mathrm{spec}}\gt {u}_{1}\end{array}\right.,\end{eqnarray} \tag{ A2 }$

where v_x,p is the x-component of PUI velocity in the SW frame for the planar surface given by

$\begin{eqnarray}&&{v}_{{\text{}}x,p}={V}_{\mathrm{spec}}-{u}_{1}.\end{eqnarray} \tag{ A3 }$

Note that the number density of reflected PUIs when V_spec > u₁ (right panel of Figure 11) is obtained by subtracting the integral of the PUI distribution over the top white (smaller) part of the shell from the total number density of PUIs ( ${n}_{\mathrm{PUI}}\,={\int }_{0}^{2\pi }{\int }_{0}^{\pi }{\int }_{0}^{{u}_{1}}f(c){d}^{3}c$ ). We evaluate the reflected PUI number density in the two cases separately.

1.
Case I: V_spec ≤ u₁. On assuming $\theta ^{\prime} =\pi -\theta$ , we can write the number density of reflected PUIs for this case as
$\begin{eqnarray}&&{n}_{{\rm{p}}}^{\mathrm{ref}}=\displaystyle \frac{3{n}_{\mathrm{PUI}}}{8\pi {u}_{1}^{3/2}}{\int }_{\phi =0}^{2\pi }{\int }_{\theta ^{\prime} ={\theta }_{0}^{{\prime} }}^{0}{\int }_{c=\tfrac{| {v}_{{\rm{x}},{\rm{p}}}| }{\cos \theta ^{\prime} }}^{{u}_{1}}{c}^{1/2}\,{dc}\,\sin \theta ^{\prime} (-d\theta ^{\prime} )\,d\phi ,\end{eqnarray} \tag{ A4 }$
where ${\theta }_{0}^{{\prime} }=(\pi -{\theta }_{0})$ and can be obtained from the trigonometric relation $\cos {\theta }_{0}^{{\prime} }=\tfrac{| {v}_{{\text{}}x,p}| }{{u}_{1}}$ (see left panel of Figure 11). We find
$\begin{eqnarray*}\begin{array}{rcl}{n}_{{\rm{p}}}^{\mathrm{ref}} & = & \displaystyle \frac{3{n}_{\mathrm{PUI}}}{4{u}_{1}^{3/2}}{\displaystyle \int }_{0}^{{\theta }_{0}^{{\prime} }}{\left[\displaystyle \frac{{c}^{3/2}}{3/2}\right]}_{\tfrac{| {v}_{{\text{}}x,p}| }{{\cos }\theta ^{\prime} }}^{{u}_{1}}\,\sin \theta ^{\prime} d\theta ^{\prime} \\ & = & \displaystyle \frac{{n}_{\mathrm{PUI}}}{2}\left[{\displaystyle \int }_{0}^{{\theta }_{0}^{{\prime} }}\sin \theta ^{\prime} d\theta ^{\prime} -{\left(\displaystyle \frac{| {v}_{{\text{}}x,p}| }{{u}_{1}}\right)}^{3/2}{\displaystyle \int }_{0}^{{\theta }_{0}^{{\prime} }}\displaystyle \frac{\sin \theta ^{\prime} d\theta ^{\prime} }{{\cos }^{3/2}\theta ^{\prime} }\right]\\ & = & \displaystyle \frac{{n}_{{PUI}}}{2}\left[\left(1-\displaystyle \frac{| {v}_{{\text{}}x,p}| }{{u}_{1}}\right)-{\left(\displaystyle \frac{| {v}_{{\text{}}x,p}| }{{u}_{1}}\right)}^{3/2}\times 2\left(\displaystyle \frac{1}{\sqrt{\tfrac{| {v}_{{\text{}}x,p}| }{{u}_{1}}}}-1\right)\right].\end{array}\end{eqnarray*}$
Using ∣v_x,p∣ = ∣V_spec − u₁∣ = u₁ − V_spec, the above expression becomes
$\begin{eqnarray}&&{n}_{{\rm{p}}}^{\mathrm{ref}}={n}_{\mathrm{PUI}}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}-{\left(1-\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}\right)}^{3/2}\left(\displaystyle \frac{1}{\sqrt{1-\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}}}-1\right)\right].\end{eqnarray} \tag{ A5 }$
2.
Case II: V_spec > u₁. The maximum azimuthal angle θ₀ for the planar top surface in the right panel of Figure 11 is given by $\cos {\theta }_{0}=\tfrac{{v}_{{\text{}}x,p}}{{u}_{1}}$ . The ϕ, c, and θ integrals are the same as before, and, after evaluating these integrals, the number density for reflected PUIs becomes
$\begin{eqnarray*}&&{n}_{{\rm{p}}}^{\mathrm{ref}}={n}_{\mathrm{PUI}}-\displaystyle \frac{{n}_{\mathrm{PUI}}}{2}\left[\left(1-\displaystyle \frac{{v}_{{\text{}}x,p}}{{u}_{1}}\right)-{\left(\displaystyle \frac{{v}_{{\text{}}x,p}}{{u}_{1}}\right)}^{3/2}\times 2\left(\displaystyle \frac{1}{\sqrt{\tfrac{{v}_{{\text{}}x,p}}{{u}_{1}}}}-1\right)\right].\end{eqnarray*}$
Again, using v_x,p = V_spec − u₁ in the above equation, the number density of reflected PUIs becomes
$\begin{eqnarray}&&{n}_{{\rm{p}}}^{\mathrm{ref}}={n}_{{\rm{PUI}}}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}+{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}-1\right)}^{3/2}\left(\displaystyle \frac{1}{\sqrt{\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}-1}}-1\right)\right].\end{eqnarray} \tag{ A6 }$

Finally, combining Equations (A5) and (A6) yields the final expression for the number density of reflected PUIs,

$\begin{eqnarray}{n}_{{\rm{p}}}^{\mathrm{ref}}=\left\{\begin{array}{ll}{n}_{\mathrm{PUI}}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}-{\left(1-\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}\right)}^{3/2}\left(\displaystyle \frac{1}{\sqrt{1-\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}}}-1\right)\right], & \mathrm{if}\ {V}_{\mathrm{spec}}\leqslant {u}_{1}\\ {n}_{\mathrm{PUI}}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}+{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}-1\right)}^{3/2}\left(\displaystyle \frac{1}{\sqrt{\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}-1}}-1\right)\right], & \mathrm{if}\ {V}_{\mathrm{spec}}\gt {u}_{1}\end{array}\right..\end{eqnarray} \tag{ A7 }$

It is important to note that Equation (A7) takes the form of the expression given in Z10 (note that there is a typo in that paper) under the condition $\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}\lt \lt 1$ , as shown below:

$\begin{eqnarray}\begin{array}{rcl}{n}_{{\rm{p}}}^{\mathrm{ref}} & = & {n}_{\mathrm{PUI}}\left[\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}-\left(1-\displaystyle \frac{3}{2}\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}+\displaystyle \frac{3}{8}{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}\right)}^{2}+O(3)\right)\right.\\ & & \left.\times \left(\left(1+\displaystyle \frac{1}{2}\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}+\displaystyle \frac{3}{8}{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{{u}_{1}}\right)}^{2}+O(3)\right)-1\right)\right],\end{array}\end{eqnarray} \tag{ A8 }$

i.e.,

$\begin{eqnarray}&&{n}_{{\rm{p}}}^{\mathrm{ref}}={n}_{\mathrm{PUI}}\displaystyle \frac{3}{2}{\left(\displaystyle \frac{{V}_{\mathrm{spec}}}{2{u}_{1}}\right)}^{2}.\end{eqnarray} \tag{ A9 }$

This expression is valid for any value of $\tfrac{{V}_{\mathrm{spec}}}{{u}_{1}}\lt 0.5$ within 5% of the general expression (Equation (A7)).

Strength of the Termination Shock Inferred from the Globally Distributed Energetic Neutral Atom Flux from IBEX

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Simulation Method