SIMULATIONS OF A DYNAMIC SOLAR CYCLE AND ITS EFFECTS ON THE INTERSTELLAR BOUNDARY EXPLORER RIBBON AND GLOBALLY DISTRIBUTED ENERGETIC NEUTRAL ATOM FLUX

For every iteration of the MHD-plasma code, the SW inner boundary conditions are updated based on a dynamic solar cycle. First, we assume that the SW boundary conditions follow an idealized solar cycle, similar to Pogorelov et al. (2009a), where the transition between fast and slow SW varies sinusoidally with an 11 yr period, the transition approaches ±35° at solar minimum and ±80° at solar maximum, the current sheet tilt angle goes from 0° at minimum to 90° at maximum, the magnetic field polarity reverses at solar maximum, and the fast and slow SW speeds (800 and 400 km s⁻¹, respectively) and densities (3.6 and 8 cm⁻³, respectively) are held constant. Between each iteration of the MHD-plasma code, the kinetic-neutral code is also run, coupling the neutral H atoms to the evolving plasma through source terms (see Appendix). The coupled MHD/kinetic code is run for a sufficiently long time to produce a full, time-dependent solution of the heliosphere and to fill the IHS with material from previous solar cycles. In practice, the simulation runs for more than 500 yr, such that the slow LISM plasma and neutrals have sufficient time to interact with the heliosphere and propagate through the simulation domain.

For the last 28 yr of the simulation, we use new, time-dependent values for the fast and slow SW speeds and densities. These values, beginning at solar maximum in 1990.5 (see Figure 2), are introduced into the simulation initially at a time corresponding to solar maximum in the simulation. The SW boundary conditions are then updated at every iteration of the simulation for the next 28 yr, as shown in Figure 2. These boundary conditions are calculated from time- and latitude-dependent functions of SW speed and density from Sokół et al. (2013), which are based on line-integrated interplanetary scintillation and in situ measurements of the SW, scaled to 1 AU. We take values for the SW speed and density at heliolatitudes of 0° (equator) and +90° (north pole) using Equations (3) and (4) from Sokół et al. (2013), as a function of time, and use these for the SW boundary conditions for the slow and fast SW regimes, respectively (see Figure 2). We assume that the fast and slow SW still maintain uniform speeds and densities at the inner boundary, in their respective regions, although varying as a function of time. The transition between fast and slow SW, the current sheet tilt angle, and magnetic field polarity still vary according to the method from Pogorelov et al. (2009a). The MHD-plasma and kinetic-neutral codes continue to iterate in order to couple the neutral H with the plasma. These last 28 yr correspond to times between 1990.5 and 2017.5, during which we compute H ENA fluxes at 1 AU in the post-process step (Figure 1, red boxes), which is described in Section 2.2.

Zoom In Zoom Out Reset image size

It would be ideal to incorporate the full latitudinal profile of Equations (3) and (4) from Sokół et al. (2013) as boundary conditions in our MHD/kinetic simulation. However, since we connect the last 28 yr of the simulation (with dynamic boundary conditions; see Figure 2) to the simplified solar cycle with a fast/slow SW transition that varies sinusoidally with an 11 yr period, we resort to the boundary conditions discussed above in order to preserve this frequency. This method, however, may not provide the most realistic distribution of ENAs produced in the supersonic SW, which is important for the secondary ENA mechanism. Therefore, we can improve the ENA distribution in the supersonic SW by scaling the speeds of those ENAs generated in our kinetic-neutral module by the latitude-dependent speeds from Equation (3) of Sokół et al. (2013). For example, an ENA generated in the supersonic SW in our kinetic-neutral module is given a speed calculated by $v_{{\rm EQ}3}^{1{\rm AU}}(\theta )\times {{v}_{{\rm MHD}}}(r,\theta )/v_{{\rm MHD}}^{1{\rm AU}}(\theta )$, where $v_{{\rm EQ}3}^{1{\rm AU}}(\theta )$ is the speed determined by Equation (3) from Sokół et al. (2013), ${{v}_{{\rm MHD}}}(r,\theta )$ is the local speed at position r determined by charge exchange with the MHD plasma, and $v_{{\rm MHD}}^{1{\rm AU}}(\theta )$ is the speed determined by charge exchange with the MHD plasma at 1 AU, all at the same latitude θ. This method was also used by Heerikhuisen et al. (2010a, 2012, 2014) and Zirnstein et al. (2013), providing a simple way to improve the latitudinal spectrum of the secondary ENA ribbon, while also taking into account the slowing of the supersonic SW due to the incorporation of PUIs. The probability for the charge-exchange event to occur, however, is still self-consistently determined by the local plasma properties from the MHD/kinetic solution of the heliosphere. Elsewhere in the heliosphere (IHS, OHS, and LISM), the speeds of newly created ENAs are self-consistently determined from the local plasma properties.

The data taken from Sokół et al. (2013) only span from 1990.5 to 2011.5; therefore, we must estimate the SW boundary conditions beyond 2011.5 in order to provide us with a background heliosphere for which to self-consistently calculate ENA flux from, beyond 2011.5. First we notice that the slow SW speed has, on average, remained fairly constant since 1990.5 (although with quasi-periodic oscillations around the mean). Therefore, for the slow SW speed beyond 2011.5 we reverse the profile from 2011.5 to 2006 and substitute it for the profile for 2011.5 to 2017.5. For example, the slow SW speed in 2011.75 is the same as 2011.25, 2012 is the same as 2011, etc. Although the slow SW density has decreased from 1990.5 to ∼2005, it has remained fairly constant since then. Therefore, we do a similar reversal for the slow SW speed.

The fast SW speed and density, however, are more complicated. As the solar cycle approaches maximum, the speed (density) near the poles decreases (increases) closer to the slow SW value as the coronal holes close. However, we again choose to reverse the profile of fast SW speed and density beyond 2011.5, for the following reasons. First, since the latest solar maximum was significantly different from the previous maximum, we can not accurately predict the fast SW behavior. Second, the maximum (minimum) in fast SW density (speed) during the previous solar maximum occurred in ∼2000.5, which is a full cycle before 2011.5, suggesting that a reversal at 2011.5 would be sufficient. Third, while we may assume some form of the fast SW profiles beyond 2011.5, we do not need accurate SW boundary conditions beyond 2011.5 to predict the ribbon flux for the next ∼3.5–9.5 yr, which will be discussed in more detail in Section 3.3. Therefore, for the sake of simplicity we reverse the profiles in 2011.5 as was done for the slow SW. The point of reversal for slow and fast SW in 2011.5 is indicated by the vertical dashed lines in Figure 2.

We note that it takes time for the new boundary conditions taken from Sokół et al. (2013) to propagate through the simulation domain, and this may have adverse effects on our ENA flux results. Therefore, we consider its potential impact in a few simplified cases. The earliest time that we simulate flux at 1 AU is 2005.5 (see Figures 15, 16, 19, 20). For ENAs in the lowest IBEX-Hi energy passband, with energy ∼0.71 keV and speed ∼369 km s⁻¹ (the center of passband 2), the time it takes to propagate from 1000 to 1 AU is ∼12.9 yr, from 200 to 1 AU is ∼2.6 yr, and from 100 to 1 AU is ∼1.3 yr. Assuming reasonable values for the SW speed (400 km s⁻¹ upstream, 100 km s⁻¹ downstream in the tail direction, and 50 km/s downstream in the nose direction) and TS distance (90 AU upwind, and 140 AU downwind), the SW plasma takes ∼42 yr to propagate from 1 AU down the heliotail to the simulation boundary, and ∼4.9 yr to propagate to the HP in the upwind direction (∼130 AU). Therefore, the plasma in the heliotail is not completely filled with the new SW conditions from Sokół et al. (2013) before we begin computing ENA fluxes at 1 AU (in 2005.5). Near the front of the heliosphere, however, there is a sufficient amount of time to fill the IHS before we begin H ENA flux computations.

