Constraints on the IBEX Ribbon's Origin from Its Evolution over a Solar Cycle

As with our previous studies on modeling the IBEX ribbon (Zirnstein et al. 2013, 2015, 2018, 2019a, 2019b, 2020a, 2021a), we utilize results from a three-dimensional, global magnetohydrodynamic (MHD) simulation of the SW-LISM interaction based on MS-FLUKSS (Pogorelov et al. 2008; Heerikhuisen et al. 2015), where we simulate secondary ENA fluxes from outside the heliopause using different assumptions for PUI pitch angle scattering: weak pitch angle scattering (Zirnstein et al. 2018) or strong pitch angle scattering using the spatial retention model (Schwadron & McComas 2013; Zirnstein et al. 2019a). The heliosphere simulation solves ideal MHD equations for a single-fluid plasma and the Boltzmann equation for propagation of neutral H atoms. The plasma and neutrals are coupled by charge-exchange source terms (Heerikhuisen et al. 2015).

The SW boundary conditions at 1 au for the steady-state heliosphere simulation are chosen to be the same as those used by Zirnstein et al. (2021b). At low latitudes (<∣ ±37°∣, from the OMNI database in 2004–2009), the SW speed is 449 km s⁻¹, plasma density is 6.53 cm⁻³, temperature is 1.02 × 10⁵ K, and radial magnetic field component is 37.4 μG. At high latitudes (>∣ ±37°∣, from Ulysses' third polar scan (McComas et al. 2008)), SW speed is 743 km s⁻¹, density is 2.23 cm⁻³, temperature is 2.98 × 10⁵ K, and radial magnetic field component is 34.7 μG.

The LISM boundary conditions are applied at 2500 au from the Sun, which we note is larger than the 1000 au outer boundary we have used in previous studies. The radial grid resolution in the ribbon source region ranges from Δr = 0.9 to 1.3 au for r = 120–200 au. The angular resolution is (Δθ, Δφ) = (15, 20), which is sufficiently small to capture the structure of the modeled ribbon (whose width is >10°). The LISM boundary conditions are chosen based on IBEX, Voyager, and New Horizons (NH) SWAP observations. The pristine magnetic field vector is pointed toward (22728, 3462) in ecliptic J2000 with magnitude 2.93 μG (Zirnstein et al. 2016b). The interstellar plasma number density is 0.09 cm⁻³, and the neutral H number density is 0.17 cm⁻³. Our steady-state model overestimates the thickness of the heliosheath, such that we cannot simultaneously reproduce both the heliospheric termination shock (HTS) and heliopause distances observed by the Voyagers. Therefore, this combination of LISM number densities is designed such that the middle of the simulated heliosheath, averaged in the Voyager 1 and 2 directions, approximately agrees with the average of Voyager 1 and 2's observations of the middle of the heliosheath. Finally, similarly to Zirnstein et al. (2021b), we assume that 4% of the plasma density in the SW consists of alpha particles and 10% of the plasma density in the LISM consists of He⁺ ions, such that those percentages of the ions are not included in the H–H⁺ charge-exchange source terms of the MHD-kinetic model of the heliosphere. This ad hoc modification of the plasma density in the charge exchange source terms yields a filtrated interstellar neutral H density at the upwind HTS of ∼0.13 cm⁻³, in agreement with extrapolations from NH SWAP observations (Swaczyna et al. 2020).

The neutral particles that are the parent population for PUIs outside the heliopause originate from (1) protons in the supersonic SW, i.e., the core SW protons, (2) interstellar H⁺ PUIs embedded in the outward-propagating SW, and (3) hot ions in the IHS that were accelerated at the HTS. The neutralized SW source of the ribbon has been included in all past secondary ENA models (Chalov et al. 2010; Heerikhuisen et al. 2010; Schwadron & McComas 2013; Zirnstein et al. 2013). Neutralized hot heliosheath ions have been included in several studies of the "partial shell" model of the ribbon using kinetically driven neutral source terms (Heerikhuisen et al. 2010; Zirnstein et al. 2013, 2016a), but these studies suggest heliosheath neutrals do not contribute significantly to the ENA fluxes measured by IBEX-Hi below ∼6 keV, and therefore we do not include them. The inclusion of neutralized interstellar H⁺ PUIs in our current study, however, is motivated by the analysis of Schwadron & McComas (2019), who showed that these particles may produce the majority of ENAs in the IBEX-Hi ESA 4–6 passbands (∼2–6 keV). As we will show in our study, the neutralized core SW component (hereafter called NSW) generally dominates the ribbon signal at energies ≲2 keV, and the neutralized interstellar H⁺ PUIs in the SW (hereafter called NPI) dominates the source of the ribbon signal at energies >2 keV. The latitudinal distribution of their contributions to the ribbon are also distinctly different, due to their vastly different thermal energies in the SW frame.

One main difference in this study compared to our previous work is that we calculate NSW and NPI fluxes for individual years (time t) based on the IPS-based empirical model of Sokół et al. (2020). Because IPS-based SW speeds generally overestimate near-ecliptic in situ measurements (Tokumaru et al. 2021), the SW speed is subtracted by a constant value uniformly at all latitudes to match OMNI at low latitudes. However, it is unclear whether IPS-derived speeds should be shifted uniformly in speed at all latitudes, or differently in the northern and southern hemispheres, or not at all except near the ecliptic plane. Therefore, we also present results where no shifting is performed in order to understand the effects of this assumption on our ribbon model results. The SW density is determined by taking measurements from the OMNI database over a similar time period, namely the central time from the IPS model ±0.375 yr and calculating the density as a function of latitude, assuming the dynamic pressure is constant with latitude. We also assume the SW temperature is constant over latitude and time (see Section 2.2.1, Equation (4)). While this is not strictly true, the effects of SW temperature on the model results are small compared to the radial spread of the distribution. The PUI temperature is modeled more accurately using the filled-shell distribution as a function of distance from the Sun (see also Sections 2.2.2 and 2.2.3).

2.2.1. Neutral SW H⁺ Flux

Building upon our previous work (Swaczyna et al. 2016; Zirnstein et al. 2019a, 2019b), we model the NSW source term by integrating the flux of neutralized SW H⁺ from 1 au to the HTS for SW conditions derived from IPS observations by Sokół et al. (2020). We calculate the NSW flux for each annual time period, denoted by time t_j , where j is the index corresponding to years shown in Table A1 of Sokół et al. (2020):

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{NSW}}\left({r}_{\mathrm{HTS}},\theta ,\phi ,v,{t}_{j}\right)\\ \ =\,{\displaystyle \int }_{{r}_{0}}^{{r}_{\mathrm{HTS}}}\left[{n}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right){u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right)\right]{\left(\displaystyle \frac{{r}_{0}}{{r}_{\mathrm{HTS}}}\right)}^{2}\\ \ \ \ \times \,{{\rm{e}}}^{-\left\langle {n}_{{\rm{H}}}\left({r}_{0}:r,\theta ,\varphi \right)\right\rangle {\sigma }_{\mathrm{ex}}\left(v\right)r}\left[{n}_{{\rm{H}}}\left(r,\theta ,\phi \right){\sigma }_{\mathrm{ex}}\left(v\right)\right]\\ \ \ \ \times \,{N}_{\mathrm{MB}}\left(r,\theta ,\phi ,v,{t}_{j}\right)\displaystyle \frac{1}{{m}_{{\rm{p}}}v}{dr},\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray*}&&\begin{array}{l}{N}_{\mathrm{MB}}\left(r,\theta ,\phi ,v,{t}_{j}\right)=\displaystyle \frac{1}{\sqrt{2\pi }\delta v\left(r,\theta ,{t}_{j}\right)}\exp \\ \times \,\left[\displaystyle \frac{-{\left(v-{u}_{\mathrm{SW}}\left(r,\theta ,\phi ,{t}_{j}\right)\right)}^{2}}{2{\left(\delta v\left(r,\theta ,{t}_{j}\right)\right)}^{2}}\right],\\ {u}_{\mathrm{SW}}\left(r,\theta ,\phi ,{t}_{j}\right)={u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right)\\ \times \,\left[1-\left(1-\displaystyle \frac{1}{2}\displaystyle \frac{\gamma -1}{2\gamma -1}\right)\displaystyle \frac{r}{{\lambda }_{\mathrm{ml}}\left(r,\theta ,\phi ,{t}_{j}\right)}\right],\\ {\lambda }_{\mathrm{ml}}\left(r,\theta ,\phi ,{t}_{j}\right)=\left[{\sigma }_{{ex}}\left({u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right)\right)\left\langle {n}_{H}\left({r}_{0}:r,\theta ,\phi \right)\right\rangle \right.\\ +\,{\left({n}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right){u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right)\right)}^{-1}\left({\nu }_{{\rm{H}}}\left({r}_{0}\right)\right.\\ \times \,{\left.\left.\left\langle {n}_{{\rm{H}}}\left({r}_{0}:r,\theta ,\phi \right)\right\rangle +4{\nu }_{\mathrm{He}}\left({r}_{0}\right){n}_{\mathrm{He}}\right)\right]}^{-1},\end{array}\end{eqnarray*}$$

