NEUTRAL INTERSTELLAR HELIUM PARAMETERS BASED ON IBEX-Lo OBSERVATIONS AND TEST PARTICLE CALCULATIONS

To achieve the highest possible realism and fidelity of the simulations, the simulation pipeline accurately addressed all relevant geometry, timing, and instrumental aspects including the following:

• 1.  
  The Visible sky Strip. IBEX-Lo observes a strip on the sky almost exactly perpendicular to the IBEX rotation axis. The field of view (FOV) of IBEX-Lo was defined in the simulation program according to the FOV specification for IBEX-Lo (Fuselier et al. 2009), and in a separate study (Hłond et al. 2012), it was verified that the pointing of IBEX-Lo indeed agrees with its specified pointing in the spacecraft system to better than 015. The spin axis of IBEX is close to the ecliptic during science operations and pointed <1° above (Figure 1). For each orbit the Visible Strip of sky viewed by the sensor was calculated based on the exact pointing of the IBEX spin axis determined by the IBEX Science Operations Center (ISOC) (Schwadron et al. 2009) and illustrated in Figure 1.
• 2.  
  Collimator shape and transmission function. Transmission function T(ρ, θ) of the collimator was adopted from pre-launch calibration (Fuselier et al. 2009); its shape is shown in Figure 2. The value of the transmission function at a given location within the FOV of the sensor, described by the polar coordinates (ρ, θ) relative to the boresight axis of the collimator, corresponds to the percentage of the flux that is able to enter the sensor at a given area element sin ρ dρ dθ. The FOV of the low angular resolution portion of the collimator of IBEX-Lo is hexagonal in shape (Fuselier et al. 2009), and its arrangement relative to the sky strip scanned during spacecraft rotation is shown in Figure 3. The profiles of the transmission function from the boresight to the corner and to the side of the hexagonal base were fitted with third- and second-order polynomials, respectively. Both the shape of the collimator transmission function and its arrangement in the spacecraft reference system were exactly simulated in the program.
• 3.  
  Positions and velocity of Earth relative to the Sun. We used the ephemeris obtained from the SPICE-based program developed and operated by the ISOC (Acton 1996; Schwadron et al. 2009). These include actual solar distances and ecliptic longitudes of Earth as well as Earth velocity vectors relative to the Sun for all dates for which simulations were performed.
• 4.  
  Velocity vectors of the IBEX satellite relative to Earth. They were taken from the same SPICE-based program source as for the Earth orbit; together with the Earth information these spacecraft velocity vectors were used to calculate the state vectors of IBEX relative to the Sun for the simulations.
• 5.  
  Selection of observations. We use the observations selected from the IBEX-Lo data set with data dropouts, spacecraft pointing knowledge problems, and other spacecraft and sensor conditions that affect overall flux and direction removed (see Möbius et al. 2012 for a detailed description of the select ISM flow observation times).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulation pipeline accepts the following as input: parameters of the NISHe gas in the LIC, energy limits of the incoming atoms in the spacecraft inertial frame to be adopted as flux integration boundaries, parameters that describe heliospheric conditions (time series of the photoionization rate and solar wind density and velocity averaged over Carrington rotations), the number of the orbit for which the simulation is to be performed (i.e., dates and times of the simulation), the spin-axis pointing for the orbit, the list of select ISM flow observation times for the orbit from Möbius et al. (2012), and the state vectors of the IBEX spacecraft relative to the Sun for the observations. It returns collimator-averaged fluxes of the NISHe gas as a function of IBEX spin angle averaged by the selected times and the collimator transmission function. The simulation pipeline products can be directly compared with the observed count rates for the given orbit after linear scaling. In the simulations carried out for this study, the integration boundaries were adopted from zero to infinity, so effectively the integration was over the full energy range of the incoming NISHe atoms. As discussed by Möbius et al. (2012), such an approach is valid because IBEX-Lo actually does not measure incoming He atoms directly; rather, it detects H, O, and C atoms sputtered off the conversion surface, so He atoms of all relevant energies contribute significantly to the sputtered H signal collected by the energy steps 1, 2, and 3 (energy passbands between 0.01 and ∼0.075 keV) of IBEX-Lo. Details of calibration of the IBEX-Lo instrument for detection of H, He, and O atoms are provided by Saul et al. (2012) and Bochsler et al. (2012).

The core of the simulation program calculates the flux of NISHe relative to the IBEX spacecraft located at a point $\bm r$ relative to the Sun, traveling at a velocity $\bm v$ at a time t for a line of sight determined by ecliptic coordinates (λ_LOS, β_LOS). This part of the simulation set is described in the following section. Here we discuss simulations of the NISHe flux averaged over the IBEX-Lo collimator FOV and select ISM flow observation times in a given orbit.

The simulation pipeline is organized as follows. With the select ISM flow observation times transformed into Julian days, a series of dates at halves of full Julian days that straddle and fill in the selected intervals is determined. Subsequently, the Visible Strip is determined based on the spin-axis pointing for the given orbit, and coordinates of its boundaries in the ecliptic reference system are calculated. The Visible Strip is then transformed into a heliographic reference system (HGI; Fränz & Harper 2002) and mapped on a grid of equal-area, equidistant pixels whose boundaries and centers in the heliographic coordinates are adopted following the HealPix scheme (Górski et al. 2005) with the resolution parameter N = 64, which corresponds to 49,152 pixels for the whole sky. Thus, the angular resolution of the Visible Strip coverage is better than 1 deg². Centers of these pixels make the simulations mesh (λ_LOS, β_LOS).