However, plasma flowing toward the front of the heliosphere is slowed and diverted near the HP and takes a longer amount of time to flow to the poles and flanks and down the heliotail. Based on estimates of the transit times of SW through the IHS up to ENA detection from the poles (Dayeh et al. 2014; Siewert et al. 2014), we estimate a maximum transit time of ∼15 yr in the front half of the heliosphere. Since we compute flux after 2005.5 (at least 15 yr after introducing the new SW boundary conditions), our IHS flux results from the front of the heliosphere remain largely unaffected by the initial SW boundary conditions. However, near the downwind direction, our results may potentially include some ENAs generated by plasma from before the introduction of the SW boundary conditions from Sokół et al. (2013).

For the majority of the ribbon flux, the total time it takes to detect a secondary ENA from past SW conditions ranges between ∼3.5 and 9.5 yr, depending on the distance to the ribbon source and the ENA energy (see Section 2.2.2 and Table 2). Also, the amount of time it takes for the SW to fill the region upstream of the TS, which produces the majority of primary ENAs that fuel the secondary ENA ribbon, is approximately 1 yr. Therefore, we conclude that our simulated ribbon results, beginning in 2005.5 (15 yr after the introduction of SW boundary conditions from Sokół et al. 2013), are generally unaffected by the change in SW boundary conditions.

As discussed in Section 2.1.1, ENAs generated in the supersonic SW in the kinetic-neutral model during the last 28 yr of the MHD/kinetic simulation are given speeds scaled by the latitude and time-dependent results from Sokół et al. (2013). While this method is not initially self-consistent with the MHD/kinetic simulation, the resulting ENAs are still coupled to the MHD plasma by charge exchange. This method provides (1) a more realistic outward-propagating, neutralized, supersonic SW distribution; (2) more realistic heating in the OHS plasma by charge exchange; and (3) more realistic secondary ENA spectra. However, we note that this method does not include the convection time of SW from 1 AU to the location at which a new ENA is generated in the supersonic SW, which may introduce (1) a ≲0.5 yr (on average) discrepancy in our results, and (2) the ratio ${{v}_{{\rm MHD}}}(r,\theta )/v_{{\rm MHD}}^{1{\rm AU}}(\theta )$ may vary between ∼0.5 and 2 when computed near the fast/slow SW transition (e.g., when the SW at distance r from the Sun is slow [fast], and the SW at 1 AU, at the same point in time, is fast [slow]). While the effects of these simplifications have not been quantified, we note that (1) the resolution of our ENA flux results presented in Section 3 is 1 yr, and thus we do not expect a significant uncertainty when analyzing the evolution of ENA flux over multiple years; and (2) the fast/slow SW transition region is narrow most of the time and likely only affects a small portion of the outward-propagating, supersonic ENA flux.

Since there is a considerable delay in time between parent ENA creation in the supersonic SW and secondary ENA measurement at 1 AU (between ∼3.5 and 9.5 yr; see Table 2), we are able to predict ribbon flux consistently up until ∼2015–2021, depending on the energy of the ENA and the direction. Flux from the IHS, however, cannot be accurately predicted more than a few years into the future. Therefore, one must keep in mind that our results presented later in the paper for the IHS flux depend strongly on our assumptions of the SW boundary conditions beyond the currently available data.

The LISM boundary conditions for our time-dependent MHD/kinetic simulation are shown in Table 1, which we assume do not vary over the timescale of the simulation. The LISM plasma temperature (T), velocity (u), and magnetic field strength (B) were chosen based on consensus values derived from Schwadron et al. (2011), Möbius et al. (2012), Bzowski et al. (2012), and McComas et al. (2012a). The plasma and neutral densities were chosen in order to constrain the TS position to approximately agree with Voyager observations (although our simulation is not accurate enough to reproduce the time-dependent asymmetry of the V1 and V2 crossings), as well as to produce an H atom density at the nose of the TS of ∼0.09 cm⁻³.

Table 1.  LISM Boundary Conditions

Parameter	Value	Reference
T (K)	6200	1,2
np (cm−3)	0.07	⋯
nH (cm−3)	0.18	⋯
u (km s−1)	23.2	1–3
B (μG)	3	4
${\boldsymbol{u}}$ (ecliptic J2000)	(79°, −5°)	1–3
${\boldsymbol{B}}$ (ecliptic J2000)	(45°, −44°)	⋯

References. (1) Möbius et al. (2012); (2) Bzowski et al. (2012); (3) McComas et al. (2012a); (4) Schwadron et al. (2011).

Download table as:  ASCIITypeset image

From the mathematical perspective, we use characteristics-based boundary conditions at the outer boundary with additional modifications described by Pogorelov & Semenov (1997).

A recent study that looked at the trend of the interstellar He flow direction through the heliosphere, inferred from multiple sources of measurements over the past 30 yr, reported that a linear fit to the flow's ecliptic longitude direction shows a statistically significant increase over time (Frisch et al. 2013). However, we cannot accurately predict the evolution of the LISM conditions, and we assume that they vary slowly enough to not affect H ENA fluxes over the timescales presented in this paper.

There have been studies (see, e.g., Pogorelov et al. 2007, 2008a, 2009b, 2013; Richardson et al. 2008; Izmodenov 2009; Opher et al. 2009; Washimi et al. 2011) that have used modeling to constrain the interstellar magnetic field strength and direction needed to explain the asymmetry of the TS inferred by the Voyager spacecraft crossings (Stone et al. 2005, 2008), while other models have included a fit to the direction toward ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ outside the HP to the IBEX ribbon (see, e.g., Grygorczuk et al. 2011; Strumik et al. 2011; Ben-Jaffel & Ratkiewicz 2012; Ratkiewicz et al. 2012; Ben-Jaffel et al. 2013). The significance of the ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ surface in the generation of the IBEX ribbon was identified during initial IBEX observations (McComas et al. 2009a; Schwadron et al. 2009a), based on correlations between the ribbon direction and the ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ surface outside the HP (Pogorelov et al. 2009b, 2009c). In this paper we choose the strength and direction for ${{{\boldsymbol{B}} }_{{\rm LISM}}}$ based on recent analyses that fit simulations of the IBEX ribbon flux using the secondary ENA mechanism to the observed ribbon (Heerikhuisen & Pogorelov 2011; Funsten et al. 2013; Heerikhuisen et al. 2014), which are in general agreement with other studies that fit the IBEX ribbon to directions near ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$.