where N_MB is the Maxwell–Boltzmann speed distribution, u_SW is the SW speed slowing with distance from the Sun, and λ_ml is the mean free path length for neutralization of SW protons. The term $\left\langle {n}_{{\rm{H}}}\left({r}_{0}:r,\theta ,\phi \right)\right\rangle$ in the exponential function is the arithmetic mean of neutral H density from r₀ to r in direction $\left(\theta ,\phi \right)$. Because the rates of photoionization for H and He are small compared to H–H⁺ charge exchange, we assume they are constant in time with values at r₀ = 1 au of ν_H = 1.3 × 10⁻⁷ s⁻¹ and ν_He = 1 × 10⁻⁷ s⁻¹ (Sokół et al. 2020). We integrate the SW flux (first term), as it adiabatically expands to the HTS and attenuates due to charge exchange losses (exponential term), multiplied by the probability of charge exchange with neutral H, and multiplied by a one-dimensional Maxwell–Boltzmann speed distribution to account for variations in SW speed that occur over approximately ∼1 yr intervals (i.e., the cadence of the IPS-derived SW speed model). To do this, first we find the standard deviation of SW speed at 1 au from the OMNI database corresponding to the IPS measurement times (which is typically δ v_u ∼ 100 km s⁻¹ at 1 au) and use New Horizons SWAP observations of SW speeds between 20 and 40 au from the Sun to determine how the speed deviation δ v_u changes with distance. We find from SWAP observations that the SW speed spread drops significantly at large distances (becoming δ v_u ∼ 25 km s⁻¹ halfway to the HTS). We find a best-fit power-law relation between OMNI and SWAP observations of the SW speed standard deviation (testing bin sizes between 0.5 and 2 yr, which asymptotically yields a power-law slope of −0.37 for 2 yr bins), such that

$$\begin{eqnarray}&&\delta {v}_{u}\left(r,\theta ,{t}_{j}\right)=\delta {v}_{u,\mathrm{OMNI}}\left({t}_{j}\right)\left[\displaystyle \frac{{u}_{\mathrm{IPS}}\left(\theta ,{t}_{j}\right)}{{u}_{\mathrm{OMNI}}\left({t}_{j}\right)}\right]{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-0.37}.\end{eqnarray} \tag{ 2 }$$

Because we do not know the variation of SW speeds at high latitudes except from Ulysses observations before 2009, we make a simplified assumption that the SW speed standard deviation scales linearly with SW speed. This is motivated by the fact that, as the Sun rotates over time, the mid-latitude SW streams create stream interactions and thus larger speed variations. This assumption does break down at high latitudes during solar minimum when the SW speed appears steady (McComas et al. 2008), likely affecting the northern section of the ribbon model above +40° latitude, i.e., the ribbon "knot."

We also include a term for the thermal motion of the SW,

$$\begin{eqnarray}&&\delta {v}_{T}\left(r,\theta ,{t}_{j}\right)=\delta {v}_{T,\mathrm{OMNI}}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-0.67},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\delta {v}_{T,\mathrm{OMNI}}=\sqrt{{k}_{{\rm{B}}}\left(9\times {10}^{4}\ {\rm{K}}\right)/{m}_{p}},\end{eqnarray} \tag{ 4 }$$

where we have assumed a constant SW proton temperature of 9 × 10⁴ K at 1 au, for simplicity. The exponent in Equation (3) is also derived by fitting a power-law relation between OMNI and 1 yr binned SWAP observations of SW proton temperature. This yields the total speed deviation,

$$\begin{eqnarray}&&\delta v\left(r,\theta ,{t}_{j}\right)=\sqrt{{\left[\delta {v}_{u}\left(r,\theta ,{t}_{j}\right)\right]}^{2}+{\left[\delta {v}_{T}\left(r,\theta ,{t}_{j}\right)\right]}^{2}}.\end{eqnarray} \tag{ 5 }$$

We note that the second term in Equation (5) is small compared to the first term and can be neglected.

The NSW flux decreases with distance beyond the HTS due to adiabatic expansion and due to losses by charge exchange outside the heliopause (where plasma densities are much higher than inside the heliopause), yielding (Zirnstein et al. 2021a)

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{NSW}}\left(r,\theta ,\phi ,v,{t}_{j}\right)={I}_{\mathrm{NSW}}\left({r}_{\mathrm{HTS}},\theta ,\phi ,v,{t}_{j}\right)\\ \times {\left(\displaystyle \frac{{r}_{\mathrm{HTS}}\left(\theta ,\phi \right)}{r}\right)}^{2}{{\rm{e}}}^{-{\displaystyle \int }_{{r}_{\mathrm{HP}}}^{r}{n}_{{\rm{p}}}\left(r,\theta ,\phi \right){\sigma }_{\mathrm{ex}}{{dr}}^{{\prime} }}.\end{array}\end{eqnarray} \tag{ 6 }$$