The Visible Strip does not change during one orbit, so during all select ISM flow observation times in a given orbit the instrument is looking at the same portion of the sky. The simulations are carried out for all pixel centers within the Visible Strip for all select ISM flow observation times in a given orbit in the inertial frame of the IBEX spacecraft. The inertial frame is determined by the IBEX velocity vector $\bm v$ relative to the Sun, which is obtained from the ISOC.

With the detailed map of the NISHe flux within the Visible Strip for a given day, we calculate the flux transmitted through the collimator. We do so by sliding the collimator boresight along the spin angle ψ in 1° steps (see Figure 3), integrating the flux as convolution of the transmission function T(ρ, θ) with the flux F_He(ρ(ψ), θ(ψ)):

$$\begin{eqnarray} &&F_{\mathrm{He,\,coll}}(\psi)=\int \nolimits _{\theta =0}^{2 \pi }\int \nolimits _{\rho =0}^{\rho _{1}(\theta)}\! T(\rho, \theta)F_{\mathrm{He}}(\rho (\psi), \theta (\psi)) \sin \rho \,{d} \rho \,{d} \theta,\nonumber\\ \end{eqnarray} \tag{ 1 }$$

where F_He(ρ(ψ), θ(ψ)) is the flux calculated at (ρ, θ) for a given spin phase angle ψ of the collimator boresight and ρ₁(θ) describes the hexagonal boundary of the FOV. For a different boresight ψ the same flux element will be located differently relative to the boresight direction and consequently will contribute to the collimator-averaged flux with a different weight. In practice, the collimator FOV was divided into regions of approximately equal areas distributed symmetrically around the boresight at a series of (ρ_i, θ_j). The integration over the FOV was in fact a summation of the flux with appropriate weights:

$$\begin{eqnarray} &&F_{\mathrm{He,\,coll}}(\psi)={\sum _{i=1}^{N_i}}{\sum _{j=1}^{N_j}}{\sum _{k=1}^{N_{ij}}}T(\rho _i, \theta _j) F_{\mathrm{He}}(\rho _{ik}, \theta _{jk})s_{ij}/N_{ij},\quad \end{eqnarray} \tag{ 2 }$$

where i marks the radial and j the azimuthal index of the mesh, s_ij is the unity-normalized area of the i, j field, and k counts from 1 to N_ij the pixels at (ρ_ik, θ_jk) in the (i, j)th field, in which the field-averaged F_He flux is calculated. Since the s_ij fields are equal area, the number of sky pixels per integration field is approximately constant, which adds to the numerical stability of the calculation scheme.

Following the procedure described in the preceding paragraphs, we obtain a series of collimator-integrated fluxes for given days, which subsequently are time averaged over the select ISM flow observation times. The result of this averaging is taken as the simulation result for a given set of parameters of NISHe gas for a given orbit. The procedure of calculating the collimator- and orbit-averaged flux was repeated for all orbits within the 2009 and 2010 observation seasons.

In the inertial frame of IBEX, the flux of NISHe gas $F_{\mathrm{He}}\left(\lambda, \beta, \bm r,t\right)$ that goes into the ecliptic-coordinates direction (λ, β) at the location described by the heliocentric vector $\bm r$ at a time t is calculated by

$$\begin{eqnarray} F_{\mathrm{He}}(\lambda, \beta, \bm r, t)&=&\int \nolimits _{0}^{\infty }\! v_{\mathrm{He,\,sc}}f_{\mathrm{He}} (\bm v_\mathrm{{He,\,ecl}}, \bm r, t)\hat{\bm e}(\lambda, \beta )\nonumber\\ &&\times v_{\mathrm{He,\,sc}}^2 {d}v_{\mathrm{He,\,sc}}, \end{eqnarray} \tag{ 3 }$$

where v_{He, sc} is the magnitude of the He atom velocity vector $\bm v_{\mathrm{He,sc}}$ in the inertial frame of the spacecraft, $\hat{\bm e}\left(\lambda, \beta \right)$ is the unity vector pointing toward (λ, β), and $f_{\mathrm{He}}\left(\bm v_{\mathrm{He,ecl}},\bm r,t\right)$ is the distribution function of the NISHe gas for time t and solar-frame velocity $\bm v_{\mathrm{He,ecl}}$ at the location specified by the solar-frame radius vector $\bm r$. Assuming that the flow of the NISHe gas in the LIC is constant, the distribution function f_He at $\bm r$ is time dependent only because of variations in the helium ionization rate.

The transition from the spacecraft inertial frame to the solar inertial frame is done by a simple vector subtraction: with the IBEX solar-inertial velocity $\bm v_{\it{IBEX}}\left(t\right)$ the relation between the IBEX-inertial $\bm v_{\mathrm{He,sc}}$ and solar-inertial $\bm v_\mathrm{He,ecl}$ velocities is

$$\begin{equation} \bm v_{\mathrm{He,ecl}}\left(t\right)=\bm v_{\mathrm{He,sc}}-\bm v_{\mathrm{}{\it IBEX}}\left(t\right). \end{equation} \tag{ 4 }$$