In these simulations we use an older estimate for the interstellar magnetic field direction, coming from (225°, 44°), rather than aligned with the more recent value from Funsten et al. (2013). Since it will take a considerable amount of time for any changes to the LISM boundary conditions to propagate through the entire computational grid, and the simulations require a significant amount of computational time, we use this older estimate for the LISM magnetic field in this paper, rather than assuming that the center of the ribbon is aligned with the interstellar field direction. Therefore, the direction of maximum flux from our simulated ribbon is slightly different from the observed ribbon, partially due to this difference. Heerikhuisen et al. (2014) also showed that the simulated ribbon center is offset from the simulated, pristine interstellar field direction by a few degrees, depending on the strength of the interstellar magnetic field (also see Grygorczuk et al. 2011). However, for the purposes of this paper, we only require an approximate direction for the interstellar magnetic field, such that the directions toward ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$, and thus the simulated ribbon, approximately coincide with the IBEX ribbon.

Recent V1 measurements of anomalous and galactic cosmic rays (Webber & McDonald 2013; Krimigis et al. 2013; Stone et al. 2013), as well as (indirectly) plasma density (Gurnett et al. 2013), have indicated that the spacecraft is outside the HP. However, initial measurements of the interstellar magnetic field orientation showed little change from a Parker spiral (Burlaga et al. 2013), although recent measurements of the field direction indicate a larger deviation (Burlaga & Ness 2014). Many studies have followed to explain these observations (Krimigis et al. 2013; Swisdak et al. 2013; Fisk & Gloeckler 2013, 2014; Florinski et al. 2013; Opher & Drake 2013; Schwadron & McComas 2013a; Borovikov & Pogorelov 2014; Gloeckler & Fisk 2014; Strumik et al. 2014). While the thickness of the heliosheath suggested by the recent crossing is smaller than in our simulation, the evolutionary trends in our ribbon results will be largely unaffected. We will discuss this further in Section 4.

In Figure 3 we show cross sections of plasma density computed using our time-dependent, 3D MHD/kinetic code, taken during solar minimum (left) and solar maximum (right). Notice the features similar to Pogorelov et al. (2009a), such as the transition between fast and slow SW near ±35° at solar minimum and ±80° at solar maximum, the sectored fast and slow SW densities in the heliotail, and the fluctuations in plasma density propagating through the OHS. We will discuss the effects of these fluctuations on the simulated ribbon flux in Section 4.

Zoom In Zoom Out Reset image size

The flux simulation starts at the time of measurement at 1 AU, ${{t}_{{\rm m}}}$, and integrates backward in time along the H ENA trajectory based on the speed of the H ENA we wish to simulate. The time of creation of an H ENA that will later be detected by the spacecraft is determined by the propagation time between the creation point and the spacecraft position. A primary H ENA from the IHS that is detected at 1 AU will travel anywhere from ∼90 (i.e., from the nose of the TS) to >1000 AU (i.e., from the extended heliotail), depending on the source location. For example, primary H ENAs with energies between ∼0.71 and 4.29 keV (corresponding to the nominal energies for IBEX-Hi passbands 2 and 6, respectively) will take between 1.16 and 0.47 yr for 90 AU and 12.9 and 5.23 yr for 1000 AU (we stress that these are only examples, not constraints in our simulations). However, as the SW flows away from the TS into the distant heliotail, high-energy PUIs that form a significant parent population for the ENA flux measured by IBEX-Hi will likely charge exchange within a few hundred AU from the TS (see, e.g., Zirnstein et al. 2014, Figure 3), thus reducing the effective time of propagation for the majority of high-energy ENAs from the IHS that are detected by IBEX. The time it takes the parent proton to advect with the plasma from the Sun to the position in the IHS before charge exchange occurs is an important characteristic for H ENA measurement and modeling (see, e.g., Allegrini et al. 2012; Reisenfeld et al. 2012; Dayeh et al. 2014; Siewert et al. 2014), and it is self-consistently accounted for in the MHD/kinetic simulation.

Secondary H ENAs with energies between ∼0.71 and 4.29 keV that propagate from the OHS (∼200 AU) to the spacecraft will take between ∼2.6 and 1 yr to travel to the spacecraft, respectively. However, secondary ENAs originally come from primary ENAs that crossed the HP, produced by charge exchange in the SW upstream or downstream of the TS. The propagation of these primary ENAs to the OHS is also taken into account in the MHD/kinetic simulation, and the propagation of secondary ENAs to 1 AU is taken into account in the post-process simulation.

As stated earlier, primary and secondary ENAs both have intrinsic propagation times that depend on the distance between their source location and the detector. However, we must also take into account the delay time between the creation of PUIs in the OHS and the creation of the secondary H ENAs (which form the ribbon). This is a considerable amount of time for all energies relevant for IBEX.

In Table 2 we show estimates of the times relevant to the production and detection of secondary ENAs from outside the HP, for the benefit of the reader. We stress that for our simulations we directly calculate the time-dependent flux as a function of space and time and do not use the values shown in Table 2. For the nominal energy in each IBEX-Hi passband we provide estimates of the mean free path of ENAs traveling through the OHS, the charge-exchange time for PUIs at the same energy, and the total time between parent ENA creation (i.e., approximately the time in the solar cycle that created the parent ENA) in the supersonic SW and secondary ENA detection at 1 AU.

As one can see, the time in the past that determined the energy and distribution of secondary ENAs ranges between ∼3.5 and 9.5 yr (for the majority of flux), largely depending on the energy of the ENA/PUI, the distance to the HP, and the OHS plasma and neutral properties. Although the parameters shown in Table 2 are not used in our simulations and are actually more complicated (for a space- and time-dependent heliosphere), these estimates are meant to assist the reader later in the paper.

While the source of PUIs that create secondary ENAs at a particular point in space and time extends from 0 to $\infty$ into the past, this is computationally unfeasible to simulate. However, within three times the mean charge-exchange time, approximately 95% of the source of PUIs are accounted for. Therefore, the flux of secondary ENAs created at a particular point in space and time is approximately determined by a source of PUIs from a time 0 to 3$\tau _{{\rm ex}}^{{\rm m}}$ in the past, where

$$\begin{eqnarray}&&\tau _{{\rm ex}}^{{\rm m}}={{({{n}_{{\rm H}}}v{{\sigma }_{{\rm ex}}}(v))}^{-1}}\end{eqnarray} \tag{ 1 }$$

is the mean charge-exchange time (when only 1/e particles have not experienced charge exchange yet), n_H is the background H density, v is the PUI velocity (i.e., the ENA velocity, which we assume is much greater than the bulk neutral speed), and ${{\sigma }_{{\rm ex}}}$ is the energy-dependent, charge-exchange cross section (Lindsay & Stebbings 2005).

Figure 4 shows the (normalized) probability of PUIs, as a function of time in the past, to be the source of secondary ENAs, for energies 0.71 and 4.29 keV, assuming that the flux of secondary ENAs is created at time 0. We calculate the total flux of secondary ENAs at time 0, at a specific location in the OHS, by integrating the flux of secondary ENAs created over the interval 0 to 3$\tau _{{\rm ex}}^{{\rm m}}$ in the past (i.e., creating secondary ENAs from parent primary ENAs over a range of time in the past), and we weight the contribution of flux over the interval by an exponential function of the mean charge-exchange time (see Equation (4)). The majority of secondary ENAs created at time 0 will originate from parent primary ENAs that charge exchanged within three times the mean charge-exchange time in the past (vertical dashed lines), where the mean charge-exchange time is a function of energy and neutral density (assumed to be 0.2 cm⁻³ for the example shown in Figure 4). For smaller energies, the charge-exchange time is longer, and thus we integrate over a longer time. Higher-energy PUIs charge exchange more frequently, therefore only requiring a shorter integration time.