2.2.2. Neutral Interstellar PUI H⁺ Flux

In contrast to SW protons, which can be approximately described with a Maxwell–Boltzmann speed distribution, interstellar PUIs are characterized by a nonthermalized, filled-shell distribution, due to the process in which they are injected into the SW and scatter in pitch angle (Vasyliunas & Siscoe 1976), in addition to experiencing nonadiabatic cooling with distance from the Sun (McComas et al. 2021). The NPI flux source term is therefore given as

$$\begin{eqnarray}&&\begin{array}{l}{I}_{\mathrm{NPI}}\left({r}_{\mathrm{HTS}},\theta ,\phi ,v,{t}_{j}\right)\\ ={\displaystyle \int }_{{r}_{0}}^{{r}_{\mathrm{HTS}}}\left[{n}_{\mathrm{PI}}\left(r,\theta ,\phi ,{t}_{j}\right){u}_{\mathrm{SW}}\left(r,\theta ,\phi ,{t}_{j}\right)\right]{\left(\displaystyle \frac{r}{{r}_{\mathrm{HTS}}}\right)}^{2}\\ \times \left[{n}_{{\rm{H}}}\left({{rw}}^{\alpha },\theta ,\phi \right){\sigma }_{\mathrm{ex}}\left(v\right)\right]{N}_{\mathrm{FS}}\left(r,\theta ,\phi ,v,{t}_{j},\alpha \right)\displaystyle \frac{1}{{m}_{p}v}{dr},\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray*}&&\begin{array}{l}{N}_{\mathrm{FS}}\left(r,\theta ,\phi ,v,{t}_{j},\alpha \right)=\displaystyle \frac{\alpha }{2\left(\alpha -1\right){V}_{\mathrm{inj}}}\left[1-{w}^{\alpha -1}\right]\ \ \mathrm{for}\ w\lt 1,\\ {V}_{\mathrm{inj}}=\left|{u}_{\mathrm{SW}}\left(r,\theta ,\phi ,{t}_{j}\right)-{u}_{{\rm{H}}}\cos \beta \right|,\\ w=\left|v-{u}_{\mathrm{SW}}\left(r,\theta ,\phi ,{t}_{j}\right)\right|/{V}_{\mathrm{inj}}.\end{array}\end{eqnarray*}$$

The function N_FS is based on the generalized filled-shell distribution (Chen et al. 2013; Swaczyna et al. 2020; McComas et al. 2021) integrated over the speed components perpendicular to the radial direction. The PUI cooling index is α ≅ 2.1 halfway to the HTS (McComas et al. 2021), therefore we assume α = 2.1 is constant with distance r in Equation (7). β is the angle between the SW radial vector and the LISM H inflow vector assumed to be moving at u_H = 22 km s⁻¹ (Lallement et al. 2010). The neutral H density term in Equation (7) is evaluated at distance rw^α , to compensate for the ionization cavity effect for neutral H close to the Sun. The interstellar PUI density ${n}_{\mathrm{PI}}\left(r,\theta ,\phi ,{t}_{j}\right)$ is calculated from first-order solutions to fluid transport equations from, e.g., Lee et al. (2009):

$$\begin{eqnarray}&&\begin{array}{l}{n}_{\mathrm{PI}}\left(r,\theta ,\phi ,{t}_{j}\right)=\displaystyle \frac{{r}_{0}^{2}}{{u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right)r}\left\langle {n}_{{\rm{H}}}\left({r}_{0}:r,\theta ,\phi \right)\right\rangle \\ \ \times \left[{\nu }_{{\rm{H}},0}+{n}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right){u}_{\mathrm{SW}}\left({r}_{0},\theta ,{t}_{j}\right){\sigma }_{\mathrm{ex}}\left({u}_{\mathrm{SW}}\right)\right],\end{array}\end{eqnarray} \tag{ 8 }$$

Similarly to the NSW flux (Section 2.2.1, Equation (6)), the NPI flux decreases in intensity beyond the HTS, due to adiabatic expansion and charge-exchange losses.

Examples of the NSW and NPI flux at the upwind HTS are shown in Figure 1. In general, the NSW flux has a narrower, more intense distribution centered near the SW speed at a given latitude (blue curves). The NPI fluxes are lower in intensity near the bulk SW speed but are more thermally broad, extending to lower and higher energies. The NPI flux does not have a sharp cutoff like the PUI distribution (McComas et al. 2021), because the radially integrated NPI flux includes PUIs created at different distances from the Sun, which were scattered onto shells with different cutoff speeds, due to the slowing of the SW with distance from the Sun. The NPI flux also has a minimum at the bulk SW speed (e.g., 0.3 keV in Figure 1(a)), corresponding to PUIs with very low speeds in the SW frame—which, however, could not be created, due to significant losses of neutral H inside the ionization cavity close to the Sun (Sokół et al. 2019).

Zoom In Zoom Out Reset image size

We show examples of NSW and NPI fluxes in 2009 and 2014, which reflect solar minimum and maximum conditions, respectively, in solar cycle 24. At low latitudes (Figures 1(a) and (d)), the neutral fluxes are similar at both times in the cycle. At 30° latitude, which is approximately the latitude at which the polar coronal holes (PCH) reach their lowest point during solar minimum, the NSW and NPI fluxes shift to higher energies in 2009. At 60° latitude, well within the northern PCH and at a latitude similar to the ribbon knot seen at 2.73 keV, the NSW and NPI fluxes reach their highest energies in 2009, due to fast SW from the PCH. Moreover, the NSW distribution is more flat-topped, which appears to be a result of a larger spread of SW speeds as defined in Equation (2), which increases with higher SW speeds. As we stated earlier, it is not clear if this assumption is correct at high latitudes, but we believe it is reasonable at latitudes <∣ ± 60°∣ near solar minimum and all latitudes near solar maximum.

We note that, in this study, we only consider neutral sources inside the HTS, which have been shown to be an essential component of the ribbon source. However, neutralized PUIs from the IHS will also contribute to the ribbon source signal, although we expect this to occur at low intensities from ∼0.5 to 3 keV. Above ∼3 keV, hot heliosheath neutrals may be a significant source for the ribbon signal (Zirnstein et al. 2016a; Schwadron & McComas 2019). This component will be analyzed and included in a future study.

2.2.3. Neutral SW and PUI Source Transverse Components