The conversion to the solar-inertial frame during the integration specified in Equation (3) is done separately for each value of v_{He, sc}, and the calculation of the local distribution function f_He is performed in the solar-inertial frame.

The model of neutral interstellar gas flow in the inner heliosphere, used to calculate the maps of NISHe flux at Earth orbit, is a derivative of the Warsaw Test Particle Model developed since the mid-1990s (Ruciński & Bzowski 1995). Previous versions, as well as its development history, are found in Tarnopolski & Bzowski (2008b). Recent applications of this code in interpreting measurements of neutral interstellar hydrogen (NISH) in the inner heliosphere are discussed by Bzowski et al. (2008, 2009), and its use in interpreting interstellar helium measurements by Gloeckler et al. (2004). Predictions of neutral interstellar deuterium flux at IBEX, obtained using the model, can be found in Tarnopolski & Bzowski (2008a). Details of test particle calculations of NISHe in the inner heliosphere are in Ruciński et al. (2003). The model was used by Möbius et al. (2009a) to verify the detection by IBEX of the NISHe atoms.

In order to be used in the determination of the flow parameters of the NISHe gas from IBEX-Lo observations, the model had to be modified. Modification to the model was done in three main areas: (1) atom dynamics, (2) inertial frame, and (3) ionization rate as a function of time and location in the heliosphere.

The first modification was the most straightforward: since the resonance radiation pressure force acting on the neutral He atoms in the heliosphere is practically negligible, the radiation pressure module in the code could be switched off. The atoms now move solely under the 1/r² solar gravity force. This change greatly simplified the requirements for the atom tracking module. Nevertheless, this module still had to maintain its ability to accurately link the time on orbit with the locus on orbit, and the current sophisticated Runge-Kutta tracking scheme was not replaced to save on the development time and maintain sufficient homogeneity of the code in view of planned future applications of the model to interstellar hydrogen analysis. A version of the code with the full radiation pressure module installed was used to calculate the predictions of the neutral interstellar H signal discussed later on in the paper.

Since the calculation of the NISHe flux needs to be done in the IBEX spacecraft inertial frame, the input direction in space and speed of the atom are formulated in the moving frame and transformed to the solar frame. Therefore, initial values of the atom velocity are taken relative to the Sun, not to the spacecraft. The integration over speed, which returns the flux relative to the spacecraft from a given direction in space, is performed in the spacecraft reference frame, but parameters of the integrand function are converted to the solar inertial frame in the heliographic reference system. This change to the solar HGI frame is needed because the ionization rate model, which is used to calculate the survival probability of the atom, uses the solar equator as the natural reference plane.

The transformation from the spacecraft-inertial frame to the solar-inertial frame requires only specifying the velocity vector of the spacecraft at the desired moment of time. No further assumptions need to be made, which facilitates adoption of various spacecraft velocity vectors in the calculation scheme.

The ionization rate, which is discussed in greater detail in the following section, is time dependent. We determine all the quantities relevant for the calculation of the net ionization rate as a function of time by interpolating between Carrington-period averaged quantities. Thus, the model is fully time dependent and uses current best parameters obtained directly from measurements, which adds to the accuracy of the results.

Helium has the highest first ionization potential of all elements (27.587 eV), and hence the ionization losses of the NISHe gas in the heliosphere are relatively low. Where IBEX makes its measurements (at 1 AU), as much as 70% of the atoms from the original population are able to survive (Ruciński et al. 2003). Nevertheless, ionization has to be taken into account in the analysis because it modifies the shape of the observed helium beam. Ionization changes the apparent velocity distribution of the NISHe beam because it more readily removes slower atoms from the ensemble than faster ones, and thus the mean velocity vector of the remaining distribution differs from the conditions when no ionization is operating (this effect is much more pronounced for hydrogen and was discussed in this context by Lallement et al. 1985 and Bzowski et al. 1997). The selective ionization results in a change in the ecliptic longitude at which the maximum of the NISHe beam is observed by a few tenths of a degree and, if unaccounted for, biases the derived speed and longitude of the flow direction. Similarly, this effect reduces the width of the beam somewhat, which, if neglected, leads to an underestimation of the temperature.

Heliospheric conditions that affect the flow of the NISHe gas in the inner heliosphere were extensively discussed by McMullin et al. (2004). The dominant ionization process is solar photoionization, which varies throughout the solar cycle from about 5.5 × 10⁻⁸ s⁻¹ at minimum to ∼1.5 × 10⁻⁷ s⁻¹ at maximum. In the present study, following P. Bochsler et al. (2012, in preparation), we adopted the cross section after Samson et al. (1994) and Verner et al. (1996) and directly integrated the spectra obtained from TIMED/SEE (Woods et al. 2005). We verified the agreement of the results with the measurements from CELIAS/SEM (Judge et al. 1998). As seen in Figure 4, measurements of the NISHe flow in the 2009 season followed a prolonged period of very low solar activity and very stable photoionization rate. In contrast, measurements in the 2010 season occurred during a period of increasing activity, with the photoionization rate higher by 15% than during the preceding measurement season.

Zoom In Zoom Out Reset image size

As pointed out by Auchère et al. (2005a, 2005b), the photoionization rate of helium appears to vary weakly with heliolatitude, with the polar rate probably being about 80%–85% of the equatorial rate. This latitudinal variation was accounted for by implementing the following relation:

$$\begin{equation} \beta _{\mathrm{ph}}\left(\phi \right)=\beta _{\mathrm{ph}}\left(0\right)\sqrt{a_{\beta _{\mathrm{ph}}}\sin ^{2}\phi +\cos ^{2}\phi }, \end{equation} \tag{ 5 }$$

where $a_{\beta _{\mathrm{ph}}}$ is the latitudinal "flattening" parameter, adopted to be 0.8. In test simulations we verified, however, that this flattening has a small effect on the expected NISHe flux in the heliospheric tail region and practically no effect at the interval of ecliptic longitudes where IBEX measurements were taken. The weakness of this effect can be easily explained by the fact that the trajectories of NISHe atoms detected by IBEX-Lo remain close to the ecliptic throughout their travel from the LIC to Earth orbit and therefore never experience the ionization rates relevant for higher latitudes.

Another ionization process of neutral helium is ionization by impact of solar wind electrons. The importance of this ionization process for NISHe in the heliosphere was first pointed out by Ruciński & Fahr (1989, 1991). As discussed by McMullin et al. (2004), who used more recent measurements of solar wind electrons, this rate close to the ecliptic plane at 1 AU from the Sun is equal to about 2 × 10⁻⁸ s⁻¹, i.e., it is at an appreciable level of ∼30% of the photoionization rate, but owing to the rapid cooling of the solar wind electrons, it falls off with solar distance much faster than 1/r², i.e., faster than the drop-off of the photoionization rate.

We expanded the electron-ionization model used by McMullin et al. (2004) assuming the thermal behavior of solar wind electrons as conforming to the core+halo model (Pilipp et al. 1987). Following the approach adopted by Bzowski (2008) to develop an electron ionization rate model for hydrogen, we used the solar wind electron temperature and density measurements by Scime et al. (1994), Issautier et al. (1998), and Maksimovic et al. (2000) and implemented the model by Ruciński & Fahr (1989, 1991), where the cross section for electron-impact ionization by Lotz (1967) is convolved with the Maxwellian distribution function separately for the core and halo temperatures, assuming the radial dependence of the temperatures and the proportions between the core and halo population densities as compiled by Bzowski (2008). The total electron density was implemented as tied to the density of solar wind protons (enhanced by the doubled average alpha particle abundance). The radial behavior of thereby obtained ionization rates for the 2009 and 2010 seasons is presented in Figure 5. The electron-impact ionization is important only in the final phase of an NISHe atom flight before detection by IBEX, when its distance from the Sun is close to 1 AU. One has to note, however, that because IBEX measures only atoms near their perihelia, i.e., those which travel nearly tangentially to the 1 AU circle around the Sun, the influence of electron-impact ionization is stronger than when they are observed at ecliptic longitudes in the upwind hemisphere.

Zoom In Zoom Out Reset image size

The least significant ionization process for neutral helium in the inner heliosphere is charge exchange with solar wind particles: protons and alphas (Ruciński et al. 1996, 1998; McMullin et al. 2004). While significantly less intense, we include this process for completeness. The instantaneous charge exchange rate is defined by McMullin et al. (2004) in their Equations (2)–(4), from the formula

$$\begin{eqnarray} \beta _{\mathrm{He,\,cx}}(t)&=&n_{\mathrm{p}}(t)|\bm v_{\mathrm{HeENA}}-\bm v_{\mathrm{SW}}| [2\alpha _{\mathrm{\alpha }}\sigma _{\mathrm{He,\alpha }}(|\bm v_{\mathrm{HeENA}}-\bm v_{\rm SW}|) \nonumber \\ && +\,\sigma _{\mathrm{He,p}}(|\bm v_{\mathrm{HeENA}}-\bm v_{\mathrm{SW}}|)], \end{eqnarray} \tag{ 6 }$$

where $\left|\bm v_{\mathrm{HeENA}}-\bm v_{\mathrm{SW}}\right|$ is the relative speed between a He atom at $\bm v_{\mathrm{HeENA}}$ and the radially expanding solar wind at $\bm v_{\rm SW}(t)$, α_alpha ≈ 0.04 is a typical abundance of solar wind alphas relative to protons, n_p(t)) is the local proton density taken from the OMNI-2 compilation of solar wind observations (King & Papitashvili 2005), and σ_{He, p}, σ_{He, α} are the charge exchange cross sections for the reaction given by Equation (2) and a sum of reactions given by Equations (3) and (4) by McMullin et al. (2004). The net rate from the three charge exchange processes that were taken into account is ∼2.6 × 10⁻⁹ s⁻¹ regardless of the activity phase, which is of the order of 4% of the typical photoionization rate. Thus, typically the rate of charge exchange losses is less than the uncertainty in the photoionization rate. We implemented it only to make sure that we do not miss a sudden increase in total ionization rate due to possible high-flux events in the solar wind (like coronal mass ejections), when the solar wind density may increase by an order of magnitude.