Zoom In Zoom Out Reset image size

In order to test whether only accounting for PUIs within three times the mean charge-exchange time in the past is sufficient, we compared our results to simulations of ribbon flux at 2009.5 that accounted for up to only twice the mean charge-exchange time for each energy. In directions toward the ribbon, i.e., near ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$, and in particular the directions where we focus our analysis of the ribbon (see Figure 14, north, middle, and south), the percent difference between these two simulations was approximately within ±10%. In certain regions of the sky away from the ribbon, the deviations reached up to ±20–30% depending on the energy and direction, which is likely due to (1) areas of low accumulated flux, (2) rapid variations over time, or (3) low OHS neutral H densities (which will increase the mean charge-exchange time). Unfortunately, computational limitations prevent us from simulating the case for four times the mean charge-exchange time, since this amount of time will likely exceed the maximum time range of plasma-neutral results allocated in our post-process simulation. However, since the errors of our current simulations compared to a more accurate one can only be smaller than the errors reported above, we do not expect our approximations to significantly affect the results presented in this paper.

When we calculate the mean charge-exchange time, the variables from Equation (1) are determined at the time of secondary ENA creation, $t={{t}_{{\rm p}}}-{{t}_{{\rm m}}}$ (i.e., at time 0 in Figure 4), although it should be integrated as a function of time into the past, from t = 0 to $3\tau _{{\rm ex}}^{{\rm m}}$. However, first we assume in Equation (1) that the relative speed of charge exchange is dominated by the PUI speed such that u_H and ${{T}_{{\rm H}}}\simeq 0$, and that the PUI speed remains constant. Second, we assume that the bulk neutral density outside the HP remains nearly constant over the charge-exchange time period. Although not shown, we computed the deviation in the bulk neutral density in the OHS, over an 11 yr period near the middle of the ribbon, and found that the maximum deviation is approximately ±2% and drops to less than 1% over a 2 yr period. Thus, this assumption will have little impact on our results. Therefore, the mean charge-exchange time remains nearly constant over these periods of time, outside the HP, and only varies as a function of space.

Before charge exchange occurs, PUIs may stream along the interstellar magnetic field, depending on their initial pitch angle and the level of scattering. However, we ignore the effect of PUI streaming. Close to the ribbon, the parallel velocity (along the magnetic field) of PUIs is relatively small. For example, for a 1 keV PUI with an initial pitch angle of 80°, the parallel velocity (ignoring pitch-angle scattering) is ∼76 km s⁻¹. Over the course of 2 yr before charge exchange, the PUI will move ∼32 AU, which corresponds to a maximum angular spread of ∼9° (depending on ${{{\boldsymbol{B}} }_{{\rm LISM}}}$) as viewed from 1 to 200 AU, which is only slightly larger than the IBEX field of view (∼7°). If we include pitch-angle scattering, this may inhibit streaming if they scatter across 90°. Therefore, within the majority of the ribbon width of 20°, the movement of PUIs parallel to the interstellar field is not noticeably visible. However, this will be larger for faster PUIs, farther away from the ribbon, or if PUIs scatter to larger pitch angles, which may explain the large ribbon width at high energies (see, e.g., McComas et al. 2012b; Schwadron & McComas 2013b). On the other hand, PUIs with large pitch angles (assuming they do not scatter) will not be visible by IBEX, and therefore should not significantly affect our simulation results.

In their simulation of the ribbon, Chalov et al. (2010) assumed the no-scattering limit, allowing PUIs to stream along field lines until charge exchange occurred. In this limit, the magnetic moment is conserved such that the streaming of PUIs will slow near strong interstellar fields, allowing a significant number of PUIs to remain in the region of ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ before charge exchange occurs. As we showed in Table 2, secondary ENAs detected from different directions will have different time delays. If a significant number of PUIs that produced secondary ENAs near ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ actually originated far from that location, then the time evolution of ribbon flux will be much more complicated than our simulations can produce. However, while including this effect may produce a small, but noticeable, increase in flux near ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ (Chalov et al. 2010), most ribbon ENAs visible by IBEX will originate near ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ before charge exchange, suggesting that this effect will not be noticeable in the evolution of the ribbon flux and that our assumptions are reasonable.

We also do not take into account the convection of PUIs with the bulk plasma flow perpendicular to the interstellar field direction (see, e.g., Möbius et al. 2013). Using similar parameters as Möbius et al., the bulk motion of the PUIs would be ∼16.4 km s⁻¹, which over 2 yr moves ∼3.5 AU. Close to the HP, the bulk plasma flow would be even smaller; therefore, the effects of this motion are negligible for the purpose of this paper.

Using the characteristics of time delay for primary and secondary ENAs as described in the previous sections (see examples in Table 2), we can modify existing time-independent flux equations to calculate LOS integrated flux as a function of space (${\boldsymbol{r}}$) and time (t), in units (cm² s sr keV)⁻¹. For primary ENAs created inside the IHS, we use the differential primary ENA flux equation from Zirnstein et al. (2012, 2013), modified as a function of time, which is given by

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

where m_H is the mass of H, ${{f}_{{\rm p}}}$ is the (isotropic) proton phase space distribution in the plasma frame, ${{{\boldsymbol{v}} }_{{\rm p}}}$ is the parent proton (and resulting ENA) velocity in the plasma frame (i.e., ${{{\boldsymbol{v}} }_{{\rm p}}}={\boldsymbol{v}} -{{{\boldsymbol{u}} }_{{\rm p}}}$, where ${{{\boldsymbol{u}} }_{{\rm p}}}$ is the bulk plasma velocity), ${\boldsymbol{v}}$ is the ENA velocity in the Sun frame, n_H is the background H density, v_rel is the relative velocity between parent proton and background H distribution (see, e.g., Ripken & Fahr 1983; however, notice the corrected version in Heerikhuisen et al. 2006a), ${{{\boldsymbol{u}} }_{{\rm H}}}$ is the bulk H velocity, ${{\sigma }_{{\rm ex}}}$ is the charge-exchange cross section (Lindsay & Stebbings 2005), and ${\Delta }t$ is the integration, or flux accumulation, time step. For calculating primary ENA flux, we assume that that the background H distribution is Maxwellian, which is sufficient for the purposes of this paper (also see Zirnstein et al. 2013). The survival probability P is integrated as a function of space and time, where we only begin integrating losses to H ENAs outside 100 AU from the Sun. All of the plasma and neutral parameters are found at a time $t={{t}_{{\rm m}}}-{{t}_{{\rm p}}}$ in the past, where ${{t}_{{\rm m}}}$ is the time of measurement and ${{t}_{{\rm p}}}$ is the propagation time for ENAs from their source to the detector.