As we discuss the NSW and NPI distributions beyond the HTS, first it is important to note that, from hereafter, we refer to the local time variable as t, such that $t\left({\boldsymbol{r}}\right)={t}_{j}\left({{\boldsymbol{r}}}_{0}\right)-l\left({\boldsymbol{r}}\right)/v$, which is the time in the past when the ENA traveling at speed v was created at LOS length distance l from IBEX.

In this study, we model the ribbon in the limits of weak pitch angle scattering (Zirnstein et al. 2018) and strong pitch angle scattering (Schwadron & McComas 2013; Zirnstein et al. 2019a). In the limit of weak pitch angle scattering, we assume the NSW and NPI distributions outside the heliopause have thermal spreads transverse to the radial direction, with different temperatures. The transverse component of Equations (1) and (7) are given as (Zirnstein et al. 2021a)

$$\begin{eqnarray}&&\begin{array}{c}{N}_{\perp ,i}\left(r,\theta ,\phi ,v,t\right)=\frac{1}{\pi \delta {v}_{\perp ,i}^{2}\left(r,\theta ,\phi ,t\right)}\\ \ \times \exp \left[-\frac{{v}_{\perp }^{2}}{\delta {v}_{\perp ,i}^{2}\left(r,\theta ,\phi ,t\right)}\right]H\left(90^\circ -{\cos }^{-1}\left({\hat{{\boldsymbol{v}}}}_{p}\cdot \hat{{\boldsymbol{r}}}\right)\right),\end{array}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray*}&&{v}_{\perp }=v\sin \left({\cos }^{-1}\left({\widehat{{\boldsymbol{v}}}}_{p}\cdot \widehat{{\boldsymbol{r}}}\right)\right),\end{eqnarray*}$$

$$\begin{eqnarray*}&&\delta {v}_{\perp ,i}\left(r,\theta ,\phi ,t\right)=\sqrt{\displaystyle \frac{2{k}_{B}{T}_{\perp ,i,\mathrm{HTS}}}{{m}_{p}}{\left(\displaystyle \frac{{r}_{\mathrm{HTS}}\left(\theta ,\phi \right)}{r}\right)}^{2}},\end{eqnarray*}$$

where H is the Heaviside step function, and ${{\boldsymbol{v}}}_{p}={{\boldsymbol{v}}}_{p}\left(\omega ,\varphi \right)$ is the PUI velocity vector with the same magnitude as the ENA $\left(v\right)$ with pitch angle φ and gyro-phase ω. The transverse temperature at the HTS, T_⊥,i,HTS, for neutral component i (NSW or NPI) is

$$\begin{eqnarray}&&{T}_{\perp ,\mathrm{NSW},\mathrm{HTS}}=5\times {10}^{3}\ {\rm{K}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\begin{array}{l}{T}_{\perp ,\mathrm{NPI},\mathrm{HTS}}={\displaystyle \int }_{{r}_{0}}^{{r}_{\mathrm{HTS}}}{T}_{\perp ,\mathrm{PI}}\left(r,\theta ,\phi ,t\right)\\ \quad \times \,{\left(\displaystyle \frac{r}{{r}_{\mathrm{HTS}}\left(\theta ,\phi \right)}\right)}^{2}{W}_{\mathrm{NPI}}\left(r,\theta ,\phi ,t,\alpha \right){dr},\\ {T}_{\perp ,\mathrm{PI}}\left(r,\theta ,\phi ,t\right)=\displaystyle \frac{{m}_{p}}{3{k}_{B}}{u}_{\mathrm{SW}}^{2}\left(r,\theta ,\phi ,t\right)\\ \quad \times \,\displaystyle \frac{{E}_{2+2/\alpha }\left({\rm{\Lambda }}\right)}{\exp \left(-{\rm{\Lambda }}\right)-{\rm{\Lambda }}{E}_{1}\left({\rm{\Lambda }}\right)},\\ {W}_{\mathrm{NPI}}\left(r,\theta ,\phi ,t,\alpha \right)\\ =\,\displaystyle \frac{{n}_{\mathrm{PI}}\left(r,\theta ,\phi ,t\right){\left(\tfrac{r}{{r}_{{HTS}}\left(\theta ,\phi \right)}\right)}^{2}{n}_{{\rm{H}}}\left({{rw}}^{\alpha },\theta ,\phi \right)}{{\displaystyle \int }_{{r}_{0}}^{{r}_{\mathrm{HTS}}}{W}_{\mathrm{NPI}}\left({r}^{{\prime} }\right){{dr}}^{{\prime} }}.\end{array}\end{eqnarray} \tag{ 11 }$$

For simplicity, we assume the NSW temperature at the HTS is constant, a justifiable assumption due its small thermal spread (Florinski & Heerikhuisen 2017) compared to the width of the IBEX ribbon. The NPI temperature is much larger, however, and therefore we approximate its temperature using a weighted average of the interstellar PUI temperature from r₀ to r_HTS, weighted by the local PUI and neutral H densities. The interstellar PUIs, once neutralized, will have a temperature that decreases due to adiabatic expansion, hence the r² term in the first part of Equation (11). The interstellar PUI temperature T_⊥,PI is taken from Zirnstein et al. (2022). The weighting term W_NPI accounts for the number density of PUIs that, once neutralized by charge exchange with the local neutral H density, will also adiabatically expand with distance to the HTS.

The NPI distribution is much hotter than the NSW, and because of our approach to modeling the ribbon, our methodology does not account for the fact that we might produce PUIs whose initial pitch angle vectors do not originate from within the HTS. Therefore, we add one final criterion to the NPI source, by excluding any PUI whose initial velocity vector angle from the radial direction is larger than ${\tan }^{-1}\left(90\ \mathrm{au}/r\right)$, where 90 au roughly approximates the average HTS distance in the nose-ward hemisphere (where most of the ribbon appears in the sky). We note, however, that this exclusion yields a negligible change in our results.

The time-dependent, differential ENA flux at 1 au is calculated as

$$\begin{eqnarray}&&\begin{array}{c}J\left(\theta ,\phi ,v,{t}_{0}\right)={\int }_{{r}_{\mathrm{HP}}}^{{r}_{\mathrm{OB}}}\frac{1}{{m}_{p}}{f}_{0}\left(r,\theta ,\phi ,v,t\right){v}^{2}\\ \ \times {n}_{{\rm{H}}}\left(r,\theta ,\phi \right){\sigma }_{\mathrm{ex}}P\left(r,\theta ,\phi ,v,t\right){dr},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${f}_{0}\left(r,\theta ,\phi ,v,t\right)$ is the PUI distribution function along the IBEX line of sight (LOS) vector direction (θ, ϕ) for the weak or strong scattering models of the ribbon, following the work of Zirnstein et al. (2021a) (see their Equations (4), (6), and (9)), and $P\left(r,\theta ,\phi ,v,t\right)$ is the probability of ENA survival from their point of creation to the HTS (it should be noted that we compare to IBEX data corrected for ENA losses from 100 au, i.e., near the HTS, to 1 au). The integration is performed along the IBEX LOS from the heliopause (r_HP) to the outer boundary set at 700 au from the Sun (r_OB). As stated in Section 2.2.3, the time variable is t = t₀ − l/v, where l is the distance traveled by the ENA with speed v backward in time along the IBEX LOS. The variable t₀ is chosen as a function of longitude and latitude to coincide with IBEX's measurement times as it takes observations over the course of a full year in the "ram" frame.

2.3.1. Weak Scattering Model ENA Source

Because this study simulates time-dependent ENA flux, for completeness, we rewrite the equations for f₀ from Zirnstein et al. (2021a). The time-dependent PUI distribution in the weak scattering model is the same for both the NSW and NPI sources:

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0,i}(r,\theta ,\phi ,v,t)={\displaystyle \int }_{0}^{{t}_{\max }}W({t}^{{\prime} })\left[\displaystyle \frac{1}{{\rm{\Delta }}\varphi }{\displaystyle \int }_{-{\rm{\Delta }}\varphi /2}^{+{\rm{\Delta }}\varphi /2}\right.\\ \ \left.\times \,\left(\displaystyle \frac{1}{2\pi }{\displaystyle \int }_{0}^{2\pi }{S}_{i}({{\boldsymbol{r}}}^{{\prime} },\omega ,\varphi ,v,{t}^{{\prime} }){\tau }_{\mathrm{ex}}d\omega \right)d\varphi \right]{{dt}}^{{\prime} },\\ W\left({t}^{{\prime} }\right)=\displaystyle \frac{\exp \left(-{t}^{{\prime} }/{\tau }_{\mathrm{ex}}\right)}{{\displaystyle \int }_{0}^{{t}_{\max }}\exp \left(-{t}^{{\prime\prime} }/{\tau }_{\mathrm{ex}}\right){{dt}}^{{\prime\prime} }},\end{array}\end{eqnarray} \tag{ 13 }$$