Since IBEX-Lo is able to observe helium only indirectly via sputtering products from the conversion surface, which include H atoms (Möbius et al. 2009a, 2009b, 2012; Saul et al. 2012), we determined from the simulation in which orbits the flux expected from the NISHe flow should exceed the flux expected from NISH. We compared collimator-averaged total NISHe fluxes expected assuming the prior consensus NISHe flow parameters (Möbius et al. 2004), which are very close to the parameters obtained by Witte (2004) from Ulysses, with the NISH flux in energy step 2 (center energy 27 eV) of the IBEX-Lo detector (Fuselier et al. 2009). For this comparison, we assumed that the population of interstellar hydrogen at IBEX is a mixture of the primary population of interstellar hydrogen and a secondary component due to charge exchange with the heated and compressed plasma in front of the heliopause (Malama et al. 2006). We used the parameters of the two populations as determined by Bzowski et al. (2008) based on pickup ion measurements on Ulysses (Gloeckler et al. 2008).

As shown in Figure 6, in orbit 11 the helium signal dominates and the only appreciable NISH flux (two orders of magnitude lower than the He flux) is from the secondary population. The dominance of He over H increases from orbit 11 to 17, but the intensity of the primary H population gradually increases and in orbit 17 it exceeds the peak intensity of the secondary hydrogen. However, He fluxes are still significantly higher than the combined hydrogen fluxes. Starting in orbit 14, the wings of the H signal become wider than the wings from He, but these wings are more than three orders of magnitude lower than the peak of the He flux. The situation changes in orbit 20, when the hydrogen primary population is only a few times weaker than He and thus might appear as an extra component in the total signal. The secondary hydrogen wings should be at a level of ∼1% of the He peak. In orbit 21 (when IBEX is viewing the nose of the heliosphere) H exceeds He, and in orbit 22 H becomes dominant. The observation of interstellar H is discussed by Saul et al. (2012).

Zoom In Zoom Out Reset image size

Even though a change in the solar wind or interstellar parameters may change details, the basic conclusion is that the best orbits to study the NISHe flow are orbits 13 through 19–20 and their equivalent during the second ISN season for IBEX (see Figure 1; further justification is provided in the data selection section). Since the NISHe population is highly peaked and at the peaks it exceeds the H populations by more than three orders of magnitude, it is appropriate to analyze the Gaussian cores of the signal as due solely to the NISHe flow. Since the NISHe flow should be mostly visible at the wings and since it is expected to consist of at least two populations, making the signal fairly complex, we decided to remove these non-Gaussian wings from the NISHe analysis.

To differentiate the signal from background, secondary populations, and other potential biasing, we investigated in greater detail how the collimator-averaged signal would appear if we assume no background or secondaries and further assume that the NISHe gas distribution function in the LIC is the purely Maxwellian function:

$$\begin{equation} f_{\mathrm{He,\,Maxw}}\left(\bm v\right)=n_{0}\left(\frac{m_{\mathrm{He}}}{2\pi kT}\right)^{3/2}\mathrm{exp}\left[-\frac{m_{\mathrm{He}}}{2kT}\left(\bm v-\bm v_\mathrm{B}\right)^2\right], \end{equation} \tag{ 7 }$$

with the density n₀, temperature T, and a shift in phase space by the bulk velocity $\bm v_\mathrm{B}$. We further assume that instantaneous observations with high spin-phase resolution are performed during various orbits in one observation season at the moments when the ecliptic longitude of the spin axis of the IBEX satellite is precisely equal to the ecliptic longitude of the Sun. We will refer to such conditions as the Exact Sun-Pointing (ES) conditions.

Simulations performed for a number of parameter sets that covered the expected range of the parameters of the NISHe gas in the LIC suggest that at the orbits where the helium signal is expected to be the strongest (i.e., from orbits 13–20 and the equivalent ones during the later seasons) the observed count rate as a function of spin angle ψ can be approximated by a Gaussian core:

$$\begin{equation} F_{\mathrm{obs}}\left(\psi \right)=f_{0}\,\exp \left[-\left(\frac{\psi -\psi _{0}}{\sigma }\right)^2\right], \end{equation} \tag{ 8 }$$

with elevated wings. This is illustrated in Figure 7 for three selected orbits and three different parameter sets. The parameters of the Gaussians (peak height f₀, peak width σ, and spin angle of the peak ψ₀) depend on the choice of parameters of the NISHe gas in the LIC, but the feature of a Gaussian core and elevated non-Gaussian wings is always present. The Gaussian core is a result of convolution of the true Gaussian signal with the near-Gaussian transmission function of the collimator. Fits of the Gaussian function to the simulation results showed that residuals of the fits within the Gaussian core region were below 1%. Outside the Gaussian core region, whose span in the spin angle varied with assumed bulk velocity and temperature, the elevated non-Gaussian wings were visible in the residuals as power-law increase in the residual magnitudes. They were present for the collimator-integrated flux values F_He(ψ) ≲ 0.01F_He(ψ_max), where ψ_max is the spin phase angle of the peak flux, as illustrated in Figure 7.

Zoom In Zoom Out Reset image size

We also found that for orbits earlier than 12 (and equivalent in 2010) the simulated collimator-averaged signal increasingly deviates from the Gaussian shape with the decrease of Earth's ecliptic longitude. The flux profiles as a function of spin angle become increasingly asymmetric relative to the peak even though the Maxwellian distribution function in the LIC is assumed, as shown by the blue line in the upper left panel of Figure 6 (orbit 11). Despite the non-Gaussianity of the profiles, their peaks are well defined and can be easily compared with observations. Such comparisons were in fact done and used as a basis to formulate the hypothesis that the excess signal observed by IBEX at this Earth longitude interval is due to an additional population of neutral He in or near the heliosphere.