For simplicity, when calculating the H ENA flux we assume that the IHS plasma is represented by three Maxwellian distributions, one for SW ions, transmitted PUIs, and reflected PUIs (Zank et al. 2010). We assume constant relative proton densities and temperatures, similar to Zirnstein et al. (2013), using parameters from the V1 direction in Desai et al. (2014), where the relative densities are 80%, 17.5%, and 2.5% and the relative energy partitions are 3.93%, 40.82%, and 55.25% for SW ions, transmitted PUIs, and reflected PUIs, respectively. The IHS plasma is likely more complicated (see, e.g., Chalov et al. 2003; Malama et al. 2006; Wu et al. 2010; Zirnstein et al. 2014); however, we currently cannot accurately simulate the time dependence of PUIs in the IHS. Nevertheless, the point of this paper is to demonstrate the global effect of a time-dependent solar cycle on H ENA flux; therefore, we use this as an approximation. Losses to H ENAs traveling through the IHS are also self-consistently computed while assuming that the background IHS plasma is represented by a three-Maxwellian distribution, where we ignore losses inside 100 AU of the Sun and only take into account losses due to charge exchange outside 100 AU of the Sun. Losses by electron-impact ionization may be important in the heliosheath (Scherer et al. 2014), depending on the as- yet-unknown electron temperature; however, they are currently ignored in this study.

We also simulate time-dependent flux from outside the HP using the secondary ENA mechanism. Using the partial-shell method from Heerikhuisen et al. (2010b), the differential flux of ribbon ENAs (Zirnstein et al. 2012, 2013), modified as a function of time, is given by

$$\begin{eqnarray}\begin{array}{rcl} {\Delta }{{J}_{{\rm S}}}({\boldsymbol{r}} ,t) & = & \frac{1}{4\pi {{m}_{{\rm H}}}}{{n}_{{\rm p}}}({\boldsymbol{r}} ,t){{v}^{2}}P({\boldsymbol{v}} ,{\boldsymbol{r}} ,t) \\ {} & {} & \times \left( \int _{0}^{3\tau _{{\rm ex}}^{{\rm m}}}W({{\tau }_{{\rm ex}}}){{\int }_{{\Omega }}}\frac{{{f}_{{\rm H}}}({{{\boldsymbol{v}} }_{{\rm H}}},{\boldsymbol{r}} ,t)}{{{{\hat{v}}}_{{\rm H}}}\cdot {{{\hat{B}}}_{{\rm LISM}}}({\boldsymbol{r}} ,t)} \right. \\ {} & {} & \times {{v}_{{\rm rel}}}(|{{{\boldsymbol{v}} }_{{\rm H}}}-{{{\boldsymbol{u}} }_{{\rm p}}}({\boldsymbol{r}} ,t)|){{\sigma }_{{\rm ex}}}({{v}_{{\rm rel}}}({\boldsymbol{r}} ,t)) \\ {} & {} & \times \left. d{{{\Omega }}_{v}}d{{\tau }_{{\rm ex}}} \right){\Delta }t, \\ \end{array}\end{eqnarray} \tag{ 3 }$$

where ${{n}_{{\rm p}}}$ is the proton density, f_H is the distribution of outward-propagating ENAs outside the HP in the Sun frame (Heerikhuisen & Pogorelov 2010), integrated over velocity space, ${{{\boldsymbol{v}} }_{{\rm H}}}$ is the parent ENA velocity in the Sun frame (where ${{v}_{{\rm H}}}=v$), ${{{\boldsymbol{B}} }_{{\rm LISM}}}$ is the LISM magnetic field, v_rel is the relative velocity between parent ENA and background plasma distribution, and $d{{{\Omega }}_{v}}$ is the differential solid angle in velocity space. PUIs can only contribute to the ENA flux visible at 1 AU if their partial-shell distribution intersects the LOS; otherwise, the flux is not counted.

All of the plasma and neutral parameters are found at a time $t={{t}_{{\rm m}}}-{{t}_{{\rm p}}}-{{\tau }_{{\rm ex}}}$ in the past (except for n_H in Equation (1); see previous discussion), where ${{\tau }_{{\rm ex}}}$ is the charge-exchange time between PUI creation and secondary ENA creation (ranging between 0 and 3$\tau _{{\rm ex}}^{{\rm m}}$). We integrate over the contribution of flux from PUIs between 0 and 3${{\tau }_{{\rm ex}}}$ in the past, multiplying by the weighting $W({{\tau }_{{\rm ex}}})$ to normalize the total integral to 1. The weighting is given by

$$\begin{eqnarray}&&W({{\tau }_{{\rm ex}}})=\frac{{\rm exp} (-{{\tau }_{{\rm ex}}}/\tau _{{\rm ex}}^{{\rm m}})}{\int _{0}^{3\tau _{{\rm ex}}^{{\rm m}}}{\rm exp} (-\tau _{{\rm ex}}^{\prime }/\tau _{{\rm ex}}^{{\rm m}})d\tau _{{\rm ex}}^{\prime }}.\end{eqnarray} \tag{ 4 }$$

The survival probability for secondary H ENAs is computed similarly to ENAs from the IHS, except we assume that the background plasma distribution outside the HP is represented by a single Maxwellian distribution. While the total energy of the OHS plasma includes the effect of heating by PUIs from ENAs originating within the heliosphere (Heerikhuisen et al. 2008), the distribution of PUIs outside the HP is more complicated (Zirnstein et al. 2014), and a global simulation of their distribution is left for future studies.

Since the first release of IBEX observations (McComas et al. 2009a), several different models for simulating the ribbon flux from the secondary ENA mechanism have been developed (Heerikhuisen et al. 2010b; Chalov et al. 2010; Schwadron & McComas 2013b; Möbius et al. 2013; Isenberg 2014). Each model relies on an assumption of the evolution of the PUI distribution in the OHS, where Heerikhuisen et al. (2010b) assumed that PUIs scatter onto a partial shell, Chalov et al. (2010) and Möbius et al. (2013) assumed the no-scattering limit, and Schwadron & McComas (2013b) and Isenberg (2014) assumed that PUIs are spatially confined by wave turbulence and scatter onto nearly isotropic shells. Assuming that the energies of the PUIs do not change appreciably before charge exchange occurs and the PUIs do not move far from their initial positions, the evolution of secondary ENA flux will be qualitatively similar for all of these models. Therefore, the results of this paper apply to the general secondary ENA mechanism.

Figure 7 shows the percent change in simulated flux for passband 2 (∼0.71 keV). Since 2009.5, the simulated ribbon flux has been decreasing significantly at mid-latitudes. The time in the solar cycle that fueled the 0.71 keV ribbon is approximately 6.5–9 yr earlier, depending on the direction (see Table 2). Therefore, after 2010.5, the ribbon was fueled by primary ENAs created between ∼2001.5 and 2004, which was near the beginning of the decline of the previous solar maximum. During this phase of the solar cycle, the coronal holes are opening up at lower latitudes, increasing in fast ($\geqslant$2 keV) SW (and thus fast neutral SW), and decreasing in slow ($\leqslant$2 keV) SW. This causes the decline in ribbon flux at passband 2. Behind the ribbon, the IHS flux also shows variations every year, although partially drowned out by the decrease in ribbon flux. There is a general decrease in flux near the poles, fluctuations near the front and flanks, and little change in the heliotail.

Zoom In Zoom Out Reset image size