where ${\tau }_{\mathrm{ex}}=1/\left({n}_{{\rm{H}}}{\sigma }_{\mathrm{ex}}{v}_{\mathrm{rel}}\right)$ is the 1/e charge-exchange timescale for neutralizing protons at relative interaction speed v_rel. The neutral source function is

$$\begin{eqnarray}&&{S}_{i}\left({{\boldsymbol{r}}}^{{\prime} },\omega ,\varphi ,v,{t}^{{\prime} }\right)={I}_{i}\left({{\boldsymbol{r}}}^{{\prime} },v,{t}^{{\prime} }\right){m}_{{\rm{p}}}{n}_{{\rm{p}}}\left({{\boldsymbol{r}}}^{{\prime} }\right){\sigma }_{\mathrm{ex}}{{vN}}_{\perp ,i}\left({{\boldsymbol{r}}}^{{\prime} },v,{t}^{{\prime} }\right).\end{eqnarray} \tag{ 14 }$$

Similarly to our previous work (Zirnstein et al. 2020b), we assume the initial PUI ring beam scatters onto a narrow partial shell with width Δφ = 5° centered on the initial pitch angle φ, and we integrate over gyro-phase ω assuming a uniform distribution over the partial shell. Our assumption of Δφ = 5°, as we tested in Zirnstein et al. (2020b), does not significantly affect our ribbon intensities compared to assuming 2° or 10°. We expect even less effect on the ribbon geometry because of the symmetry of Δφ around the initial pitch angle. The position vector ${{\boldsymbol{r}}}^{{\prime} }={r}^{{\prime} }\left({\theta }^{{\prime} },{\phi }^{{\prime} }\right)$ signifies the position of the PUI traced backward in time, using guiding center motion, created at time ${t}^{{\prime} }=t-{\int }_{0}^{{t}_{g}}{{dt}}^{{\prime\prime} }$, where t_g is the elapsed time for the guiding center motion of the PUI. The main difference in Equation (13) compared to Zirnstein et al. (2021a) is that the integration backward in time solves for a time-dependent source for the NSW or NPI flux, S_i , rather than assuming it is constant over time. The total PUI distribution is given as ${f}_{0}\left(r,\theta ,\phi ,v,t\right)={f}_{0,\mathrm{NSW}}\left(r,\theta ,\phi ,v,t\right)+{f}_{0,\mathrm{NPI}}(r,\theta ,\phi ,v,t)$.

2.3.2. Strong Scattering Model ENA Source

We simulate the ribbon in the strong scattering limit using the spatial retention mechanism (Schwadron & McComas 2013), separately for the NSW and NPI sources. For the case where PUIs are created from the NSW, where the effects of the small thermal spread of the NSW can be neglected (Schwadron & McComas 2013), the time-dependent PUI distribution inside the retention region is

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0,\mathrm{NSW}}\left(r,\theta ,\phi ,v,t\right)={\displaystyle \int }_{0}^{{t}_{\max }}W({t}^{{\prime} })[{S}_{\mathrm{NSW}}\left(r,\theta ,\phi ,v,{t}^{{\prime} }\right){\tau }_{\mathrm{ex}}\\ +\,f\left(z\leqslant {z}_{A}\right)]{{dt}}^{{\prime} },\mathrm{for}\ {v}_{| | }\leqslant {v}_{A},\end{array}\end{eqnarray} \tag{ 15 }$$

and outside the retention region it is

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0,\mathrm{NSW}}\left(r,\theta ,\phi ,v,t\right)={\displaystyle \int }_{0}^{{t}_{\max }}W({t}^{{\prime} })\\ \ \times \left[{S}_{\mathrm{NSW}}\left(r,\theta ,\phi ,v,{t}^{{\prime} }\right)\left[{\tau }_{\mathrm{ex}}-{\tau }_{1}\right]\right]{{dt}}^{{\prime} },\\ \quad \mathrm{for}\,{v}_{| | }\gt {v}_{A},\end{array}\end{eqnarray} \tag{ 16 }$$

where $z=r\ \tan \left(\pi /2-{\cos }^{-1}\left(\left|{{\boldsymbol{v}}}_{\mathrm{NSW}}\cdot {\boldsymbol{B}}\right|/\left|{{\boldsymbol{v}}}_{\mathrm{NSW}}\right|\left|{\boldsymbol{B}}\right|\right)\right)$ is the distance along the local ISMF line, ${\tau }_{1}=1/\left[2/{\tau }_{\mathrm{sc}}+1/{\tau }_{\mathrm{ex}}\right]$, and τ_sc = λ_∥/v is the PUI pitch angle scattering timescale with mean free path λ_∥. The NSW source function is

