© ESO, 2013

1. Introduction

Launched in late 2008, the Interstellar Boundary EXplorer (IBEX) is designed to detect energetic neutral atoms (ENAs) generated inside, near and beyond the heliospheric boundary layer, i.e. the heliopause (McComas et al. 2004, 2006, 2009b). Because it observes the entire sky in only six months, the mission is not restricted to detecting energy-resolved fluxes of neutral atoms from different directions, but it may, by temporal resolution, also enable studies of the dynamic behaviour of the heliosphere on timescales greater than half a year (McComas et al. 2010, and even smaller timescales towards the poles (see Reisenfeld et al. 2009), offering an enhanced time resolution of solar cycle-induced effects.

The first half-year results were published in late 2009 (McComas et al. 2009a; Schwadron et al. 2009; Funsten et al. 2009; Fuselier et al. 2009), which clearly proved the presence of a completely unexpected phenomenon, a distinct ENA emissive feature at the sky, which has become known as the “IBEX ribbon”. This ribbon, which consists of a special discrete region in the sky observed in essentially all energy bands measured by both IBEX ENA detectors from  ≈ 100 eV − 6 keV, is motivating the development of new ideas and models aiming to explain this unexpected feature (see e.g. McComas et al. 2009a, 2010; Schwadron et al. 2009; Heerikhuisen et al. 2010; Grzedzielski et al. 2010; Chalov et al. 2010; Fahr et al. 2011; Siewert et al. 2012). One of the most challenging properties of the ribbon is that, while generally exhibiting a power-law characteristic, the power index fitting to the flux data does not result turn out to be a special feature in the ribbon (McComas et al. 2009a), which suggests that the physical processes generating the ribbon – besides of their absolute rates – are qualitatively and physically similar to those operating in the remaining regions of the sky, thus resulting in the same qualitative spectral shape, and that the only difference between both regions might be the spectral intensity of ENA sources. Nevertheless, a recent study by Schwadron et al. (2011) in fact identifies three different regions on the skymaps characterised by different spectral power indices, but that do not correlate with the ribbon feature. Until now, this spectral ordering is another feature that needs to be explained by a unified model.

Deriving an improved ENA flux model taking these properties of the data into account is one of the most ambitious research topics in heliospheric astrophysics today (McComas et al. 2010). Possible explanations that are currently under investigation range from unusual magnetic field behaviours in the local interstellar medium (LISM; Prested et al. 2010; Heerikhuisen et al. 2010; Chalov et al. 2010) or properties of the interstellar interface (Grzedzielski et al. 2010), modified thermal pressures (McComas et al. 2009a; Funsten et al. 2009; Krimigis et al. 2009a) to a delicate combination of processes in the local medium, e.g. due to pitch-angle anisotropic pick-up ions (PUIs) in the outer heliosphere or the outer heliosheath. However, without special ad-hoc assumptions concerning the responsible kappa-character of the ion distribution functions, or vanishing pitch-angle scattering, none of these explanations is able to represent high enough spectral ENA-flux intensities at keV energies.

In two recent publications (Fahr et al. 2011; Siewert et al. 2012), we have studied two different sources of ENAs at 1 keV. The resulting ENA fluxes as a function of celestial position were in good agreement, both in intensity and spatial shape, with those flux intensities observed by IBEX. The absolute model fluxes were found to match the IBEX flux intensities within a factor of 2, while the spatial shape was found to produce a narrow ring-like region of enhanced ENA fluxes, reaching from the heliospheric flanks to the poles.

In addition, Schwadron et al. (2011) have presented a detailed study of spectral properties of the IBEX data, finding energy power-laws ( ∝ E ^− k_E) with a variable range of spectral indices, but sorted according to the position on the sky. Specifically, an energy power index of k_E ≃ 1.6 is found close to the ecliptical plane in the upwind (nose) direction; in the tail, the spectral index grows to k_E > 2, while towards the poles, a remarkably harder spectrum is observed (k_E ≃ 1.3). In this study, we examine the spectral properties of ENA fluxes calculated from the two modelled sources by Fahr et al. (2011) and Siewert et al. (2012), especially in the emergence of power laws, their power indices, and possible factors influencing them.

2. ENA sources

2.1. ENA sources in the PUI regime

Perhaps the most important mechanism for generating heliospheric ENAs is resonant charge exchange between energetic ions and cold neutral atoms entering the solar system from the LISM. Since this fundamental process is based on the exchange of a single electron, the resulting ENAs essentially take over kinetic energies and momenta from their parent ions; i.e. in first order approximation, the resulting ENA spectra do reflect those spectra of the source ions before charge exchange. In the following, we study two individual contributions to the energetic protons that represent the source of the ENAs. First, there is the regular PUI source that follows from processing the initial PUI ring distribution (which we call the primary source), and second, there is an additional contribution from ACRs after cooling down from MeV energies to keV energies (Fahr et al. 2009).

The total energy regime covered by both IBEX detectors (IBEX-Lo + IBEX-Hi) ranges from 0.01 to 6 keV (McComas et al. 2009b), while the ribbon has been observed from  ≈ 100 eV to 6 keV. For this reason, we focussed our efforts on this energy range. Considering that Krimigis et al. (2009b) have detected a similar feature at higher energies with INCA-Cassini, which, however, does not perfectly coincide with the placement of the IBEX ribbon, exhibiting a different orientation by  ≈ 30°, we emphasize energies higher than 6 keV as well (i.e. up to 10 − 15 keV). The relevant energy region is dominated by PUIs, which also result from the charge exchange processes. In the rest frame co-convected with the solar wind (with a speed of v_sw ≃ 400 km s^-1), a thermal proton is, actually, replaced by a non-thermal proton with an absolute velocity of v_pui,0 ≃ v_sw, leading to a ring distribution function with a considerably higher thermal energy than the solar wind background. These ions are unstable with respect to driving MHD waves; however, the precise mechanism driving the further development of the ring distribution is still widely discussed (see e.g. Fisk & Gloeckler 2006, 2007; Fahr 2007; Fahr & Fichtner 2013).

In this study, we adopt the model by Fahr (2007), which was also used by Fahr et al. (2011) and Siewert et al. (2012), and which is based on a drift process known as magnetic cooling (Fahr 2007; Fahr & Siewert 2008). This process then results in a velocity power-law (i.e. f ∝ v ^− k_v) with a power index of k_v = 5 of the PUIs at velocities below the solar wind speed (Siewert & Fahr 2008). The physics behind this effect is as follows. Over large regions, especially at large solar distances, the heliosphere can be described using MHD models, and the presence of a frozen-in magnetic field suggests that, in a consistent description, a pressure-anisotropic formulation of the MHD (or kinetic) equations should be used. From anisotropic MHD, it is known that this requirement leads to the conservation of two invariants identified by Chew et al. (1956), one of which is connected with the gyrating ion’s magnetic moment. These invariants can then be transformed into an energy loss term associated with a changing frozen-in magnetic field magnitude, which can therefore be called “magnetic” cooling (contrasting with the more widely used adiabatic cooling). The validity of this assumption is discussed in greater detail at the end of Sect. 2.2.

In addition to this primary PUI contribution, Fahr et al. (2009) studied a secondary mechanism that extends the ion spectrum to much higher energies, up to 100 keV. This region, which is still experimentally inaccessible at the present time, can be filled by injecting ions from an even more energetic regime, e.g. the anomalous cosmic rays (ACRs). As demonstrated by these authors, ACRs generated near the termination shock (TS) and diffusing inwards from the shock, while at the same time being modulated in energy due to adiabatic cooling processes, result in a specific particle energy loss rate. This loss causes the ACRs to lose their diffusive mobility until they are finally again co-convected with the solar wind, analogously to PUIs. Within the framework of conventional ACR modulation theory, this process leads to local sinks in the ACR population, as it is obvious from the divergence of the ACR streaming. Since it is not possible to distinguish between an unusually fast PUI and an unusually slow ACR proton of identical energy, it is necessary to introduce a semi-arbitrary value (adopted by the authors as 10² keV), where the lowest energy population of ACR ions is taken as an additional source for the PUI population. The evolution of this additional source can be derived using the same method sketched above, solving the appropriate Boltzmann-type equation. Then the ACR-induced PUI injection term acting at energies much higher than the conventional PUI injection term, and undergoing the same magnetic cooling processes, results again in a power law, but with a slightly different power index and extending to energies much higher than the PUI threshold of v_sw.

The modified PUI distribution function obtained from these two processes, which we adopt in our study, is thus given in the form (Fahr et al. 2009) (1)where n_p,E is the solar wind proton density at Earth, r and r_E the solar distances of Earth and the point where the distribution function is taken, p_pui is the PUI density and Λ and are normalising factors. Finally, H is the Heavyside step function with H(x > 0) ≡ 1 and H(x < 0) ≡ 0. A typical representation of this distribution function is given by Fig. 1. The main difference to, e.g., Vasyliunas & Siscoe (1976) is found in the power index, where we are obtaining a PUI power index of  − 5 instead of the more commonly quoted value of  − 3/2. This is a major difference that will affect the derived ENA fluxes. As we demonstrate in the following sections, this function is modulated by several additional effects, some of which have, to our best knowledge, not been explicitly studied before.

Fig. 1 Representative plot of the source PUI function (at 70 AU) used in this study, using U = 1 as the initial PUI velocity. Depending on the normalisations of the primary and secondary contributions, a different behaviour is possible. From Fahr et al. (2009).

2.2. Shock-processed pick-up ions

The distribution function presented in Eq. (1) does not consider effects of the solar wind termination shock. However, the shock heating and deceleration of the bulk plasma flow mean that the downstream ions will contribute more strongly to the final ENA fluxes than the upstream ions. Of these two major modifications of the PUI protons caused by a shock, only the second one, the deceleration of the bulk plasma flow, can be described easily. Using the mass flow conservation and introducing the MHD compression ratio s = ρ₂/ρ₁, the corresponding deceleration is described by U₂ = U₁/s. Quantifying the second major shock effect, the kinetic behaviour of the individual particles at the termination shock (i.e. the per-particle “heating” contribution) is far from trivial, considering that MHD does not provide information on individual ions. In the past years, a new model was developed that, based on the conservation of two kinetic invariants, is able to predict the modification of the individual ion velocities and the resulting distribution function (Siewert & Fahr 2008; Fahr & Siewert 2010). In the formulation applied here, we require two shock-related parameters to describe the microphysics of the shock and the downstream distribution function, i.e. the MHD compression ratio at the TS (s = ρ₂/ρ₁) and the magnetic field tilt angle, θ = ∠(n_TS,B), between the shock normal n_TS and the upstream magnetic field vector B. Introducing the functions the velocities transform according to (5)and the downstream distribution function is given by (6)This result is valid as long as the magnetic moment of the individual particles is conserved and no significant energy and momentum exchange with plasma waves occurs during the shock transition.

Obviously, the validity of this approximation strongly depends on the conservation of the magnetic moment. To our best knowledge, this has not yet been discussed in enough detail for non-ideal systems. Fahr & Siewert (2008) discussed this aspect for a magnetic field that changes its magnitude with time when seen in a comoving solar wind bulk frame. It is a very good approximation of the inner heliosphere at solar distances beyond 5 − 10 AU. Studying the TS is complicated by the presence of plasma waves, turbulence, and other kinds of dynamic behaviour (see e.g. Fisk 2005; Jokipii 2008, and accompanying publications), so no clear answer can be given for whether the magnetic moment is conserved at the TS (or other astrophysical shock waves). However, Fahr & Siewert (2008) suggest that a full and consistent electrodynamic treatment is required to explain the behaviour of the magnetic moment, and it is an unfortunate fact that most numerical simulation studies of the TS do not take this into account due to computational limits. Therefore, for the sake of consistency with our adopted (pure MHD) shock description, we assume that the magnetic moment is conserved from the upstream to the downsream side. We are currently preparing a separate publication addressing the question of the conservation of the magnetic moment in greater detail; however, there are hints that the magnetic moment should be conserved as long as it is possible to assume a gyrolimit.

2.3. The frame-of-reference effect

Since ion distributions are described as isotropic in the rest frame convected with the solar wind, the two “direct” shock effects discussed above are not the only effects that modify the final ENA flux spectra. There is also an “indirect” contribution from a combination of shock deceleration and a required conversion between the (moving) solar wind frame and the observer’s rest frame (i.e. the IBEX frame that is, to first order, the solar frame, i.e. non-moving). The relation between the velocity of an individual ion/ENA in the solar wind frame (v_ena), the velocity of the same particle in the observer’s rest frame (v_obs), and the radial solar wind speed (U_r) is given by (see also Fig. 2) (7)In the most general case, this equation is a vector-valued expression, with v_obs oriented towards the (arbitrary) position of the detector. In the case of IBEX, this direction is as a good approximation parallel to the radial direction in heliocentric spherical coordinates, i.e. (8)where t indicates those components perpendicular to the radial direction. Now, the detector can, along its line of sight, only observe those transcharged ions whose velocity vector is oriented perfectly to the observer, i.e. those ENAs with v_obs,t = 0 (as v_obs is the velocity vector in the observer’s rest frame). To make this vector component v_obs,t vanish, we therefore need (9)Obviously, as long as the subsonic solar wind is flowing radially outwards (U_t ≃ 0), this equation therefore simplifies to the scalar expression (10)For our further calculations we now need an adequate model of the heliosheath plasma flow downstream of the termination shock. Unfortunately, however, there is a complete lack of an apropriate flow model from MHD-simulation codes. None of the existing simulation models is able to represent flow features that have been found by Voyager-1 and -2, such as an absence of tangential bulk flow components at decreasing radial bulk components (see Krimigis et al. 2011), or a constant radial flow component as seen with Voyager-2 (see Richardson 2012). This is why we decided here to use a global flow model, based on analytical representation of a 3D-flow vector at every place in the heliosphere.

Fig. 2 Schematic representation: vectors used in the velocity transformation, not to scale.

The model for the plasma flow in the inner heliosheath adopted by us in this study (Fahr et al. 1988) is characterised by the property that the solar wind outflow is radial except for a small region immediately before the heliopause, and therefore, Eq. (10) is a good approximation of the vector-valued expression. If there were observations hinting that the region characterised by tangential plasma flow components may be larger than expected in our basic hydrodynamic model (Fahr et al. 1988), we estimate what kind of further corrections one may expect from this tangential flow. Defining θ as the angle between v_ena and v_ena,r, it follows from Eq. (9) that (11)and the radial velocity component in the direction of the observer can be written as (12)Now, the downstream solar wind speed is approximately U_t ≲ U ≃ 100 km s^-1, while the speed of 1 keV ions is given by v_1   keV ≃ 450 km s^-1, and due to the frame-of-reference transformation, this value must be even higher to produce observed 1 keV ENAs. Therefore, the correction factor in Eq. (12) is (13)so the additional correction can be expected to be small for 1 keV ENAs, and is even smaller for higher values. Only for lower energies is it possible that this contribution becomes more pronounced, and even then, the overall contribution to the line-of-sight integrated integral remains small for our adopted flow profile, which is why we keep using Eq. (10) in this study.

Of course, this approximation no longer holds for different flow line profiles. If, for example, one modelled the extended region of a tangential flow in the heliosheath observed by Krimigis et al. (2011), then this contribution would become more pronounced at low energies owinf to the extended contribution to the line-of-sight integral. Plus, for a different ENA detector, e.g. the INCA/CASSINI mission that is bound to Saturn, the approximation that the detector is co-aligned with the Sun may be no longer applicable. Due to this conversion, it is not possible to produce ENAs originating in normal PUIs (which are injected at  − U and then cooled down) in the inner heliosphere on the upstream side of the TS, where the solar wind outflow is radial, since v_ena,max = v_pui,max = U_r = U.

Some conditions do, however, help overcome this effect on the downstream side of the shock: the increase in v_ena and the decrease in U_r, both of which are realised at the TS. First, it is known from the theory of shock waves (and verified by experimental observations, see e.g. Richardson et al. 2008) that the downstream solar wind is slower than the upstream solar wind, resulting in U_r,2 < U_r,1 (or, respectively, s > 1), which allows observing ENAs from normal PUIs after their convection across the shock. In addition to this, v_ena is larger on the downstream side due to thermal heating, leading to the other condition, and further enhancing the range of observable ENAs. Taken together, these two effects may result in a broad ENA spectrum above o(1 keV) in the observers rest frame.

However, this frame-of-reference effect may also result in deviations from the original power laws of the source population that describe the ENAs in the solar wind frame. Converting a pure power law into the observers frame results in (14)i.e. only for v_obs ≫ U_r a power law in the solar wind frame convert into an approximate power law in the observer’s frame, while for an ENA speed comparable to the solar wind speed, observable derivations from a pure power law must be expected. The additional power of 2 appears thanks to the transformation of the differential phase space volume (instead of “just” the distribution function), which will become more important in the next section.

Fig. 3 Correction factor d(r) (see Eq. (15)) to a pure power law due to the frame-of-reference effect, assuming a radial solar wind speed of Ur = 100 km s-1.

To estimate the importance of this effect, we present the deviation in terms of the factor (15)in Fig. 3 for various velocity power indices. This figure demonstrates that, depending on the power index, the differences may grow significantly large, with slow velocities being visibly underrepresented. If we adopt a frame-of-reference correction of less than 10% deviation (i.e. d(v_c) ≥ 0.9) as an indication that the source velocity power index is reproduced, we see that v_c ≳ 10U_r is required to see the original source spectrum. Considering that 1 keV ENAs are only slightly faster than the solar wind, we must expect significant deviations from the source velocity power index. As an additional effect, this behaviour allows us in principle to estimate the power-law characteristics of the proton distribution function undergoing charge exchange processes.

3. The ENA flux at arbitrary energies

Since IBEX observes ENAs originating along a line-of-sight in the inner heliosheath, fluxes observed close to Earth are inevitably line-of-sight integrated quantities, i.e. one cannot directly infer local details of the structure of the inner heliosphere. On the other hand, neutral atoms possess much longer mean free paths than ions of comparable energies, essentially only undergoing rare reionisation during their propagation. Fahr et al. (2007) and also Bzowski (2008) demonstrate that the mean free paths of ENAs in the energy ranges of interest are large enough to ignore such losses in first order approximations (but see Sect. 3.5 for a discussion on this).

The ENA flux produced along a line of sight in the absence of reionisation is defined by (Lee et al. 2009; Siewert et al. 2012) (16)where Π(v_obs) is the ENA production rate as a function of the observed velocity v_obs. This production rate, parameterised by the ENA velocity and energy in the wind’s frame, can be written in the form (17)where σ_ex is the resonant charge exchange crosss section, j_p is the observer-frame proton flux, parameterised by the energy in the solar wind frame that is converted into ENAs, and n_H is the local interstellar hydrogen density. In this equation, we have selected an energy and velocity subscript ENA, indicating that this quantity is taken in the reference frame convected with the solar wind, which is related to the detector-relevant velocity v_obs by Eq. (7). This transformation, which has been discussed extensively in the previous sections, in general, can only be treated numerically. Finally, (α,δ) is the celestial line-of-sight orientation, which – in our model – influences the source proton flows exclusively; the hydrogen density is actually constant.

In this study, we are working exclusively with proton fluxes as a function of the ENA energy (i.e. j_p(E_ENA)), since it simplifies many of the following calculations, in addition to being more compatible with the IBEX data, which is denoted as a flux spectrum as well. In addition, we now consider that there are two different contributions to the source ion distribution function f, the primary PUIs (subsequently called PUIs), and secondary ions injected from ACRs cooling down (subsequently denoted as PUIs^∗) (see also Sect. 2.1). Using this notation, the total ENA fluxes following from these two sources can be written in the form (18)

To simplify the following calculations, we collect those factors that are constant and can be pulled before the integral, (19a)and (19b)Here, j₀ and σ₀ are constant normalisation factors, while and absorb the energy-dependent contributions. We derive the specific form of these factors in the following sections.

Before we do this, we still need to derive one additional equation related to the ENA fluxes. In the following sections, we derive the source proton fluxes based on distribution functions, so we need the standard relation between the differential flux and the distribution function, (20)The earlier study of ACR-induced energetic PUIs (Fahr et al. 2011) was restricted to calculating 1 keV-ENA fluxes with normalisation to Voyager-2 ACR measurements presented by Decker et al. (2008). Now, to calculate spectral ENA fluxes from the downstream side over a broader range of energies, we need to derive a relation connecting upstream and downstream proton fluxes at arbitrary energies using the kinetic model introduced in Sect. 2.2. On the downstream side of the shock, the proton flux transforms according to (see Eqs. (5) and (6)) (21)For a power law, i.e. f₁(E) ∝ v ^− k ∝ E ^− k/2, the factor C^-1 in the energy argument can be pulled out f₁, and then this expression reduces further to (22)For the secondary, ACR component (k = 4), this transforms into (23)while the primary, regular PUI component results in (24)

3.1. The resonant charge exchange cross section

Since the charge exchange cross section is a function of E_ENA, it may in principle modify the functional form of the integral. In earlier studies, only 1 keV ENAs (in the observers frame) were considered, and the corresponding cross section reduced to a constant value for this energy (Fahr et al. 2011; Siewert et al. 2012). Now, however, we need to parameterise the cross section as a function of the relative kinetic energy E_p ≃ E_ENA.

Various analytical approximations to the charge exchange cross section can be found in the literature (for an overview, see Fahr et al. 2007). In this study, we use a simplification of the underlying cross section that is valid in the energy range of (o(1 − 10      keV)). This was derived by McNutt et al. (1999), who found (25)where g is the value of the relative speed in km s^-1 and a₀ is the Bohr Radius. Transforming the factor g = v/1 km s^-1 into an energy ratio using 1 keV as a reference value, we obtain the parameterisation (26)with σ₀ = 2.44    ×    10^-15 cm ². This σ₀ is the same value as appears in Eq. (19), and the energy-dependent term appearing below the integral is (27)This approximation is not the most recent one (see Lindsay & Stebbings 2005, for a more up-to-date parameterisation). Using the more recent parameterisation, we obtain σ_this(1      keV) = 1.5σ_LS(1      keV) and σ_this(10      keV) = 2.2σ_LS. To estimate the impact of this correction, we first assume that the cross section is parameterised by (σ_LS(E) ≃ σ_this(E)E^k_σ). Then, using the values given above, we can estimate the correction k_σ to the resulting total power index, resulting in (28)We briefly discuss the uncertainties introduced due to this difference in Sect. 3.5.

3.2. The shape of the heliospheric boundary layer

The line-of-sight integration boundaries r₀ = r_s and r₁ = r_h are determined by the geometry of the heliospheric boundary layers, the termination shock (TS), and the heliopause (HP). The TS distance can be easily derived from a heliospheric TS model, which also allows us to derive the magnetic field tilt angle θ = ∠(B,n) required to describe the reaction of the proton distribution function to the TS. Since the heliosphere is a complicated dynamic system, the TS geometry and the resulting angles θ_Bn are not reliably well predicted from theoretical and numerical considerations alone. Due to a lack of observational data, we are restricted to magneto-kinetic/-hydrodynamic models based on varying degree of physical sophistications. One of the most recent theoretical models for the geometry of the TS has been derived by Borrmann & Fichtner (2005); in this study, we do, however, make use of the analytic approximation by Fahr et al. (2008), who describe the TS as an ellipsoid with two identical short (104.7 AU) and one long (133.5 AU) semi-major axes, which best fit the shape of the termination shock obtained by Borrmann & Fichtner (2005). This model allows an analytical representation of both the TS normal and distance, considerably simplifying further calculations.

To explain the lack of unfolding of the ACR spectrum as the Voyager spacecraft approached and crossed the termination shock, McComas & Schwadron (2006) propose that a blunt termination shock can naturally explain these observations. In particular, this geometry, characterised by the nose being significantly closer to the Sun than other parts of the termination shock, produces multiple shock crossing points of the swept-out Archimedian magnetic field, which accelerates ACRs up to their highest energies along the flanks and tail, where the field lines have been connected to the shock for roughly a year. Subsequent calculations by Schwadron et al. (2008) quantify these effects and show that the blunt shock geometry can in fact reproduce the observations. In this study we first use the model of Borrmann & Fichtner (2005) with a three-axis ellipsoidal fit to the shock surface using an elongated ellipsoid with axes a > b,c. Then, to assess the effects of the blunt shock geometry, we also examine a flattened (i.e. prolate) ellipsoid with a < b,c in Sect. 4.2.3.

The precise structure of the outer heliosphere, i.e. beyond the TS, is even less understood than the geometry of the termination shock, since the most distant missions to probe this region, i.e. the Voyagers, have not yet reached the heliopause. There are hints in the latest findings by Voyager-1 (Krimigis et al. 2011) that the heliopause might be closer to us than expected from earlier theoretical modellings (130 AU, in contrast to the 150 AU used here). However, it is also possible that the heliopause is not a narrow surface (i.e. a thin layer), but made by an extended pre- and post-HP region, as predicted theoretically by Fahr & Neutsch (1983), among others. Until existing models have been modified for recent Voyager measurements, we ignore those recent indications.

The position of the heliopause can be estimated based on existing inner heliosheath models, such as the streamline model developed by Fahr et al. (1988). Using this model, the location of the heliopause here can be calculated as (29)or (30)for a solar system moving at constant velocity through the interstellar medium, without any acceleration or deceleration. Here, L is the stand-off distance of the heliopause nose, and the angle Ψ of the line-of-sight with respect to the symmetry axis is defined by . (For a graphical representation and definition of the coordinate system, see Fig. 4.)

Fig. 4 Schematic representation: the coordinate system used in the potential flow model for the plasma stream lines, not to scale.

This model describes an open heliotail, where the outer distances in the downwind tail approach infinity, which might result in unphysical results since the energetic ENAs only have a finite lifetime before reionisation, which is not considered in this study. On the other hand, the TS in this region is significantly more distant than in the nose region, and the PUI density will be accordingly smaller there as well (due to dropping off according to r^-2). Due to these effects, the ENA fluxes in the tail region are expected to be significantly smaller, and we continue using this model even in the heliotail, cutting off the line-of-sight integral at approximately the ENA extinction length, which we estimate to be 300 AU, i.e. twice the nose distance to the HP. Interestingly, significant modifications due to the increased distances have been hinted by Schwadron et al. (2011), who found that the ENA spectra seem to possess systematically softer spectra in the heliotail. This effect should be reproducible based on transport effects and will be the focus of a future study. The present study is based on the approximation that absorption processes may be ignored (see also Fahr et al. 2007).

We also need to point out that our representation of the heliospheric boundary layers (i.e. the TS and the HP geometry) are based on hydrodynamic models. An improved, magnetohdyrodynamic description, including more realistic representations due to tilted LISM fields, may certainly change the outcome, including (but not limited to) breaking the axisymmetry assumed in the present study. This may increase the overall match with IBEX data, and will be investigated in a future study.

3.3. The radial component of the solar wind speed in the inner heliosheath

Next, we need to consider that the subsonic solar wind in the inner heliosheath is deflected from a perfect radial direction when approaching the heliopause, and U_r is therefore not constant here, but a function of solar distance. Therefore, an individual line-of-sight will inevitably cross many different, nonradial plasma flow lines downstream of the TS (see Fig. 4), and U_r(r,α,δ) needs to be parameterised accordingly. The factor in front of the integrals in Eq. (18) is defined by all contributions that remain constant along any individual streamline due to the incompressibility of the downstream medium, while the modified radial speed will result in additional derivations due to the frame-of-reference effect. Depending on the specific, model-dependent, streamline velocity profile in the inner heliosheath, this may, again, result in a different modulation of the final observed ENA spectrum at Earth. Such a possible modulation might be an explanation of the direction-dependent ordering in the hardness of ENA flux spectra pointed out by Schwadron et al. (2011) and others.

To derive an expression for the solar wind speed in the direction of a line of sight that is compatible with the model for the geometry of the inner heliosheath, we adopt the model by Fahr & Fichtner (1991), who described the inner heliosheath flow using a velocity potential field, which allows deriving the radial velocity component by (31)Here, L ≃ 150 AU is the heliopause stand-off distance (i.e. the distance between the sun and the heliopause in the nose direction), the normalised axisymmetrical coordinates are denoted by (for the distance perpendicular to the symmetry axis) and (for the distance along the symmetry axis), and Ψ is the angle introduced in the previous section (see Fig. 4). The parameters A = 1 and B = 0.5 are model parameters that take these given values for a solar system moving freely, without any acceleration or deceleration, through an unchanging LISM. The factor Γ must be selected in a way that, at the TS, the velocity is exactly the solar wind velocity downstream of the shock, which then results in the normalisation (Siewert et al. 2012) (32)Using this parameterisation, a power law in the solar wind frame is expressed in the Sun’s frame by (33)where we have introduced the abbreviations (34)and (35)This allows the basic energy-dependent part of the distribution functions to be parameterised (i.e. Eq. (1)), leaving only the normalisation factors as unknown parameters.

Since we are working with energies instead of velocities, we now derive an equivalent formulation for energy power laws. Beginning with (36)with (37)one immediately sees that the correction terms on the right-handed side are formally independent of specific model details on the inner heliosheath streamline profile. This allows us to study the frame-of-reference effect just by studying E_corr, and the initial study performed in in Sect. 2.3 for velocities can be directly applied to the energy formulation as well.

Interestingly, an analytical solution of the line-of-sight integral (see Eqs. (16) and (20)) for power-law distribution functions in the solar wind frame exists. However, the resulting expression involves non-elementary (hypergeometric) functions, which can be difficult to evaluate numerically; for this reason, we will perform the line-of-sight integrals numerically.

3.4. The source proton flux distributions

The final remaining parameters that needs to be described in the line-of-sight integral are the normalisation factors , after splitting them from the energy-dependent power law. Both of the source functions considered in this study have already been studied earlier in the context of ENA production (Fahr et al. 2011; Siewert et al. 2012); in this subsection, we therefore summarise these earlier results and derive the correct formulation in the adopted unit system (i.e. cgs units, with the exception of energies given in keV, and power laws using energy units instead of velocities as well).

3.4.1. The PUI distribution function

First, we derive the formulation of the regular PUIs, which are characterised by f(v) = f₀(v/v₀)^-5. Using this velocity-dependent parameterisation, we found that (Siewert et al. 2012, Eq. (19)), (38)on the downstream side of the shock (using v₀ = U₁). Here, s is the compression ratio, n_pui,1 is the upstream PUI number density, U₁ the upstream solar wind speed at the shock, and C a factor related to thermal heating (see Eq. (4)).

Inserting Eq. (38) into (20) then results in (39)Using the notation introduced in Eq. 19a, we obtain a constant factor of (40)and an energy-dependent part (41)The parameter n_pui,1, i.e. the PUI density, requires some additional work. We begin by introducing the PUI abundance ξ = n_pui/n_p, which as derived by Fahr & Ruciński (1999), is independent on the solar wind velocity U, (42)Here, β is the angle between the radius vector and the upwind direction with the relation cosβ = cosδcosα, the upwind direction given by {α = 0;δ = 0}. The angles α and δ are the ecliptic longitude and latitude. This allows us to derive a simple representation of the PUI density, (43)where n_p,E ≃ 4 cm^-3 is the solar wind plasma density near Earth. Depending on solar activity, local events such as coronal mass ejections, and even on the reference being quoted, this value can vary significantly; however, scaling this value by a factor of x results in propagating this factor to the final ENA spectral intensity, while the spectral properties remain unchanged. Since the PUI abundance at the shock position essentially turns out to be independent of the direction, we may simply use ξ(α,δ) ≃ ξ(0,0) ≃ 0.15, and Eq. (40) transforms into (44)Here, it must be taken into account that, under solar minimum conditions, the solar wind speed U₁ has been observationally found to be a function of the heliospheric latitude (McComas et al. 2000).

3.4.2. The ACR/energetic PUI distribution function

The ACR-induced component (i.e. the parameters and j^′ ∗ (E)) has not been studied in a fully analytical form yet; instead, Fahr et al. (2011) based their study on an extrapolation of the (more realistic) flux values observed by the Voyagers, obtaining an ACR-induced flux of Φ(E_obs = 1.12      keV) = 3.92    [cm^-2   s^-1      keV^-1   ster^-1]  in the nose region ((α,δ) ≃ (0,0)). Following this, it was possible to extrapolate the PUI^∗s to different celestial positions due to the form of the ACR injection function, (45)where θ = ∠(B,n). The distribution function was then found to have the functional form (46)i.e. the main difference between the reference value and other positions on the sky is defined by a ratio of injection functions ϵ.

Fig. 5 Schematic representation: The behaviour of streamlines characterised by the solar wind speed for solar minimum conditions, not to scale. Obviously, a modification can be expected in the heliotail region, where lines-of-sight cross deflected streamlines.

Considering that the ENA flux is, to first order, directly proportional to the distribution function, one can trivially derive that the differential proton flux due to energetic PUIs is given by (47)Adopting a reference energy of 1 keV, and making use of the v^-4 = E^-2 spectral shape of this ion component (see Eq. (1)), one easily obtains the final form (48)This expression now needs to be divided into a constant term  and an energy-dependent term j^′ ∗ (E), according to the parameterisation introduced in Eq. (19b). This is slightly complicated by the streamline behaviour in the inner heliosheath, where streamlines are deflected from a perfect radial outflow (see also Fig. 4). As a result, any line-of-sight may inevitably cross streamlines of strongly varying injection efficiency, i.e. ϵ(α,δ) → ϵ(α,δ,r) (see also Fig. 5), and ϵ(α,δ) can in principle not be absorbed in the constant factor. However, for the specific inner heliosheath model adopted here (see Sect. 3.3), it turns out that the plasma outflow is essentially radial until very close to the heliopause, i.e. the varying-speed effect can be ignored in good approximation, and we can represent the ACR-induced ENA flux component using the parameters (49)and (50)The final missing parameter in these equations is the reference flux , which can be derived following the formalism used by Fahr et al. (2011). We begin with the mostly constant flux of low-level energetic ions observed by the Voyagers between 40 and 53 keV (Decker et al. 2008) far ahead of the shock, which was found to be (51)Using Eqs. (49) and (50), it is possible to derive the reference value j₀ at 1 keV, (52)This value is slightly different from the adopted value by Fahr et al. (2011), who were using a slightly higher reference energy. The correction due to the Voyager position that does not perfectly match the nose position is negligible (Fahr et al. 2011).

This value is based on the upstream solar wind flow and needs to be adapted for the downstream side, which results in the most significant contribution to the ENA flux. According to our adopted kinetic model (Eq. (21)), this results in the additional factor (53)Taking the power-law characteristics into account, the downstream velocity conversion introduces an additional factor of C², and we obtain the final form (see Eq. (23)) (54)In the following, we focus on this downstream (i.e. inner heliosheath) component.

3.5. Survival probabilities and other propagation effects

There are more aspects that have to be considered to arrive at a largely consistent description of ENA flows. Perhaps the most important aspect that we have not discussed yet is the survival probability ξ, which considers that ENAs generated at distances of more than 100 AU may undergo reionisation processes during their propagation towards the detector (i.e. f_ENA,obs = ξf_ENA,source).

Fig. 6 ENA flow skymaps at Eobs = 1 keV for both source components. Highlighted points are used as reference positions for the spectral study in Sect. 4.2.

A consistent description of the ENA propagation requires detailed knowledge of the solar wind mass density in the inner heliosphere. However, most existing (analytical and numerical) modelling attempts are based on global descriptions and do not consider local effects, such as CMEs propagating towards the outer heliosphere. However, a detailed description of these aspects would easily be worth multiple separate studies, which is why we do not go into detail here.

Early order-of-magnitude estimates have suggested that about 90% of ENAs generated in the outer heliosphere may reach an earthbound observer (Fahr et al. 2007; Lee et al. 2009). More recent studies (see e.g. Bzowski 2008; Schwadron et al. 2011) demonstrate that this effect must be considered in greater detail, with survival probabilities reaching values of about 30% at 0.1 keV, which should definitely affect the lower end of our spectrum. To estimate the impact of this correction, we first assume that the survival probability is parameterised by (ξ(E) ≃ ξ₀E^k_ξ). Then, using the values ξ(1   keV) ≃ 0.3 and ξ(6   keV) ≃ 0.8, we can estimate the correction k_ξ to the resulting total power index (as derived in Sect. 4) (55)suggesting that the spectral power index of the observed ENA spectrum should harden by about k_ξ = 0.24. Considering that the majority of reionisation events occur within 10 AU from the sun, it can be expected that the correction factor due to this effect can be pulled before the integral, and all our results should react in the same way.

However, a similar correction was found during the discussion of the charge exchange cross section (Sect. 3.1), where it was found that, when using a different parameterisation of the cross section, the spectrum was found to soften by k_LS =  −0.17. Combining both effects, one easily sees that the modification of the spectral shape almost perfectly cancels out, leaving only an error of Δk_tot = o(0.1). Corrections of this order of magnitude can be associated to many (theoretical, numerical, or even instrumental) uncertainties, which is why we do not consider these aspects in more detail here.

Other than this, these corrections will influence the absolute flux magnitudes, which however are not the focus of this paper. A first estimate suggests that, after including these effects, the fluxes derived here overestimate more detailed results by a factor of 2 − 3, which is negligible for order-of-magnitude estimates. A more detailed study of propagation effects (e.g. extinction, but also interaction with local density fluctuations in the solar wind) is currently under preparation, but also see the conclusions.

4. ENA spectra and skymaps for both source processes combined

4.1. Skymaps

Following the introduction of our models, we now evaluate actual ENA flux spectra. We begin with a merged skymap, following the initial studies by Fahr et al. (2011) and Siewert et al. (2012), that demonstrates the general behaviour of the resulting ENA fluxes over the entire sky. There are a few differences between the model used in the two initial studies and the one used here, mainly the approach to the cross section (which was assumed to be constant earlier), and some simplifications which were not perfectly compatible with each other between both previous studies. Unless noted otherwise, and not to overly complicate the following study, we are using solar maximum conditions.

In Fig. 6, we present ENA fluxes at 1 keV, based on the models introduced above, separately for the two individual contributions and the merged view, where both components are present at the same time. These figures easily demonstrate various properties of our model. First, one immediately sees a symmetric ring-like structure encompassing the entire flank of the heliosphere (i.e. Ψ ≃ 90°, see Fig. 4 and Eq. (30)). This dominating feature is a direct result of the PUI component, which was found to be approximately proportional to the line-of-sight length in the inner heliosheath, and also proportional to the inverse square of the termination shock distance r_TS (see Siewert et al. 2012). There two concurring effects both start to become significant in the flanks of the heliosphere, with a narrow region where the first effect (which enhances ENA fluxes) dominates the second effect (which reduces ENA fluxes, eventually overcompensating for the initial gain towards the tail region).

The second contribution to the ENA fluxes, the ACR-induced component, does not seem to significantly modify the main component. Focussing exclusively on two isolated regions close to the ecliptic plane, the maximum ENA fluxes obtained from the secondary ACR-induced component do possess an interesting overlap with the primary ring feature; however, the absolute values (at 1 keV) are too low to result in a significant enhancement. It should be noted, though, that the local enhancement inferred by this secondary contribution possesses some similarity to the IBEX data, where it was found that two regions in the ribbon-like shape, most often called “knots” in the literature, seem to exhibit a local enhancement when compared to the remaining sections of the ribbon. In light of these experimental observations, the secondary ACR component might in fact explain this behaviour.

4.2. Flux energy spectra

After deriving all-sky maps of the ENA fluxes resulting from this model, we now present ENA energy flux spectra, i.e. fluxes as a function of the ENA energy observed at Earth. Using the general properties of the skymaps presented in the previous section, we select a representative sample of points on the sky for which we derive detailed spectra. All of the following points are presented, graphically, in Figs. 6a − c.

First, we select several points located in the ribbon, i.e. our main area of interest. The first reference point has been selected at (α,δ) = ( − 90°,0°), located directly in the ribbon, reasonably close to an ACR-enhanced region, while the second reference point uses (α,δ) = ( − 115°,0°), representing a region close to the maximum of ACR-induced ENAs. Finally, we select the north pole (α,δ) = (0°,90°) as one reference point far away from ACR-modified regions, allowing at the same time studying our model in a region where the experimental time resolution is greatly enhanced (Reisenfeld et al. 2009).

In addition to those points thath are correlated to the ring (i.e. the region of local maximum), we also select two more reference points in the background-dominated regions. First, we select the heliospheric nose (α,δ) = (0°,0°), representing ENA fluxes from a region where no significant modifications are expected because the TS and HP distances mostly remain unchanged, and second, we select a point in the heliospheric tail. For this point, we do not select a point reasonably close to the tail position since our theoretical models do not represent this section with sufficient accuracy (see also Sect. 3.2); instead, we select the point (α,δ) = ( − 150°,0°), which is associated with a modelled heliopause distance of approximately 300 AU, i.e. precisely the point where we implemented a cutoff as a first order approximation.

Fig. 7 ENA spectral flows as a function of energy at selected positions in the sky, for PUI-induced ENAs (blue line), ACR-ionduced ENAs (green line), and the sum of both (red line).

ENA spectra for both the individual and the total flux spectra are presented in Fig. 7, covering the energy range of 0.1 − 10 keV, i.e. two decades covering the most significant regions of the IBEX range. One immediately sees, again, that the ENA spectral fluxes are dominated by the PUI component. However, Fig. 7c suggests that ENA spectra connected with ACR-injected ions may become significant at higher energies, resulting in a modification of the spectral slope, similar to the “ankle” – like behaviour pointed out by Dayeh et al. (2012). The contribution of the ACR-induced component may grow significant at higher energies, contributing about 25% to the total ENA fluxes at 10 keV.

4.2.1. ENA flux power indices

Since the absolute fluxes presented in Fig. 7 cover several orders of magnitute, it is not easy to understand the impact of the various contributions. For this reason, we now focus on a modified representations of the ENA flows. First, we are interested in the effective power indices k_E of the flux spectra ( ∝ E ^− k_E). Multiplying the ENA fluxes by (E/1   keV)^1.5, we are able to convert the ENA fluxes into a more compact form, which is presented in Fig. 8.

The total (PUI + ACR) fluxes presented there easily demonstrate several points that have been predicted by our theory at earlier points in this study. First, one easily sees that the resulting power indices only asymptotically reach a value of k_E =  −1.5, which would provide a good representation of the underlying PUI spectrum. This overall behaviour was predicted by the discussion of the frame-of-reference effect in Sect. 2.3. As a noteworthy side effect, this behaviour proves that there should be no global power index characterising ENA fluxes in this energy region; depending on the accuracy of data, however, this effect might not be observable.

In addition, the Fig. 8 also demonstrates a specific behaviour that is already suggested by Fig. 7c. At high enough energies, the E^-2 behaviour of the ACR-induced ENA component starts to emerge, resulting in an asymptotic behaviour that overrides the PUI-induced E^-2.5 behaviour. This results in a hardening of the spectrum, and might provide an explanation for the “ankle”-like behaviour recently discovered in the IBEX flux data (Dayeh et al. 2012).

The recent study by Schwadron et al. (2011) demonstrated an ordering of power indices by heliospheric direction in the IBEX data. In the nose region, and close to the ecliptic plane, the energy power index was found to be k_E,nose ≃ 1.6 − 1.8, which is slightly softer than our model predictions, which are k_E,nose ≲ 1.5. In the tail, the slopes were found to be much higher, k_E,tail ≃ 2.3, which is completely off our predictions and might be a result of reionisation due to longer line-of-sight distances in this region (which we have not yet included in our model). Finally, the fluxes near the poles were found to be flatter (k_E,poles ≃ 1.4 − 1.5), which was attributed to the fast solar wind at high latitudes. Considering that our model exclusively uses fast solar wind conditions at all latitudes, the match in this region can be expected to be better than close to the ecliptic plane. Comparing the IBEX data (k_E,poles ≃ 1.4 − 1.5) does, in fact, result in a better match with our model.

4.2.2. Model predictions using solar minimum conditions

Aiming to improve our match near the ecliptic plane, we now repeat a selected subset of our study using solar minimum conditions, which naturally introduces an ordering of the solar wind speed with the latitude. In Fig. 9, we present the total ENA flux spectrum for both solar maximum and minimum conditions, which demonstrates several interesting differences for the solar minimum.

Since, under these conditions, the solar wind is slower at low latitudes (by a factor of approximately 2), it must be expected that the frame-of-reference modification of the source power-law is partially suppressed at lower energies, since reducing U_r in Eq. (15) moves the correction factor closer to unity. This effect is in fact visible in the model, resulting in a marginally steeper ENA flux spectrum. However, this steepening is not sufficient to reproduce the observed IBEX ENA spectral index of k_E ≃ 1.6 − 1,8, suggesting that an additional contribution to our model is needed.

4.2.3. The impact of the TS geometry

In this section, we present a brief estimate of how strongly the TS geometry affects the shape of the ribbon. We study two separate, different approximations to the TS, one based on a prolate ellipsoid (in contrast to our adopted triaxial ellipsoid), and the other model uses the analytic fit by McComas & Schwadron (2006).

Since our model depends on the magnetic field tilt angle, which has not yet been evaluated for the model by McComas & Schwadron (2006), we are unable to present a detailed comparison of full skymaps here. One possibility for simulating this behaviour would be to take our ellipsoid heliospheric model and select the half axes in such a way that a < b,c, which would likewise result in a strong broadening of the heliosphere. However, this approach might result in the Sun no longer being located in a focus point, and the heliotail might be modelled incorrectly as well. In addition, the (currently not well understood) microphysics resulting in such a TS surface will, most certainly, affect the flow profile model, along with the HP position. One final complication is the modelling of the heliotail. The ACR data on which the model by McComas & Schwadron (2006) is based is strongly focussed on the upwind region of the heliosphere, i.e. modifications in the flank and tail regions, where our ribbon-like emission feature dominates, would have to be modelled again using (M)HD simulations. In the absence of ACR data near the flanks of the TS, we would be forced to continue using the model on which this study is based.

Fig. 8 Modified ENA flux spextra (j(E)·E1.5) at the reference points as defined in Fig. 7, and represented graphically in Fig. 6.

Fig. 9 Modified ENA flux spextra (j(E)·E1.5) in the ring region (α,δ) = (−90°,0°) for the fast and slow solar wind regions.

On the other hand, it is possible to qualitatively predict the impact of a blunt heliosphere on our results. First, one has to consider that, in our model, approximately the entire nose region of the heliosphere ( − 90° ≤ α ≤ 90°) is, in zeroeth order, characterised by a mostly constant TS distance and that the magnetic field can be approximated as perpendicular. Only in the flanks does this approximation break down strongly enough so that this model’s ring feature emerges. Now, a modification of the shape to a blunter form would immediately result in an inclined magnetic field. However, from our general study, we do know that the enhanced ring feature emerges primarily in regions where the magnetic field-TS tilt angle reaches values of θ_Bn ≲ 70°. Therefore, the behaviour of the upwind hemisphere using a blunt TS can be estimated as follows: if the TS is sufficiently non-perpendicular due to the bluntness, it can be expected that an additional feature emerges in the nose region, possibly filling out the entire ring. If, on the other hand, the TS is still close enough to a perpendicular orientation, there should be no significant modifications. This estimate must be taken with care, though, since we do not know whether the HP position would be modified by the microphysics generating the blunt TS as well.

A more detailed study of different heliosheath geometries is clearly beyond the scope of this paper, and will be the focus of a separate study. However, we can already point out the following geometrical aspects that should influence the shape of the ribbon. First, as pointed out by Siewert et al. (2012), there are two major contributions to the ring feature emerging in our model. The first is the solar wind number density, which in the inner heliosheath drops according to r^-2 and which contributes to the amplitude of the ENA fluxes. After crossing the TS, the solar wind plasma density no longer continues to drop; i.e. the normalisation of the ENA source distribution behaves according to . The TS distance is almost constant in the upwind direction, but grows by a factor of 2 − 3 in the downwind direction, which suppresses downwind ENA flows.

The second contribution is the length of the line-of-sight integral, i.e. the distance between the TS and the HP, or the line-of-sight thickness of the inner heliosheath. Since the ENA fluxes are derived from a line-of-sight integral over just this distance, the number of ENA events observed by the detector is, in a first order approximation, a linear function of this distance. In nearly all heliospheric models, the thickness of the inner HS reaches a minimum in the upwind direction, and a maximum in the downwind direction. Taking both aspects together, one sees that there is one narrow, ring-like region in the flanks where both contributions do not start suppressing the ENA flows for the one or another reason, and where ENA fluxes can be much higher than in the remaining regions.

Using these arguments, it is easy to see why the ENA emission feature predicted by this study is strongly modified by the TS geometry. It is known from recent measurements that the IBEX ribbon approximates our ring feature quite well when the Mollweide projection is centred not on the nose, but on a different symmetry direction; “likely (the) external magnetic field direction” (McComas et al. 2012). Nevertheless, this identification of the symmetry direction with physical properties is speculative, and will be likely so for some time. The most likely physical interpretation of the shifted symmetry direction is that this region is characterised by the greatest external (magnetic+kinetic) pressure, while the geometrical shock configuration resulting under this pressure is still under discussion. Considering the fact that the inflow direction of the LISM and the associated local interstellar magnetic field (LIMF) are assumed to be disaligned in the majority of the recent literature, no clear identification of the ribbon symmetry axis with one characteristic external direction can be expected. Even though basic geometric properties of the boundary layer are not a focus of existing heliospheric models, it should now be understandable how – qualitatively – the LISM and also the LIMF affect the geometry of the heliopause and the symmetry of the resulting ENA flows derived from our model.

4.2.4. The heliotail

Another point that requires careful consideration is the modelling of the heliotail. In our adopted basic pressure-balance model (Fahr et al. 1988), the heliopause is open in the tail region, and the HP distance reaches unphysically high values, which would result in an infinitely high ENA intensity in the tail. This in turn would require careful modelling of transport effects, such as reionisation and maybe secondary re-reneutralisation as well.

A more complicated model for the heliotail developed by Jaeger & Fahr (1998) is based on additional physical effects, such as charge-exchange cooling of the tail plasma which includes a tail collapse which dominates the heliotail. This improved model arrives at a closed heliopause, characterised by another ellipsoid bubble enveloping the TS, with a tail distance of approximately r_HP,tail ≃ 6r_HP,nose ≃ 1200 AU. This distance is large enough that transport effects have to be taken into account. To estimate this effect and its impact on our ring feature, we present a modified skymap with a cutoff value r_HP ≤ 1200 AU in Fig. 10, which clearly demonstrates that the increasingly long distances in the heliotail do introduce a significant complication, that reionisation effects must, in fact, be treated in greater detail than before, and that energy-dependent modifications must be expected.

Fig. 10 Modified ENA flux spectra for a heliotail HP cutoff distance of 1200 AU as a first approximation to the model developed by Jaeger & Fahr (1998).

Our results obtained for a more realistic extinction length of 300 AU remain untouched by this result. The assumption that 1200 AU ENA sources remain fully visible on Earth is highly unphysical; following the study by Bzowski (2008), who studied the survival probability for ENAs generated at 100 AU and find a survival probability between 30% and 70%, centring on 50%, depending on the energies being studied. Therefore, even our conservative extinction length of 300 AU may already overestimate the effect, and the 1200 AU allowed in the model by Jaeger & Fahr (1998) is completely unphysical. Including these transport effects is far from trivial and beyond the scope of this current study. The first experimental studies based exclusively on modifications of the cutoff distance suggest that the overall ring-like shape of the feature remains constant.

5. Conclusions

This study presents differential ENA flux spectra developing out of a resonant charge exchange process with both the primary (PUI) and the secondary, ACR-induced component in the pick-up ion regime (Fahr 2007). Both components have already been studied at a fixed energy of  ≃ 1 keV (Fahr et al. 2011; Siewert et al. 2012), which is now complemented by a spectral energy study.

In this study, we demonstrate various important physical aspects of our model. First, we have seen that the source proton spectra undergoing charge exchange are not linearly converted into ENA spectra, which follows from a frame-of-reference effect. To be precise, any energetic ion (or ENA) distribution function convected outward with a solar wind speed v_sw will be seen by an Earth-bound observer with a reduced velocity of at most v_obs = v_ion − v_sw. For a power-law distribution function, this results in a hardening of the spectrum in regions where v_obs ≃ v_sw.

In addition in our ENA production model, one of our two considered source populations completely dominates the resulting ENA flows. This population is generated by upstream solar wind PUIs being heated and compressed during the passage of the termination shock (TS), resulting in downstream ions at high enough energies that, after neutralisation due to charge exchange processes, can be observed at Earth. The other contribution, namely ENAs generated from ions that have cooled down from the ACR regime, usually are too low in number density to contribute significantly to the system, at least for the presently assumed injection rates. However, one should also keep in mind that ACR injection rates are, essentially, a function of the magnetic field orientation at the TS (Scherer & Fahr 2009), but which also depends on several other unknown properties. In addition, the magnetic field orientation that influences the microphysics of the shock transition at the TS likewise depends on model calculations for the solar wind plasma flow, since the magnetic field tilt angle is a function of TS geometry, and also the frozen-in magnetic field, for which various competing approaches exist as well (see e.g. Parker 1958; Fisk 1996). Considering these uncertainties, it cannot be excluded that the position of the “knee”, i.e. the position where the (harder) ACR-induced ENA component starts to dominate the ENA fluxes, might be moved to energies accessible to IBEX.

Finally, we compared the resulting ENA power spectra with the recent study by Schwadron et al. (2011). This study points out that the IBEX data exhibits an ordering of the inferred ENA power-law indices, which can be roughly divided in three regions. The heliospheric nose is characterised by ENAs with a power index of k_E ≃ 1.6 − 1.8, the heliotail possess much softer power tails with a value of k_E ≃ 2.3, and finally, high latitudes, i.e. the region near the heliospheric poles, is characterised by the hardest spectra, with k_E ≃ 1.4 − 1.5. Of these three regions, our model is able to reproduce, within reasonable margins, the high latitude power indices. The low-latitude nose behaviour can be reproduced qualitatively by taking into account the ordering of solar wind speeds observed by McComas et al. (2000); the resulting power indices are higher, but they do not reach the correct order of magnitude. Finally, we are completely unable to reproduce the soft ENA behaviour in the heliotail, which is a region we obviously have not modelled with sufficient accurary yet.

Future studies will aim to improve these results further, by refining our heliospheric boundary model to allow more realistic magnetic fields, a more realistic shape of the boundary layer (which right now is approximated by a triaxial ellipsoid), and also transport effects. The last point might be able to explain the steeper drop of ENA spectra in the heliotail. At TS distances of 100 AU, is has been estimated that about 90% of ENAs generated may reach an Earthbound observer (Lee et al. 2009); however, more recent studies prove that this effect is not negligible over two decades of energy (see e.g. Bzowski 2008; Schwadron et al. 2011).

Since this effect is highly non-trivial, we are currently preparing a separate study that explicitly addresses ENA propagation (i.e. basic extinction, solar cycle effects, local effects such as CMEs and other density fluctuations). Since this is beyond the scope of this study, we conclude by noting that several aspects discussed here will also influence the anounced work. For example, ENA reionisation can be parameterised as (56)which however is an incomplete version of the radiation transport equation (the full equation also includes a source term). This equation strongly depends on the charge exchange cross section σ, which must then be used in a more specific parameterisation, and also on solar wind density fluctuations, which might result from the solar cycle and also local effects such as ENAs crossing CMEs propagating through the inner heliosphere.

Acknowledgments

M. Siewert is grateful to the Deutsche Forschungsgemeinschaft for financial support granted to him in the frame of the project Si-1550/2-2. D. J. McComas and N. A. Schwadron were supported by NASAs Explorers Program as a part of the IBEX mission.

References

All Figures

	Fig. 1 Representative plot of the source PUI function (at 70 AU) used in this study, using U = 1 as the initial PUI velocity. Depending on the normalisations of the primary and secondary contributions, a different behaviour is possible. From Fahr et al. (2009).
In the text

	Fig. 2 Schematic representation: vectors used in the velocity transformation, not to scale.
In the text

	Fig. 3 Correction factor d(r) (see Eq. (15)) to a pure power law due to the frame-of-reference effect, assuming a radial solar wind speed of Ur = 100 km s-1.
In the text

	Fig. 4 Schematic representation: the coordinate system used in the potential flow model for the plasma stream lines, not to scale.
In the text

	Fig. 5 Schematic representation: The behaviour of streamlines characterised by the solar wind speed for solar minimum conditions, not to scale. Obviously, a modification can be expected in the heliotail region, where lines-of-sight cross deflected streamlines.
In the text

	Fig. 6 ENA flow skymaps at Eobs = 1 keV for both source components. Highlighted points are used as reference positions for the spectral study in Sect. 4.2.
In the text

	Fig. 7 ENA spectral flows as a function of energy at selected positions in the sky, for PUI-induced ENAs (blue line), ACR-ionduced ENAs (green line), and the sum of both (red line).
In the text

	Fig. 8 Modified ENA flux spextra (j(E)·E1.5) at the reference points as defined in Fig. 7, and represented graphically in Fig. 6.
In the text

	Fig. 9 Modified ENA flux spextra (j(E)·E1.5) in the ring region (α,δ) = (−90°,0°) for the fast and slow solar wind regions.
In the text

	Fig. 10 Modified ENA flux spectra for a heliotail HP cutoff distance of 1200 AU as a first approximation to the model developed by Jaeger & Fahr (1998).
In the text