Also notice the bright, nearly horizontal streaks of flux at high (positive and negative) latitudes. These enhancements in flux come from ENAs produced in the IHS where the transition region between fast and slow SW creates two "cusps" of heated, subsonic plasma downstream of the TS that penetrate into the supersonic SW (see Figure 8). As the solar cycle transitions from minimum to maximum, these regions move to higher latitudes and eventually diminish (also see Figures 9–12, as flux from the cusps moves upward over time). While this feature may seem to be artificial due to our assumptions of a "sharp" transition between fast and slow SW, it is known that the fast and slow SW boundary is relatively sharp at solar minimum (see, e.g., McComas et al. 2000); therefore, as the solar cycle transitions from minimum to maximum, this feature may be visible in ENA spectra observed by IBEX, although possibly drowned out by other sources. These features are hardly visible, however, in the time-averaged results shown in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the percent change in simulated flux for passband 3 (∼1.11 keV). Similar to 0.71 keV, the ribbon flux decreases at mid- to high latitudes. Also, there is little change in ribbon flux near the ecliptic plane where the average SW speed is relatively constant. At this energy, the percent change in ribbon flux decreases more dramatically at high latitudes compared to 0.71 keV (compare bottom right of Figures 7 and 9). The IHS flux is similar to that at 0.71 keV.

Note that the ribbon shape is barely visible in these maps since the percent change of ribbon flux is nearly uniform, reflecting the nearly uniform decrease in neutralized SW as a function of latitude (also see McComas et al. 2014a). If we were to look at the evolution of absolute flux over time, the largest change would occur in the direction of the ribbon. Also, since the ribbon flux is strongest near the front of the heliosphere, where the HP is closest, we only see a decrease in ribbon flux near the front, and less change from the back.

Figure 10 shows the percent change in simulated flux for passband 4 (∼1.74 keV). At this energy, the ribbon flux begins to increase over time at mid-latitudes and slightly move to higher latitudes. With a total time delay of approximately 6 yr, the source of the ribbon flux measured between 2010.5 and 2013.5 was created between ∼2004.5 and 2007.5. Over this time, the coronal holes were still opening, producing more fast wind, and thus causing an increase in the low-latitude ribbon flux at this energy. Near the poles, there is still a decrease in flux, from both the IHS and OHS (also see Figure 13).

Zoom In Zoom Out Reset image size

Figure 11 shows the percent change in simulated flux for passband 5 (∼2.73 keV). Similarly to 1.74 keV, the ribbon flux increases as the source of fast neutral SW opens up toward solar minimum. The largest increase in flux occurs near the simulated ribbon knot, which appears in our simulation near (340°, 35°) at this energy and time. Notice that this is lower in latitude than the observed knot, possibly due to our use of simplified SW profiles that determine the outward-propagating ENA speeds. Nevertheless, the knot increases in flux over this period of time, more so than the ribbon flux at negative latitudes. The average fast SW speed is near 2.7 keV; therefore, passband 5 is a good indicator for the evolution of the fast SW at mid- to high latitudes.

Zoom In Zoom Out Reset image size

Figure 12 shows the percent change in simulated flux for passband 6 (∼4.29 keV). There is an increase in flux at high latitudes, with the largest increase where the ribbon knot (although faint; see Figure 5) is located near (0°, 45°), where 0° longitude is the fourth longitude line to the left of the nose. At this energy, the delay in time for secondary ENAs from outside the HP is approximately 4.5 yr; therefore, at a measurement time of ∼2013.5, this flux approximately reaches a maximum and will begin to decrease (see Figure 16). The evolution of the IHS flux is also interesting, where we see the two streaks increasing in flux at high positive and negative latitudes, moving to higher latitudes over time. The delay in time for the IHS flux at this energy is within a couple years (depending on the fast or slow SW and the direction of detection); therefore, the source of IHS flux between 2010.5 and 2013.5 is approaching solar maximum, where the transition between fast and slow SW creates two cusps of hot, subsonic plasma downstream of the TS that penetrate into the supersonic SW and move to higher latitudes (see Figure 8). These horizontal streaks are partially visible in Figures 7–12 but are quite distinct at this energy, when the relative flux from the ribbon is not as high.

Zoom In Zoom Out Reset image size

After initially increasing in intensity (after 2012.5 and 2013.5 for 4.29 and 2.73 keV, respectively), the knot of emission at these energies also moves to higher latitudes over time while diminishing in intensity, another signature of the evolution of the fast SW toward higher latitudes (i.e., the closing of the coronal holes approaching solar maximum).

Next, we compare the results shown in Figures 7–12 to IBEX data shown in Figure 20 of McComas et al. (2014a). While our results show a continuing decrease in flux for passbands 2 and 3 over a significant portion of the sky (Figures 7 and 9), similar to the data, at higher energies the data also show a global decrease in flux, whereas the simulations show an increase, particularly in the ribbon flux. The simulation results also show that little change occurs in flux from the heliotail, for all energies. While there appears to be a mixture of increases and decreases in flux across the heliotail region in the data (McComas et al. 2014a, Figure 20), on average there is a clear decrease in flux over time (McComas et al. 2014a, Figure 23). These comparisons clearly indicate that the simulations do not reproduce the global evolution of the data, suggesting that the SW dynamic pressure (which is not accurately reproduced in our current simulations) is an important driver of the H ENA flux evolution.

Although IBEX measurements are line integrated, meaning separate sources of flux cannot be easily distinguished from each other from a given direction, Schwadron et al. (2011) were able to approximately separate the ribbon and globally distributed flux for the first year of observations (2009), while Schwadron et al. (2014) performed a similar analysis for the first 5 yr of observations, weight-averaged over time. Although there are currently no data on the separated ribbon and globally distributed fluxes as a function of time over the last 4 yr, it is beneficial for this paper to show the percent change in simulated flux for the IHS and OHS (ribbon) separately. While we cannot show the change in flux every year, for every energy and source, in Figure 13 we show the percent change in flux in 2013.5 compared to 2009.5 for the separate sources (IHS flux in the left column, and OHS flux in the right column), at all passbands. In Table 3, we also show percent changes taken from a few directions in Figure 13, to aid the reader.

Zoom In Zoom Out Reset image size

The IHS flux both increased and decreased across the maps due to the dynamic SW conditions. At 0.71 keV, the flux fluctuated by a few percent near the nose, decreased by ∼20% approximately 60° away from the nose in the ecliptic plane, decreased in the poles (∼5%–30%), and slightly increased in the tail. The ribbon flux, on the other hand, uniformly decreased at mid-to high latitudes near the front of the heliosphere (∼50%), with less decrease at the poles and ecliptic plane, and little change from the heliotail. The IHS flux at 1.11 keV is similar to 0.71 keV, except with a slightly larger increase at the nose. The ribbon flux decreases uniformly at high latitudes, similar to 0.71 keV, although at a greater rate. The changes in the ribbon flux reflect the source of ∼1 keV, neutralized, supersonic SW that decreases at high latitudes as the last solar cycle approached the recent solar minimum. The changes in IHS flux, however, are more turbulent, reflecting the fluctuations in speed and density of the SW that propagate through the IHS, as well as variations in the thickness of the IHS as a function of space and time.