$$\begin{eqnarray}&&{S}_{\mathrm{NSW}}\left(r,\theta ,\phi ,v,{t}^{{\prime} }\right)=\displaystyle \frac{{I}_{\mathrm{NSW}}\left(r,\theta ,\phi ,v,{t}^{{\prime} }\right)}{4\pi {v}^{2}}{m}_{{\rm{p}}}{n}_{{\rm{p}}}\left(r,\theta ,\phi \right){\sigma }_{\mathrm{ex}}v.\end{eqnarray} \tag{ 17 }$$

More details on the derivation of Equations (15) and (16), specifically $f\left(z\leqslant {z}_{A}\right)$ and τ₁ (which is the effective timescale for either charge exchange or pitch angle scattering to occur), can be found in Schwadron & McComas (2013) and their application in our global MHD simulations in Zirnstein et al. (2019a, 2021a).

In our past work, we assumed that the rate of PUI pitch angle scattering, 1/τ_sc, outside the retention region was either much faster than charge exchange or much slower than charge exchange, with the former yielding a broad background of ENA flux, like the GDF, in directions far from the ribbon (see Figures 6 and 7, second rows, in Zirnstein et al. (2019a)). The intent of the current study is to compare our simulation results to the separated ribbon data from Swaczyna et al. (2022), who separated out the GDF, which may have a component originating from outside the heliopause. Therefore, we present simulation results for the case of slow pitch angle scattering outside the retention region, which does not produce a GDF-like signature outside the ribbon region. In this limit where τ_ex << τ_sc, Equation (16) approaches zero. We do not expect that our main results would change significantly if we simulated fast pitch angle scattering outside the retention region and then post-processed the modeled ribbon by removing the "background" similarly to the GDF removal by Swaczyna et al. (2022). This is because the key characteristics of the retention model are the enhancement in PUI density in the middle of the retention region and the closer proximity of the ribbon source to the heliopause than in the weak scattering case.

For the case where PUIs are created from the NPI distribution, which is much hotter than the NSW distribution, we must integrate over the NPI distribution similarly to the weak scattering case. Thus, inside the retention region,

$$\begin{eqnarray}&&\begin{array}{l}{f}_{0,\mathrm{NPI}}(r,\theta ,\phi ,v,t)={\displaystyle \int }_{0}^{{t}_{\max }}W({t}^{{\prime} })\\ \ \times \left[\displaystyle \frac{1}{4\pi }\iint \left({S}_{\mathrm{NPI}}\left({{\boldsymbol{r}}}^{{\prime} },{\rm{\Omega }},v,{t}^{{\prime} }\right){\tau }_{\mathrm{ex}}+f\left(z\leqslant {z}_{A}\right)\right)d{\rm{\Omega }}\right]{{dt}}^{{\prime} },\end{array}\end{eqnarray} \tag{ 18 }$$

where we integrate over solid angle Ω in phase space. It should be noted, again, that we exclude any creation of PUIs whose initial velocity vectors would not originate within the HTS (see the end of Section 2.2.3). The NPI source function ${S}_{\mathrm{NPI}}\left({{\boldsymbol{r}}}^{{\prime} },{\rm{\Omega }},v,{t}^{{\prime} }\right)$ is the same as Equation (14), except we combined Ω = (ω, φ) for simplicity.

The solar cycle averaged maps in Figure 2 showing reasonable similarity between the data and weak scattering model, in both intensity and structure. However, both models clearly underestimate the flux from the ribbon knot (high northern latitudes at energies >1 keV), likely due to our assumption for the SW speed variability as discussed in Section 2.2.1 and Equation (2). Both models also show relatively high intensity close to the ecliptic plane at energy ∼0.7 keV compared to mid latitudes, while the data show a more uniform distribution of flux at all latitudes at this energy. The reason for this discrepancy is unclear, but is likely related to either (1) inaccuracies in the ribbon separation and/or (2) inaccuracies in the SW speed model derived from IPS observations. First, the IBEX maps at ∼0.7 keV that include both the GDF and ribbon flux also show the highest intensities along the ribbon near the ecliptic plane, but this effect may have been removed in the ribbon-separated maps due to a localized GDF maximum near the nose coinciding with this portion of the ribbon. Second, while we shift IPS-derived SW speeds near the ecliptic plane to match OMNI, this shift is performed uniformly at all latitudes. If we trust the accuracy of the OMNI database of observations, then it is possible the model ribbon fluxes at mid-to-high latitudes are underestimated at energies near ∼0.7 keV (or SW speeds near ∼350 km s⁻¹). Either way, this is an example of the issues with understanding high-latitude SW speeds from IPS without in situ measurements.

Next, we show the results of our time-dependent ribbon models from 2009 to 2019. The results are calculated in the same frame of the spacecraft (specifically the "ram" frame) and the same measurement time per longitudinal swath as Earth, and thus IBEX, orbited the Sun each year (McComas et al. 2010). The model results are also integrated over the IBEX ESA response functions (Funsten et al. 2009). Figure 3 shows results for the ribbon at ESA 3 (∼0.8–1.6 keV) for annual ram maps from 2009 to 2019. In general, the agreement between the models and data are good, where there is a gradual decrease in ENA intensity from ∼2011 to 2016, after which the fluxes level out. Then, in 2019, the southern portion of the ribbon appears to brighten in both the data and model results. This is likely a direct response to changes in the SW dynamic pressure (both speed and density) in the middle of solar cycle 24 (McComas et al. 2020). There is not much structural difference between the two models, at least at 6° pixel resolution (Zirnstein et al. 2019a, 2019b), besides the significant difference in total intensity.

Zoom In Zoom Out Reset image size

Figure 4 shows results for the ribbon at ESA 5 (∼2.0–3.8 keV). Again, there is general agreement between the data and model results, although the difference in intensity between the spatial retention model and the data is even more pronounced. This implies that the discrepancy between the data and spatial retention model intensities is energy dependent. Notwithstanding this issue, there is clearly a decrease in ribbon intensity from ∼2012 to 2016, and a hint of brightening occurring in the weak scattering model south of the ecliptic plane but not visibly in the spatial retention model.

Zoom In Zoom Out Reset image size

We have not yet discussed how the NSW and NPI sources contribute to the overall intensity and structure of the ribbon in the models. Due to their significantly different distributions in energy (Figure 1), we should expect that they have energy-dependent effects on the ribbon. In Figure 5, we show sky maps of the weak scattering ribbon model in 2009 for each ESA, where we separate the components of the ribbon from the NSW (left column) and NPI (middle column), alongside the total (right column). From Figure 5, one can see that the NSW dominates the source of the ribbon signal at ESA 2 and 3, because the NSW flux has the highest intensity near the SW speed which, at low-to-mid latitudes, is near the ENA energies measurable by ESA 2 and 3 (∼0.5–1.6 keV). At ESA 4, we see an increased contribution from the NPI to the ribbon source almost uniformly at all latitudes, due to its large thermal spread. At ESA 5, the NSW only appears to contribute to the knot at high northern latitudes, while the NPI dominates the rest of the ribbon structure. At ESA 6, the NPI source completely dominates the ribbon signal in the sky. The highest SW speeds are typically ∼750 km s⁻¹, which lies in the middle of ESA 5. Therefore, the hotter NPI fluxes are the dominant source of the ENA signal observed in ESA 6.