We further verified that binning data into 6° bins does not remove the Gaussian character of the signal, as shown in Figure 8. Similarly, averaging of the signal over the entire duration of the select ISM flow observation times maintains the Gaussian shape, but the parameters of the Gaussians (peak height, peak width, and peak location) are changed, as illustrated in Figure 9, where simulations performed for the ES conditions are compared with simulations performed for the actual select ISM flow observation intervals.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The reason for the differences between the select ISM flow observation times and ES beam parameters is that because the spin axis of the spacecraft, which is never aligned with the Sun, does not change during an orbit (Scherrer et al. 2009; Hłond et al. 2012), the beam of the NISHe gas, which in the solar inertial frame is invariant relative to the distant stars, wanders through the FOV of the sensor, changing gradually its angular size, peak location, and height. This effect is especially visible in the orbits before or after orbits 16 and 64 and is illustrated in Figures 10–12. These figures demonstrate the importance of an exact determination of select ISM flow observation times in order to have a faithful representation of the data in the simulations. These figures also demonstrate why Möbius et al. (2012) had to extrapolate their observations to the ES conditions for the comparison with their analytic model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the select ISM flow observation times extended over the entire orbit, then IBEX-Lo would have observed daily fluxes marked by the thin lines, which, when averaged, would equal the thick blue lines. However, these times do not extend over the entire orbit. In Figure 10, the flux of the incoming interstellar He atoms is most intense during the first days of the orbit and with time the beam moves away from the FOV of the collimator. Since the select ISM flow observation times cover the last portion of the orbit, a lower average flux is observed, as illustrated with the thick red line. However, the spin axis pointed toward the Sun at the beginning of the orbit, so the flux relevant for the ES conditions, marked with the thick purple line, is much higher than the flux actually measured.

During orbit 16 (Figure 11) IBEX observed the peak NISHe flux. The select ISM flow observation times occur during the ∼3 days at the beginning of the orbit. However, since this is the peak flux and the NISHe beam is directed into the sensor, the flux varies little with time and the observed mean flux is very similar to the flux averaged over the entire Science Operations for this orbit. Thus, in Figure 11 the ES flux (purple), the observed average flux (red), and the orbit average flux (blue) are very similar.

In orbit 18 (Figure 12) IBEX is beyond the peak NISHe flux and viewing the beam edge. For this orbit, select ISM flow observation times occurred during the first ∼5 days of Science Operations. The beam moves into the FOV near the middle of the orbit, but IBEX views the beam when it is off the peak and the average flux is lower for the selected times than for the full orbit. The spin axis pointed exactly to the Sun at the beginning of the orbit, so the flux at the ES time is lower than the flux averaged over the Select Times and lower than the flux averaged over the entire orbit.

From this analysis we conclude that the portions of the observed count rates that are Gaussian in shape correspond to the NISHe population from the LIC and the portions that cannot be fitted by a Gaussian must correspond to something different, probably another source of neutral helium in or near the heliosphere. In either case, these non-Gaussian components are eliminated from the analysis of the NISHe population for now. We also conclude that care must be taken to accurately reproduce the flux observed during select ISM flow observation times, especially for orbits that are not near the peak flux.

Finally, before starting the parameter fitting procedure, we discuss miscellaneous effects that should be included in the simulations. These effects are listed at the beginning of Section 2.

The ellipticity of Earth's orbit results in a small deflection of the direction of the Earth velocity vector from the right angle to the Earth radius vector, which slightly modifies the aberration of the NISHe beam. Further, an additional change in the aberration and relative velocity of the beam and the detector is caused by the small radial component of Earth's velocity (on the order of 1 km s⁻¹). Also of the order of a few km s⁻¹ is the proper motion of IBEX relative to the Earth. In the simulations we used actual Earth ephemeris, which accounts for the Earth location. The total velocity vector of the spacecraft plus the Earth is accounted for by using the proper motion of the satellite along its orbit and the total velocity vector of the Earth.

As shown in Figure 13, the IBEX motion relative to the Earth has its strongest effect on the peak height of the observed NISHe beam. Only the peak height effect exceeds the measurement uncertainty. The effect on peak width is, understandably, negligible, and the effect on peak spin angle is comparable to the measurement uncertainty. Since the effect on the magnitude of the flux cannot be neglected, the satellite proper motion was included in the simulations. The velocity vector of the spacecraft in the inertial frame of the Sun was taken as a vector sum of the Earth velocity relative to the Sun and IBEX's velocity about the Earth and was calculated using the software developed by the ISOC (Schwadron et al. 2009) based on the SPICE toolkit (Acton 1996). It should be noted that the aberration effect is stronger during the ES time for each orbit, because that occurs during the ascent of IBEX to apogee when the spacecraft speed is still substantial, and therefore the effect is also to be taken into account in the analysis by Möbius et al. (2012).

Zoom In Zoom Out Reset image size

The simulations shown in Figure 13 were done for the NISHe flow parameters established in this paper based on fitting of the model with all the effects included. It is not surprising then that the simulations without the IBEX orbital velocity fit the data less well. Since we know that IBEX is moving in its orbit and we know from Figure 13 that the influence of this effect on the observed fluxes is small, but not negligible, we include these effects in the simulations.