An interesting feature of the time-dependent SW flow is visible in the IHS flux, particularly for 0.71, 1.11, and 1.74 keV. Although difficult to see in the absolute flux (Figure 5), the variation in flux in Figure 13 shows "rings" of decreased intensity that are symmetric about the nose of the heliosphere. At this time, there is an increase in flux from the nose due to fluctuations in the SW plasma, compared to the flanks and high latitudes. Fluctuations in the IHS plasma, due to either the 11 yr solar cycle or short-term variations, interact the strongest with the nose of the heliosphere, i.e., the closest position from the Sun, and waves propagating through the IHS and reflecting from the HP are the strongest near the nose. This creates higher intensities near the nose and later diminishes and spreads to the flanks and higher latitudes. Also, as the fast–slow SW transition region reaches solar minimum and moves to higher latitudes, it creates these horizontal streaks that appear to "complete" the symmetry.

At 1.74 keV, we see a larger decrease in IHS flux from the poles (∼20%–40%) compared to 1.11 keV, with a similar increase near the nose (∼12%), and a similar decrease in the flanks (∼10%). Again, there is little change in the heliotail; however, there is a slight decrease in flux from the tail near the ecliptic plane (a few percent). At this energy, the ribbon flux has increased at low to mid-latitudes (up to ∼90%), with a slight decrease near the poles (∼10%), and less change near the ecliptic plane. These features are similarly seen at 2.73 keV. The ribbon flux increased dramatically at this energy, which was also evident in Figure 11, with more than double the flux in some areas compared to 2009.5. Changes in the IHS flux at this energy are similar to 1.74 keV, although with a smaller decrease in flux from the flanks.

At 4.29 keV, we can more easily distinguish the high-latitude changes in IHS and OHS (ribbon) flux from Figure 12. The IHS flux is similar to 2.73 keV, although with stronger signatures of the fast–slow SW transition region and less change near the nose. The ribbon flux also more than doubles at high latitudes; however, it is slightly drowned out by the IHS flux in Figure 12.

In Figure 15 we show simulated, time-dependent spectra of IHS flux only, between 2005.5 and 2017.5, for all IBEX-Hi passbands. We take spectra from three common directions, i.e., the nose, V1, and north ecliptic pole. It is immediately evident that the IHS flux changes on timescales shorter than the 11 yr solar cycle, dominated by fluctuations in the background SW plasma, e.g., variations in the thickness of the IHS as a function of space and time. There is a consistent oscillation in flux at all energies, with a period between ∼2 and 6 yr that depends on the direction. At high energies, the oscillations in flux from the nose and V1 directions become less noticeable; however, they are stronger from the north ecliptic pole.

Reisenfeld et al. (2012) analyzed the time dependence of IBEX-Hi measurements from the north and south ecliptic poles over the first 2 yr of IBEX operation (2009–2011), showing an energy-dependent drop in flux over time. They showed that little change occurred at 0.71 keV, the most significant drop occurred at ∼1.1 keV, and variations in flux occurred with periods $\leqslant$6 months. Unfortunately, our simulations cannot reproduce variations on such short timescales. Nonetheless, it appears that our simulations indicate a small increase in flux from the north ecliptic pole direction for some energies, for the first 2 yr (2009.5–2011.5), in contrast to the observations. This is likely an indication that our simulations still do not accurately capture the decrease in SW ram pressure during the recent solar minimum. However, careful comparisons with IBEX-Hi data from the polar regions are left for future studies.

The phase of the oscillations in flux slightly depends on energy, but more obviously on the direction of detection. For instance, flux from the nose and V1 directions appears nearly in phase for all energies, although the V1 spectra are slightly behind. Flux at all energies, at low latitudes, oscillates with a period of approximately 2–4 yr, whereas at higher latitudes (i.e., north ecliptic pole), the oscillation period appears longer, although difficult to determine. These results reflect the differences in fluctuations in the background plasma as a function of latitude. Contrary to the observations, there does not appear to be a significant decrease in flux from the IHS from these directions, except for the drop between ∼2010.5 and 2013.5.

In Figure 16 we show simulated, time-dependent spectra of ribbon flux only, between 2005.5 and 2017.5, for all IBEX-Hi passbands. We take spectra from the north, middle, and south portions of the ribbon to study the evolution of the ribbon flux as a function of latitude as well as energy, similar to McComas et al. (2014a). The ribbon flux shows considerable change over time, with little evidence of small-scale fluctuations due to the turbulent solar cycle. This is likely due to the weighted contribution of PUIs averaged over 3$\tau _{{\rm ex}}^{{\rm m}}$ yr (a realistic process), as well as the long integration length through the OHS.

Our results show that the secondary ENA ribbon flux evolves with the solar cycle, with a period out of phase by the times of ENA propagation and PUI-to-ENA charge exchange. This, of course, depends on the direction of detection and the energy of the ENAs. For 0.71 keV, flux from all three directions decreases in intensity from 2009.5 until ∼2013–2015.5. The delay in time between creation of parent ENAs and measurement of secondary ENAs is approximately 6.5–9 yr; therefore, the source of parent ENAs is during the fall of the last solar maximum, which began around ∼2001.5. There is a full recovery of flux at 0.71 keV after ∼2013–2015.5, depending on the direction of detection.

As pointed out by McComas et al. (2014a), the distance to the HP will play a large role in the delay time between measurements from the north and south portions of the secondary ENA ribbon, where we should expect a recovery of flux from the south before the north. For example, toward the southern portion of our simulated ribbon, the HP is relatively close to the Sun (∼140 AU; see Figure 17) due to the interstellar magnetic pressure on the HP. However, we note that the OHS plasma and neutral densities also play a significant role. The densities are larger near the southern portion of the ribbon than the north since it is closer to the nose of the heliosphere, where the plasma experiences the most compression. These larger densities effectively reduce the mean free path of ENAs in the OHS (due to a large plasma density) and reduce the PUI-to-ENA charge-exchange time (due to a large neutral density). Therefore, the source of ENAs from this region is relatively close (i.e., short distance to the HP and short distance of propagation through the OHS), and thus they are detected at an earlier time, with a recovery in flux occurring in 2014. In the middle of the ribbon, the HP is slightly closer (∼130 AU), and the plasma and neutral densities are slightly larger (see Figure 17). While the differences in these parameters compared to the south will complicate the total time delay (e.g., larger plasma density will reduce the ENA mean free path), overall the phase of the flux from the middle of the ribbon is similar to the south, although recovering slightly sooner. Flux from the northern region originates the farthest away (∼170 AU), from a region of the OHS with significantly lower plasma and neutral densities (thus larger ENA mean free path and PUI-to-ENA charge-exchange time). These changes overall cause the flux to be delayed by ∼2 yr compared to the south (see Table 2 for approximate time delays), such that the flux continues to decrease until ∼2015.5 and then recovers. These trends are similarly seen in recent IBEX observations, where flux from the southern part of the IBEX ribbon appears to be leveling off in 2013, while flux from the north continues to decrease (see McComas et al. 2014a, Figure 28).

Zoom In Zoom Out Reset image size

Flux simulated at 1.11 keV produces a similar trend to passband 2, with the latest recovery occurring in the north after 2015, and the earliest recovery after 2013.5 in the south. Again these are largely due to the differences in the distances to the ENA source, as well as the difference in plasma and neutral densities in the OHS. The recovery times are slightly earlier at this energy since the ENAs are faster and thus are detected sooner than at 0.71 keV, particularly in the north. Flux from the middle region appears to be the most stable, with little evidence of the 11 yr solar cycle. This is likely due to the nearly constant source of slow SW near 1 keV in our simulations (although fluctuating about the mean; see Figure 2); and thus slow primary ENAs that fuel the secondary ENA flux at this energy.