Zoom In Zoom Out Reset image size

Figure 6 shows results similar to those in Figure 5, except for the spatial retention model of the ribbon. While the relative contribution of the NSW versus NPI to the source of the ribbon is similar to that in the weak scattering model, the NPI component is more broadly spread out in the sky and does not exhibit a narrow ribbon-like structure. This behavior may be consistent with part of the GDF or the wider ribbon observed at higher energies (Fuselier et al. 2009; Schwadron et al. 2011). Similarly to the weak scattering model, here the NSW dominates the ribbon knot seen in ESA 5, and the NPI component dominates the entire ribbon signal at ESA 6.

Zoom In Zoom Out Reset image size

The energy dependence of the ribbon flux with the inclusion of NPI is illustrated in more detail in Figure 7. We extracted the model ribbon fluxes from a pixel near ecliptic longitude 256° and latitude −23°, i.e., about 28° below the nose of the heliosphere. The models without NPI, shown in blue, exhibit steep drop-offs in intensity at energies above 1 keV. With the inclusion of the NPI source, the ribbon spectra flatten in the IBEX-Hi energy range, yielding a power-law slope in energy of ∼E⁻³. This illustrates the importance of the NPI component for the source of the IBEX ribbon at energies >2 keV, for any secondary ENA ribbon source mechanism.

Zoom In Zoom Out Reset image size