The small tilt of the spin axis out of the ecliptic also affects the observed flux because it excludes a small portion of the beam while accepting a different part compared to the situation when the spin axis is exactly in the ecliptic plane. This effect is especially pronounced in early orbits before the crossing of the ISM flow peak (Figure 14). Both the width and spin phase of the peak maximum are affected with offsets clearly larger than the error bars. In contrast, the peak height is only weakly affected.

Zoom In Zoom Out Reset image size

The exact magnitudes of the effects discussed in this section depend on the details and are challenging to plan in advance, i.e., on the actual select ISM flow observation times, which are determined by a combination of operational aspects, stochastic backgrounds, particle events, actual realizations of spin-axis repointing maneuvers, etc. We decided that instead of attempting to correct the observations for all of these issues, it was better to simply include them in the simulation. We emphasize that the magnitude of various effects can vary depending on the adopted parameter set and also from orbit to orbit. This supports the decision to complicate the simulation pipeline for the sake of fidelity of the model rather than try to correct the observations.

The goal of our analysis is to determine the flow direction, velocity, and temperature of the NISHe gas in the LIC ahead of the heliosphere. We accomplished this by fitting simulations of the NISHe flux to the data, with the ecliptic longitude and latitude of inflow direction, inflow speed, and gas temperature in the LIC as free parameters. Optimizing a multi-parameter model fit to data usually involves selecting a merit function whose free parameters are the fitted model parameters and finding its minimum in the multi-dimensional parameter space. Colloquially speaking, the merit function describes the "distance" of the model predictions from the data in the observation N-space, and searching for the best parameters requires finding the parameter set for which this distance is minimum.

A well-tested and widely used method of fitting parameters of a model to a data set is the maximum likelihood method. In this method, the merit function is the likelihood function. To use it, one needs to know probability distributions $f_{p,i}(x_{p,i}, d_i, \bm p)$ of all the data points d_i, parameterized by the model parameters $\bm p$. In principle, probability distributions for different data points can be described by different probability distribution functions, but in our case we assume that for all data points they are identical, i.e., for all i, $f_{p,i}(x_{p,i}, d_i, \bm p) = f_{p}(x_{p,i}, d_i, \bm p)$.

With these definitions we calculate the conditional probability P_i that if the model with a given parameter set $\bm p$ is correct, then our experiment in case i provides measurement d_i, given by the formula

$$\begin{equation} P_{i}(x_{p,i})=f_{p}(x_{p,i},d_i,\bm p). \end{equation} \tag{ 10 }$$

The series x_{p, i} is the series of model predictions of the measurements for parameters $\bm p$. The probability P that our series of N measurements returns a series of results d_i, i = 1, ..., N is, of course, a product of all N probabilities P_i:

$$\begin{equation} P(x_{p,1},\ldots,x_{p,N},d_1,\ldots,d_N)=\prod _{i=1}^{N}{P_i}=\prod _{i=1}^{N}{f_{p}(x_{p,i},d_i,\bm p)}. \end{equation} \tag{ 11 }$$

Fitting the parameters $\bm p$ is equivalent to finding the parameters $\bm p_{\mathrm{best}}$ for which absolute maximum of P is achieved. Finding this absolute maximum is the basis of the maximum likelihood method. Remaining details determine how to best accomplish the goal, and the mathematical methods to apply depend on the nature of the problem on hand.

The IBEX-Lo detector actually counts incoming NISHe atoms in 6° spin angle bins, so the number of atoms in each bin is subject to Poisson statistics. Hence, we immediately have estimates of the measurement errors according to

$$\begin{equation} \sigma _{i}=\sqrt{d_i}. \end{equation} \tag{ 12 }$$

But the counts are relatively high and in this case the Poisson statistics asymptotically transforms into the Gaussian. Thus, the likelihood function in Equation (11) becomes

$$\begin{eqnarray} &&P\left(x_{p,1},\ldots,x_{p,N},d_1,\ldots,d_N\right)\nonumber \\ &&\quad =\prod _{i=1}^{N}{\left(\frac{1}{\sqrt{\pi }\sigma _i}\exp \left[-{\left(\frac{d_i-x_{p,i}}{\sigma {i}}\right)}^2\right]\right)} \end{eqnarray} \tag{ 13 }$$

which we must maximize. Since all the probabilities are positive numbers, we can take the natural logarithm of both sides of this equation and obtain

$$\begin{eqnarray} &&\ln [P(x_{p,1},\ldots,x_{p,N},d_1,\ldots,d_N)]\nonumber \\ &&\quad = \sum _{i=1}^{N}\left[\ln \left(\frac{1}{\sqrt{\pi }\sigma _i}\right)-{\left(\frac{d_i-x_{p,i}}{\sigma {i}}^2 \right)}\right]. \end{eqnarray} \tag{ 14 }$$

Since for a given measurement series the first term under the logarithm in the sum in Equation (14) is a constant, we can remove it because our goal is to find the parameter set for which the likelihood function will be maximum, and not the maximum value itself. Thus, we define the following merit function $-L\left(\bm p\right)$:

$$\begin{equation} -L\left(\bm p\right)=-\sum _{i-1}^{N}{{\left(\frac{d_i-x_{p,i}}{\sigma _i}\right)}^2} \end{equation} \tag{ 15 }$$