It is important to note that higher-energy parent ENAs will travel farther through the OHS before charge exchange occurs, compared to lower-energy ENAs, due to the energy-dependent, charge-exchange cross section (e.g., Möbius et al. 2013). However, even though the mean free path increases with energy, the propagation and charge-exchange times decrease. Also, near the front of the heliosphere the plasma and neutral densities typically increase farther away from the HP, up to the H wall/bow wave (see Figure 17; also see Zank et al. 2013; Heerikhuisen et al. 2014). This, in effect, reduces the distance that high-energy parent ENAs can travel through the OHS before charge exchange. While the time-delay estimates shown in Table 2 may oversimplify the effective mean free path for H ENAs through the OHS, our simulation results suggest that higher-energy secondary ENAs from the ribbon are detected slightly sooner than lower-energy secondary ENAs.

At 1.74 keV, there is a significant change in the evolution of the simulated ribbon flux. At this energy, as expected we are detecting fewer secondary ENAs at lower latitudes and detecting more at higher latitudes. Again, this is due to the existence of faster SW at higher latitudes. The delay in time between parent ENA creation and secondary ENA detection is approximately 5–7 yr; therefore, the fuel for these secondary ENAs was created at a time in the solar cycle beginning in ∼2003.5, during the decline of the last solar maximum. There is no longer a steady decline in flux, like at 0.71 and 1.11 keV, but rather a slight increase, most noticeable in the north and south, and an eventual decline after ∼2013–2014.5. In the middle of the ribbon, there is less flux and only slight variability. While the evolution of flux from low latitudes appears to show slight variation with the short-term solar cycle, flux from higher latitudes is more obviously driven by the transition between fast and slow SW, i.e., the 11 yr solar cycle.

To illustrate this, in Figure 18 we show cross sections of the plasma speed in the Sun's polar plane over the course of the solar cycle. Near solar minimum, the fast SW reaches down to lower latitudes near ±35°. At this time, interstellar H atoms at mid- to high latitudes have a higher probability to survive charge exchange in the slow SW, and instead charge exchange in the fast SW, creating a fast primary ENA. Near the ecliptic plane, however, most interstellar H atoms will, on average, charge exchange in the slow SW (although the SW speed is highly variable at low latitudes; see, e.g., Sokół et al. 2013, which is not reproduced in our simulations). As the solar cycle transitions from minimum to maximum, the slow SW spreads to higher latitudes, reducing the probability of creating fast primary ENAs and thus fast secondary ENAs. This, in effect, creates a maximum in the high-energy flux when the source of secondary ENAs was created near solar minimum, and a minimum in the high-energy flux when the source of secondary ENAs was created near solar maximum, and vice versa for the low-energy flux (also see McComas et al. 2012b, 2014a, for a similar discussion). These maxima and minima are not in phase with the solar cycle, nor with each other, but are delayed by the propagation and charge-exchange delay times as a function of energy and direction. We stress, however, that this process is likely more complicated than we can currently simulate, especially considering our simplified model of the SW.

Zoom In Zoom Out Reset image size

Returning to the flux profiles in Figure 16, at higher energies, we see a global change in the simulated ribbon evolution, where we are now detecting high-energy ENAs that were created in the fast SW. At 2.73 keV, there is even less flux at lower latitudes, with little visible variation over time. At this energy, the peak in the flux in the north and south occurs sooner than the recovery of flux at 0.71 and 1.11 keV. The delay in time between parent ENA creation and secondary ENA detection is approximately 5 yr; therefore, the peak in flux corresponds to the last solar minimum, which began in ∼2007.5. The peak in flux in the south occurs first in 2012.5 and later in the north in ∼2013.5. This energy is particularly important since it is approximately equal to the average energy of the fast SW, where v = 723 km s⁻¹. Therefore, our results show the largest relative change in flux compared to all other passbands, where the transition from slow to fast SW, and thus slow to fast parent ENAs, is reflected strongly in the ribbon flux.

At 4.29 keV, the ribbon flux from all directions is considerably diminished because we are detecting ENAs at speeds near 900 km s⁻¹, which is higher than the average fast SW speed. However, the IBEX-Hi energy response has an FWHM of ∼0.6–0.7E for passbands 2–6, where E is the nominal energy for each passband (Funsten et al. 2009a; McComas et al. 2012b). This allows IBEX-Hi to detect ENAs with energies lower and higher than the nominal energy (in this case 4.29 keV). The energy response functions are applied in our simulations.

While there is little change in the middle of the ribbon at 4.29 keV, there is still an increase in flux from the southern and northern regions. The delay in time between parent ENA creation and secondary ENA detection is approximately 4.5 yr, where the peaks again approximately correspond to the last solar minimum. Surprisingly, the peaks in flux from the north and south appear at nearly the same time in 2012.5. However, this is due to our choice of the northern direction from which to show simulated flux (see Figure 14). As the solar cycle approaches maximum, the majority of fast neutralized SW only exists at higher and higher latitudes, causing the simulated flux from the ribbon at high energies (at a sufficient time later) to appear at higher and higher latitudes. The direction to ${{{\boldsymbol{B}} }_{{\rm LISM}}}\cdot {\boldsymbol{r}} \approx 0$ reaches relatively high latitudes near the north region (up to +60° latitude), but not as "high" in the south region (down to –45°). Therefore, flux from the south region will peak and quickly diminish as the source of fast neutralized SW moves to higher (negative) latitudes, resulting in the flux from the south abruptly decreasing after 2012.5. Flux from the north, however, only gradually decreases after 2012.5 since the high-energy flux from the north is still visible at 1 AU up through +60° latitude. However, our chosen "north" direction is lower than this limit, and the flux from the north actually peaks at higher latitudes approximately a year later in 2013.5.

It is important to reiterate the implications of our assumptions for the SW boundary conditions beyond 2011.5 on our simulated ribbon results. At 1.74 keV, the delay in time between parent ENA creation and secondary ENA detection is ∼6 yr (assuming a distance through the OHS within which the majority of ENA flux has been accumulated, beyond which the time is even further in the past). Therefore, a detection of ENAs at this energy in 2017.5 corresponds to a solar cycle time of ∼2011.5, which is within the range of data from Sokół et al. (2013). For 2.73 keV, the delay is ∼5 yr; thus, simulations of measurements made after 2016.5 can only be approximated by predicting the SW boundary conditions. Therefore, a simulation of flux for passbands 5 and 6 beyond 2016.5 is dependent on assumptions we make for the SW, while flux from passbands 2–4 can be reasonably simulated approximately 6 yr into the future. More recent measurements of the SW in the ecliptic plane could be extracted (see, e.g., McComas et al. 2013); however, we used the data from Sokół et al. (2013) to provide SW speeds and densities as a function of latitude to compute H ENA speeds created in the supersonic SW, which will better represent the energy distribution of the secondary H ENA flux. The IHS flux depends more strongly on SW boundary conditions closer to the time of measurement and therefore cannot be accurately predicted more than a few years ahead of the data, at least for IBEX-Hi energies.