In this section, we show the time series of the observed and modeled IBEX ribbon signal averaged over three regions of the sky. The three regions, as illustrated in Figure 8, are: (1) the "knot", where the highest intensity of the ribbon is observed at ESA 5 and 6; (2) the "ecliptic", which is centered around the ecliptic plane (which is also close to the solar equatorial plane, i.e., an angular difference of ∼1 IBEX pixel); and (3) the "south", which is the southern region of the ribbon. Each of these regions has a width (i.e., angular size perpendicular to the ribbon's circle) of ∼20°. Time series of the observed and modeled ribbon fluxes averaged over these three regions in the sky are shown in in Figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left, middle, and right columns show time series of ribbon fluxes from the knot, ecliptic, and south regions of the sky, respectively. We normalize the fluxes such that each time series (per ESA and region) has the same average over time. The factors used to scale the data and model results are shown next to the labels of the respective time series curves. We will discuss the significance of these scaling factors near the end of this section. We show modeling results where we shift the IPS-derived SW speeds by a value to match OMNI near the ecliptic plane, and this shift is performed uniformly at all latitudes. We also show results where no shift is performed.

In general, there is good agreement between the models and the data in terms of overall evolution with the solar cycle. First focusing on the ecliptic region results, in general there is a decline in ENA fluxes over all energies, with the start of a recovery beginning after 2018, depending on the ENA energy. Comparisons with data become difficult at ESA 5 and 6 in the ecliptic region, because there is little ribbon flux at low latitudes and data uncertainties of the separated ribbon are large. The models predict the greatest relative changes occur at the highest ESA, illustrating the energy-dependent influence of the SW as it changes over the solar cycle. It is important to note that there are stark contrasts in the models with and without modifying IPS-derived speeds to match OMNI. Clearly, if there is a difference between IPS speeds and in situ measurements in the ecliptic plane, then IPS speeds must be shifted to match those from in situ measurements. The OMNI-shifted modeling results compare moderately better to the data at all ESA, except possibly for ESA 2. Moreover, there is little difference between the weak scattering and strong scattering results, at least in terms of relative evolution over time. If there are differences, they are typically smaller than the uncertainties of the data and mostly in total intensity (see scaling factors).

Next, looking in the "knot" region which encompasses ecliptic latitudes 20°–70°, there is very good agreement in the temporal trend between the OMNI-shifted model results and the data, except for ESA 2. Again, it appears there is an issue in our model of the ribbon at ESA 2 energies over time, especially at high northern latitudes. The reason for this is unclear, but it is likely related to the NSW flux model. At higher energies, both models match the data very well. It also appears that shifting the IPS-derived SW speeds makes for better agreement with IBEX data. This is particularly noticeable at ESA 3, where the increase in ENA flux in the "no scale OMNI" model is clearly disjointed from the observations.

Finally, looking in the "south" region, which encompasses ecliptic latitudes −20° to −40°, there is also good agreement between the models and the data. Interestingly, the models reproduce ESA 2 observations quite well in this region, in particular for the models not adjusted to OMNI. There are more scattered IBEX data points in the south region, in particular at ESA 5 and 6, possibly due to uncertainties in the ribbon separation that were potentially unaccounted for. For example, the increase in ESA 5 ribbon flux observed in 2011 and 2012 is correlated with a significant reduction of the surrounding GDF (Swaczyna et al. 2022), which may have resulted in an overestimation of the ribbon flux. Again, it is impossible to determine which model of the ribbon compares better to the data in terms of temporal trends. Interestingly, unlike the results in the knot region, the "not shifted to OMNI" results appear to match IBEX data slightly better than the shifted results, in particular after 2016, which may be linked to the large increase in SW dynamic pressure in late 2014, which affected ENA emissions observed at 1 au starting in late 2016 (McComas et al. 2018b).

To summarize an important point: the requirement of adjusting IPS-derived SW speeds to match OMNI is obviously needed in the ecliptic region, because "ground truth," in situ measurements are more reliable. In the knot region, it appears that the shifting in speed is typically better than not shifting, at least within the data uncertainties. But in the south region, there is a stark contrast in the temporal trend of the modeling results after 2016, and comparisons with the data imply shifting the IPS-derived SW speeds makes the comparison slightly worse. This implies that the discrepancy between IPS-derived SW speeds and in situ measurements cannot be assumed to be the same at all latitudes, and in fact perhaps the derivation of SW speed from IPS data using density fluctuations (Tokumaru et al. 2021) may be different in the northern and southern hemispheres of the Sun. The physical reason for this is unknown, but it should be considered in future studies of the SW observations from OMNI and speeds derived from IPS observations.

While the normalized temporal trends of the data and modeling results generally agree overall, we note that the total intensities of the modeling results differ significantly between the weak and strong scattering models, as is clear in the labels in Figure 9. The strong scattering model fluxes are approximately a factor of 4 lower in intensity than the weak scattering model, and the discrepancy is larger at higher ENA energies. While the implementation of the spatial retention model by Schwadron & McComas (2013) showed good agreement with observations (see also Schwadron et al. 2018), our implementation of their methodology in a global MHD simulation, especially in the current study, yields significantly different ribbon intensities. The difference in intensities between the implementations in the current study (also in Zirnstein et al. 2021a) and that by Schwadron & McComas (2013) is at least partially due to the fact that our global MHD simulation yields self-consistent plasma and neutral H densities, as well as the compression and draping of the ISMF, whereas Schwadron & McComas (2013) used simplified assumptions for the interstellar densities and ENA source region thickness to compute the model intensities. However, a more important difference is likely the NSW and NPI flux source to the ribbon. With a smaller spread in SW speed with distance from the Sun in our model (typically δ v < 50 km s⁻¹ in Equation (1), versus δ v = 250 km s⁻¹ for Schwadron & McComas 2013), a larger δ v will yield higher intensities at energies >1 keV. Zirnstein et al. (2019a), who first implemented the spatial retention model in a global MHD simulation, assumed δ v = 100 km s⁻¹ at all distances from the Sun, yielding better agreement with IBEX data than the current study. This discrepancy demonstrates how crucial it is to understand the dynamic properties of the SW, as a function of both time and space, which motivated our use of NH SWAP observations to constrain the radial dependence of δ v.

One of the most important characteristics of the ribbon that not only informs us of the draping of the ISMF around the heliosphere (Zirnstein et al. 2016b) but also may discriminate between microphysics of PUIs that become secondary ribbon ENAs (Zirnstein et al. 2021a) is the ribbon's geometric shape across the sky; specifically, the ribbon's center and radius as a function of energy and time, assuming the ribbon is circular (Funsten et al. 2013; Dayeh et al. 2019). We present results of the observed and modeled ribbon's center location and radius as a function of energy and time in Figures 10 and 11, respectively. A description of the method for calculating the ribbon center and radius, as well as their uncertainties, is similar to the work by Dayeh et al. (2019) and Zirnstein et al. (2021a), and is summarized in the Appendix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown in Figure 10, both models and the data ribbon centers' exhibit a large spread in latitude of ribbon centers as a function of ENA energy, with centers at larger energies observed at lower latitudes. The centers of the ribbon-separated results from spherical harmonic decomposition (Swaczyna et al. 2022) are qualitatively similar to earlier analyses of determining the ribbon center from the total ENA flux in the sky (Funsten et al. 2013; Dayeh et al. 2019), in terms of the trend in latitude as a function of energy. It is also consistent with the theoretical model of the ribbon centers from Swaczyna et al. (2016), who showed that the latitudinal asymmetry of the SW speed likely creates the spread in latitude of the ribbon centers as a function of ENA energy.

Dayeh et al. (2019) performed the first analysis of the time-dependent evolution of the IBEX ribbon centers, which we compare to here. Again, while qualitatively similar, the change in observed ribbon centers from our current analysis over time is much more significant than the results from Dayeh et al. (2019), presumably due to the ribbon separation before performing the circle fitting. Without ribbon separation prior to fitting, Dayeh et al. (2019) show that the ribbon centers do not significantly change over time, within their 1σ uncertainties (see their Figure 6). We also note that the ribbon center uncertainties from Dayeh et al. (2019) were assumed to be the same as the radius uncertainties, which are typically larger than our center uncertainties, which are derived from the minimized chi-square deviates. In our analysis, we see more significant changes over time in both models and the data results. For example, the latitudinal spread in centers is largest in 2009–2011 and 2018–2019, and in 2014–2015 the data show a clockwise rotation in the ribbon centers compared to other years. The data in 2009–2011 and 2018–2019 reflect solar minimum conditions (with fast SW at high latitudes and slow SW at low latitudes), while in 2014–2015 the SW conditions were more symmetric in latitude. Zirnstein et al. (2016a) showed that, under the assumption of uniform SW speed in latitude, the ribbon centers do not show significant spread in latitude, but rather become more aligned with the B - V plane. This explanation appears consistent with the data results shown in Figures 10(c) and (h) in 2014–2015. Thus, it appears that the ribbon centers as a function of energy, when derived after separating the ribbon signal from the GDF, do change over time and reflect the SW latitudinal characteristics over the solar cycle.

The modeling results show similarities to the data, although there are significant departures in 2014–2017 for the weak scattering model, especially at 2.7–4.3 keV. The spatial retention model better reproduces the data, but there is an almost time-independent, systematic ∼2°–3° shift in the model centers from the data along the B – V plane. It is possible that, with the larger plasma grid size in the current simulation (with radius 2500 au compared to 1000 au in Zirnstein et al. 2016b), the draping of the ISMF around the heliosphere is slightly different. We have determined that more ENAs from the heliotail in our 2500 au simulation are able to travel upstream of the heliopause and slightly modify the plasma conditions upstream of the bow wave by slightly increasing the temperature. This may then slightly alter way the ISMF is compressed and draped in the ribbon source region. While this effect appears small (a few degrees), this will be the subject of a future investigation.

Figure 11 shows the ribbon radius as a function of ENA energy and time for the models and data. In general, the weak scattering model results are qualitatively similar to the data results, but the model systematically overestimates the observed ribbon radius by a few degrees. As explained by Zirnstein et al. (2021a), the reason why the weak scattering model produces a larger ribbon radius is because the mean ENA source distance is farther from the heliopause than in the spatial retention model. Zirnstein et al. (2021a) explained that, in order to fix this discrepancy, the pristine ISMF magnitude would need to decrease by ∼1 μG, which is a significant change and highly unlikely to be realistic. The spatial retention model compares very well to the data at all energies and times. The observed ribbon radius at 4.3 keV in 2018–2019 is smaller by almost 4° than the spatial retention model, but we note that the separated ribbon at 4.3 keV has significantly low signal, and uncertainties in the fitting procedure may not be accounted for in these plots. Therefore, we believe the spatial retention model of the ribbon radius overall agrees with the observations.

While the ribbon centers and radii show significant evolution in time for different ENA energies, Figure 12 demonstrates that the energy-averaged center and radius of the ribbon do not significantly change over time. Both the models and the data average center results lie on or close to the B – V plane, consistent with our previous analysis showing that the average ribbon center lies along the same plane defined by interstellar magnetic field B and the LISM inflow vector V far from the heliosphere (nominally chosen at 1000 au; Zirnstein et al. 2016b). There are only slight differences between the weak and strong scattering model centers, but a significant difference in their radii, as explained earlier. The observed ribbon centers do not appear to change over time, within uncertainties, and the spatial retention model results lie within the observed radius uncertainties over the solar cycle. We note that the energy average does not include ESA 6 results, because the fitting analysis of the observations did not yield meaningful results from 2012 to 2017, due to the culling criteria as described in the Appendix (which mostly affected ESA 6), and therefore we left out ESA 6 from the averages in both models and the data.

Zoom In Zoom Out Reset image size