which takes negative values. We can omit the minus signs, and then instead of maximizing the term at the right-hand side, we have to minimize it. The parameter set $\bm p$ for which function $L(\bm p)$ is minimal will not change when we divide it by the number of degrees of freedom in the problem equal to N − n_p, where n_p is the number of parameters in the parameter set $\bm p$. In our case n_p = 4. Dividing by the number of degrees of freedom converts this function into the chi-squared function and enables direct comparison of the quality of approximation between data series with different numbers of degrees of freedom. Effectively, the merit function in the form

$$\begin{equation} L\left(\bm p\right)=\frac{1}{N-n_p}\sum _{i-1}^{N}{{\left(\frac{d_i-x_{p,i}}{\sigma _i}\right)}^2} \end{equation} \tag{ 16 }$$

is the mean distance between data and simulations in the measurement N-space, normalized by the number of degrees of freedom and by the uncertainties of the measurements.

The simulations return count rates of NISHe atoms, while the observations are total counts accumulated during the select ISM flow observation times. We had to make the two quantities compatible, and the choice was either to convert the model count rate into total counts or to convert the counts into the average count rate by dividing the counts and their errors by the duration of the select ISM flow observation times. We decided to adopt the second solution for practical reasons: recalculation could be done only once, and another selection would require converting all of the simulation cases, adding an unnecessary computational burden.

Although extensive pre-launch calibrations were conducted on the sensor (Fuselier et al. 2009; Möbius et al. 2012), we chose to avoid possible systematic changes in the observation conditions due to changes in the instrument functions from year to year by comparing observations and simulations separately for the 2009 and 2010 seasons.

From the simulations, we knew that count rate profiles as a function of spin angle should be Gaussian. For each season we selected a reference orbit and fitted its data with a Gaussian function specified in Equation (8). The fitted peak height f₀ from the reference orbit was used as a scaling factor for all the data points from that season. Effectively, this returned observed count rates relative to the fitted peak value at the reference orbit. To facilitate comparison, simulated count rates were scaled using a similar procedure. As the reference orbits we selected those with the highest count rates: orbit 16 in 2009 and orbit 64 in 2010.

Having brought simulations and observations to the common scale, we could look for the minimum of the four-parameter function given by Equation (16). This function is purely numerical, because simulation results $x_i\left(\bm p\right)$ used to evaluate the merit function L are purely numerical. In such a situation, calculating derivatives in the parameter space is problematic. Hence, the numerical method used to minimize this function used finite differences instead of derivatives of the merit function in the parameter space. In addition, since individual simulations are very time consuming, the number of evaluations of the merit function had to be kept as low as possible. We decided to adapt the quasi-Newton method suggested by Press et al. (2007).

In this method, a starting point $\bm p_{\mathrm{start}}$ in the parameter space is adopted and then, using finite differences to approximate partial derivatives, the gradient of the merit function is calculated. The gradient provides the direction of the slope of the merit function in the parameter space. This direction is sampled by calculation of a number of points $\bm p_1,\ldots,\bm p_n$, and a minimum at a point $\bm p_{\mathrm{loc}}$ along this direction is estimated from a parabola fit to the sampled points. The $\bm p_{\mathrm{loc}}$ that is found by this method becomes a new starting point. This iteration continues until the change in the magnitude of the function being minimized is considered sufficiently small. The final parameter set is then retrieved as $\bm p_{\mathrm{best}}$. It is essential in this method to specify initial steps in all coordinates of the parameter space such that the magnitudes of partial derivatives are similar.

Following this method, we computed χ² (merit function) values defined in Equation (16) for various sets $\bm p$ of flow direction longitude λ, latitude β, velocity v, and temperature T, starting from the values reported by Witte (2004). Resulting from the simulations were five-dimensional "landscapes" of χ²(λ, β, v, T). To find the best-fitting set of parameters $\bm p_{\mathrm{best}}$, we performed numerical minimizations of the χ²(λ, β, v, T) function separately on observations from each of the two seasons and collectively on combined observations from both seasons.

The minimizations were performed on a mesh of flow longitudes λ, with the longitude fixed and latitude, speed, and temperature being free parameters. The values of the three free parameters in the simulations were selected by the algorithm based on the local gradient of the minimized function and hence they are not on a regular mesh. The robustness of the minima of χ² values for all the λ mesh points was checked by restarting the minimization algorithm from the parameter set the algorithm had reported as optimum, until no further improvements could be obtained. Usually, such a restart of the procedure did not result in a significant reduction of χ².

After the absolute solution $\bm p_{\mathrm{best}}$ was found, we extended the simulations on a regular grid of parameters centered at the best solution $\bm p_{\mathrm{best}}$ in order to check on the covariance of parameters and to better illustrate the acceptable parameter region. This mesh regularization resulted in a very small correction of the best-fitting parameter set. The improvement in χ² value was only ∼0.001 with the changes in parameters of ∼50 K in temperature and ∼0.1 km s⁻¹ in the flow velocity.

During the minimization process we performed simulations using a total of about 4000 sets of the NISHe gas parameters. The minimization algorithm kept track of its own steps, so that results of simulations for individual parameter sets could be used in as many minimization processes as needed.

The values of the merit function χ² minimized for the mesh values of the flow direction λ are shown in Figure 17. The results of the minimization performed separately on the data from 2009 and 2010 (red and green, respectively) are consistent with each other and with the results of the minimization performed collectively on the data from 2009 and 2010 (blue). They are presented and discussed in greater detail in the following section.

Zoom In Zoom Out Reset image size