DETERMINATION OF INTERSTELLAR O PARAMETERS USING THE FIRST TWO YEARS OF DATA FROM THE INTERSTELLAR BOUNDARY EXPLORER

Interstellar neutral (ISN) flow measurements made by the Interstellar Boundary Explorer (IBEX) mission (McComas et al. 2009b) include direct H, He, O (Möbius et al. 2009), and D (Rodríguez Moreno et al. 2013, 2014) and a determination of the local interstellar cloud Ne/O ratio (Bochsler et al. 2012; Park et al. 2014). Each of the ISN species (e.g., H, D, He, O, Ne, etc.) flows through the outer heliospheric interfaces. Neutral particles pass through the outer heliosheath (beyond the heliopause but inside the bow shock/wave) and can charge-exchange with compressed interstellar plasma leading to the subsequent removal (or filtration) of ISN atoms, and the addition of secondary neutral atoms from this plasma. The so-called primary component of the ISN atoms is comprised of the atoms that flow directly through the heliosphere. The secondary component is comprised of atoms that have interacted in the heliosheath. Primary components are only modified by filtration due to loss through ionization (charge-exchange with the plasma, photo-ionization, electron impact ionization), gravitational effects, and the primary component of H is also acted upon by radiation pressure. Therefore, each neutral species' primary component provides key information about the local interstellar velocity distribution, which is analyzed to yield the best available estimates of the local interstellar bulk flow speed, flow direction, and the temperature of the species.

Studies of interstellar parameters including temperature, flow speed, and flow direction typically rely on measurements of neutral particles, pickup ions, or Lyα emissions. Previous studies have derived considerable information by examining interstellar He and interstellar H. Abundances have been derived for ³He, ⁴He, ¹⁵N, ¹⁴N, ¹⁸O, ¹⁶O, ²²Ne, ²⁰Ne, and other species. Observations by the PLAsma and SupraThermal Ion Composition instrument on board the Solar TErrestrial RElations Observatory mission (STEREO) of He⁺, Ne⁺, and O⁺ pickup ions have been used to establish the interstellar inflow direction of these species (Drews et al. 2010, 2012). However, there are very few studies of the interstellar flow speed, flow direction, and interstellar temperature for species other than H and He based on observations of ISN atoms.

In contrast to the primary components, the secondary neutral components are created by charge-exchange between ISN atoms and the plasma in the heliosheath. For example, the "warm breeze" of interstellar He (Kubiak et al. 2014) discovered by IBEX is both substantially hotter and slower than the primary neutral He from the local interstellar medium (LISM). Similarly, Möbius et al. (2009) and Park et al. (2015) found the imprint of the secondary O component in neutral atom maps.

Interstellar atoms experience a number of effects on their trajectory through the heliosphere. All ISN atoms experience ionization, predominantly through photo-ionization loss and charge-exchange collisions (Bzowski et al. 2013a). These ionization processes cause the loss of primary neutral atoms and the creation of new pickup ions. Interstellar H atoms also experience a large force associated with radiation pressure due to resonant absorption and re-emission of Lyα. This force from radiation pressure is roughly comparable in magnitude but opposite in direction to the force of gravity (Bzowski et al. 2013b, p. 67). Deflection of the primary ISN H flow by solar radiation pressure was revealed by IBEX from in situ observations for the first time (Schwadron et al. 2013; Katushkina et al. 2015).

The measurements of ISN He provide the best existing method for characterizing the flow properties of the LISM. Compared to other ISN species with lower first ionization potentials, interstellar He atoms are relatively unaffected by charge-exchange. In addition, due to its high universal abundance (second only to H), the primary component of ISN He has a relatively high flux at 1 au, and therefore provides for statistically accurate determinations of the LISM neutral temperature and bulk velocity via three distinct observational methods: backscattering of solar EUV (Meier & Weller 1972; Lallement et al. 2004), analysis of the related pickup ion distributions (Moebius et al. 1985, 1995; Gloeckler et al. 2004), and most recently direct neutral flow imaging (Witte et al. 1996, 2004; Möbius et al. 2009).

IBEX observations are remarkably sensitive to ISN He atoms, with a signal to background ratio of >1000, enabling in-depth study of ISN He flow characteristics (Bzowski et al. 2012, 2015; Möbius et al. 2012, 2015a; McComas et al. 2012, 2015a, 2015b; Leonard et al. 2015; Schwadron et al. 2015). These recent observations provide the most detailed and accurate direct measurement of the ISN flow vector and temperature to date. In addition, these observations illuminate departures from the idealized Maxwell–Boltzmann distribution (Kubiak et al. 2014; Sokół et al. 2015).

IBEX observations pose significant analytical challenges (Leonard et al. 2015; McComas et al. 2015a; Möbius et al. 2015b). Due to the fact that IBEX observes the ISN flow in the ecliptic plane over a limited longitude range, the ISN flow properties are constrained to a four-dimensional (4D) tube of parameters including the flow longitude ${\lambda }_{\mathrm{ISN}\infty }$, latitude ${\beta }_{\mathrm{ISN}\infty }$, speed ${V}_{\mathrm{ISN}\infty }$, and temperature ${T}_{\mathrm{ISN}\infty }$ (Bzowski et al. 2012; Lee et al. 2012; McComas et al. 2012; Möbius et al. 2012). While the allowable range of parameters included the previously established ISN flow vector determined from Ulysses measurements (Möbius et al. 2004; Witte 2004; Witte et al. 2004), the interstellar temperature from IBEX measurements for the same ISN flow vector was much higher than obtained previously (Bzowski et al. 2012; Möbius et al. 2012, 2015b).

Schwadron et al. (2015) used the three-step method to find the interstellar parameters. The three-step method consists of (1) finding the ISN He peak rate in ecliptic longitude that determines uniquely a relation (as part of the tube in parameter space) between the longitude ${\lambda }_{\mathrm{ISN}\infty }$ and the speed ${V}_{\mathrm{ISN}\infty }$ of the He ISN flow at infinity; (2) comparing the ISN He peak latitude (on the great circle swept out in each spin) to simulations, and thereby deriving unique values for ${\lambda }_{\mathrm{ISN}\infty }$ and ${V}_{\mathrm{ISN}\infty }$ along the parameter tube; and (3) finding the angular width of the He flow distributions as a function of spin-phase to derive the interstellar He temperature. For simulated peak latitudes, Schwadron et al. (2015) used a relatively new analytical tool that traces He atoms from beyond the termination shock into the position of IBEX and incorporates the detailed response function of IBEX-Lo. By varying interstellar parameters along the IBEX parameter tube, Schwadron et al. (2015) were able to find the specific parameters that minimize the χ² difference between observations and simulations. The new computational tool for simulating neutral atoms through the integrated IBEX-Lo response function makes no assumptions or expansions with respect to spin axis pointing or frame of reference. Thus, Schwadron et al. (2015) were able to move beyond closed-form approximations and utilized observations of interstellar He during the complete five-year period from 2009 to 2013 when the primary component of interstellar He is most prominent. This three-step method resulted in a He ISN flow longitude 756 ± 14, latitude −512 ± 027, speed 25.4 ± 1.1 km s⁻¹, and temperature 8000 ± 1300 K, where the uncertainties are related and apply along the IBEX parameter tube (Schwadron et al. 2015).

The parameters derived by Schwadron et al. (2015) are compared to those derived by other studies in Table 1. In particular, we note the similarity in the ISN flow longitude, speed, and latitude from IBEX (Bzowski et al. 2015; Leonard et al. 2015; McComas et al. 2015a, 2015b; Schwadron et al. 2015) compared to those three parameters derived from Ulysses data (Witte et al. 2004; Bzowski et al. 2014; Wood et al. 2015). The key difference is that the temperatures derived from IBEX measurements are significantly higher than those derived from early estimates by Ulysses (e.g., Witte et al. 2004). The re-analyses of Ulysses observations (Bzowski et al. 2014; Wood et al. 2015) have lead to significantly higher temperature estimates from the Ulysses data, which are closer to but still noticeably smaller temperatures than the IBEX results.

Table 1.  ISN Flow Parameters Using Direct ISN He Neutral Flow Observations by Either the Ulysses or IBEX Spacecraft. (J2000 Coordinates Used throughout). (from Schwadron et al. 2015a)

Publication	${\lambda }_{\mathrm{ISN}\infty }$ (°)	${V}_{\mathrm{ISN}\infty }$ (km s−1)	${\beta }_{\mathrm{ISN}\infty }$ (°)	TISN (kK)	Spacecraft
Witte et al. (2004)	75.4 ± 0.5	26.3 ± 0.4	−5.2 ± 0.2	6.30 ± 0.34	Ulysses
Bzowski et al. (2014)	75.3 + 1.2(−1.1)	26.0 + 1.0(−1.5)	−6.0 ± 1.0	7.5 + 1.5(−2.0)	Ulysses
Wood et al. (2015)	75.54 ± 0.19	26.08 ± 0.21	−5.44 ± 0.24	7.26 ± 0.27	Ulysses
Leonard et al. (2015)a
(z ∼ 0, 2012–14)	74.5 ± 1.7	27.0 + 1.4(−1.3)	−5.2 ± 0.3		IBEX
McComas et al. (2015a)	∼75	∼26	∼−5	7−9.5	IBEX
Bzowski et al. (2015)a
(z ∼ 0, 2012–14)	75.3 ± 0.6	26.7 ± 0.5	−5.14 ± 0.16	8.15 ± 0.39	IBEX
Bzowski et al. (2015)a
(z, no restriction, 2009–14)	75.8 ± 0.5	25.8 ± 0.4	−5.17 ± 0.10	7.44 ± 0.26	IBEX
Schwadron et al. (2015a)a
(z ∼ 0, 2012–14)	75.8 ± 1.8	25.4 ± 1.4	−5.11 ± 0.28	7.9 ± 1.4	IBEX
Schwadron et al. (2015a)a
(z, no restriction, 2009–13)	75.6 ± 1.4	25.4 ± 1.1	−5.12 ± 0.27	8.0 ± 1.3	IBEX

Note.

^aThe total uncertainties lie along the parameter tube and are therefore dependent on one another.

Download table as:  ASCIITypeset image

The present study derives ISN flow properties from primary interstellar O atoms observed by IBEX. By observing a different primary species than He and H, we have a way to further test whether the local interstellar cloud (LIC, or more specifically, the part of the LIC that surrounds the solar system) is isothermal. These measurements provide direct insight into the nature of turbulence in the interstellar plasma that drive the state of the plasma away from a thermal equilibrium. This is of particular interest given that the solar system is near the edge of the LIC (Frisch et al. 2009) where turbulent interactions may be more significant than in other regions of the LIC.

We apply the same numerical tools as those developed by Schwadron et al. (2015), with the key difference being that we treat interstellar O instead of He. We are able to derive the combined interstellar O density times the ionization rate (referenced to 1 au). The model takes into account the effects of ionization and the survival probability through the heliosphere, inside the termination shock.

Understanding the local interstellar abundance of oxygen is also important for understanding nucleosynthesis and the oxygen content of interstellar dust grains. Specifically, oxygen and neon are mainly produced in massive stars, which have relatively short lifetimes. The Milky Way should have been rapidly enriched in oxygen and neon within the first few stellar generations. Further, neon and oxygen are produced from similar sources. This suggests that if oxygen remains predominantly within the interstellar medium as a gas (e.g., not locked away in grains), the observed neon-to-oxygen ratio in the interstellar gas should be similar to that in solar abundances. Any difference between the solar and interstellar ratios of neon to oxygen would then trace back to the poorly known oxygen content of interstellar dust grains.

Interstellar dust permeates space and flows through the heliosphere as part of the interstellar wind (Grün et al. 1993; Frisch et al. 1999; Sterken et al. 2015). Oxygen is locked up in silicates that contain a significant fraction of the interstellar dust mass, and possibly all of the dust mass for the interstellar cloud surrounding the heliosphere (Slavin & Frisch 2008; Altobelli et al. 2016). Levels of oxygen depletion onto dust grains are difficult to ascertain for low-density interstellar sightlines because the O i 1302 Angstrom absorption feature is heavily saturated and gives uncertain oxygen abundances. Measurements of O/Ne and O/He in the interstellar wind through the heliosphere provide insight into the interstellar O abundances in the gas, and therefore the depletion of O onto dust grains since Ne and He remain in gas phase. Three substantial corrections are required to obtain interstellar abundances from the IBEX O, He, and Ne data, first for the partial ionization of these elements in the interstellar gas (Slavin & Frisch 2008), second for the filtration of these elements through the heliosheath regions (Müller & Zank 2004), and finally the survival of these neutrals propagating from the termination shock to 1 au (Bochsler et al. 2012; Park et al. 2014). After correcting for these effects, the IBEX O/Ne data indicate significant amounts of oxygen depletion onto dust grains.

Direct study of interstellar oxygen abundances with IBEX have been made in several papers (Bochsler et al. 2012; Park et al. 2014). The question of oxygen abundances requires knowledge of LISM parameters, since they help determine the filtration of oxygen in the heliosheath and the ionization of oxygen in the heliosphere, inside the termination shock. Therefore, the present study into interstellar O properties should have an impact on our understanding of the interstellar O abundance.

The paper is organized as follows. Section 2 discusses the observations utilized for the study. Section 3 details our approach to analysis. Section 4 summarizes the results of the oxygen χ² analysis and Section 5 summarizes the results of the maximum likelihood analysis. Section 6 discusses the derivation of the interstellar O density. Section 7 discusses the broader implications of the interstellar oxygen results in the context of previous work. We summarize the main conclusions in Section 8. The paper also includes several appendices. Appendix A describes the geometric factors used in the analysis. Appendix B describes our general approach to statistical and propagation uncertainties. New analytical methods are detailed for minimizing the χ² (Appendix C) and maximizing likelihood (Appendix D) in solving for a model's scaling factor and the background rate within each orbit. Also detailed are propagation and statistical uncertainties from χ² fitting (Appendix E) and maximum likelihood fitting (Appendix F) in nonlinear models.

IBEX has two energetic neutral atom (ENA) sensors for remotely mapping the global heliosphere and making direct measurements of ISN atoms (McComas et al. 2009b). The IBEX-Lo sensor measures neutral atoms from ∼10 eV to 2 keV and includes a triple coincidence time-of-flight analysis to provide compositional information (Fuselier et al. 2009b; Möbius et al. 2009). The IBEX-Hi sensor measures ENAs from ∼300 eV to 6 keV (Funsten et al. 2009a).

IBEX is a Sun-pointed spinner with the sensor field of view pointing at 90° from the spin axis. The IBEX-Lo sensor sweeps out a great circle on the celestial sphere roughly every 15 s. During the season of prime interstellar O viewing in the spring of each year the Earth and thus IBEX ram into the oncoming ISN flow, which covers a limited spin-phase range close to the ecliptic. The ISN O flow rate peaks around February 8 each year.

Our analysis takes advantage of the changes through the spacecraft longitude in the oxygen rate distributions versus spin-phase. These spin-phase distributions provide the fundamental measurement that can be compared to simulations to derive unique values for the flow parameters. The observed rate distributions are accumulated over very low background periods ("good times") for each orbit in 2009 and 2010 when O can be measured well. These good times are the same as those for He documented by Möbius et al. (2012).

The IBEX-Lo entrance system accepts incoming neutral atoms through a large-area collimator with a field of view of 7° FWHM. After passing through the collimator, neutrals collide with a conversion surface where a small fraction of these incoming atoms are converted into negative ions. The negative ions are then filtered based on their energy and charge by an electrostatic analyzer (ESA). After post-acceleration to boost their energy, negative ions pass through a time-of-flight system, which, together with the energy and charge measurements, determines the mass and therefore the atomic species of these particles.

The conversion surface acts differently for different atomic species. Incoming O atoms predominantly convert into negative ions when they collide with the IBEX-Lo diamond-like conversion surface. During optimal ISN O observing periods near the beginning of each year, the motion of IBEX, which moves with Earth around the Sun at ∼30 km s⁻¹, opposes the velocity of incident neutral atoms. ISN O atoms, based on IBEX-Lo observations, move at an average speed of ∼20–30 km s⁻¹ relative to the Sun in the outer heliosphere. The O atoms that make it in to 1 au increase in kinetic energy and speed to ∼50 km s⁻¹ in the solar inertial frame due to the Sun's gravitational attraction. During the IBEX-Lo O observing periods, in the frame of the spacecraft, incident ISN O atoms have typical speeds of ∼80 km s⁻¹ (kinetic energy ∼530 eV) with respect to the IBEX spacecraft. Upon conversion, O⁻ atoms lose ∼10%–15% of their incident energy.

The converted O⁻ atoms during optimal O observing periods have energies ∼439 eV, the central energy of step 6 of the IBEX-Lo ESA. While the ISN O temperature slightly broadens the angular distribution at 1 au, the incoming ISN O distribution is still narrow and beam-like. The IBEX-Lo ESA steps admit a broad range of energies (ΔE/E ∼ 0.7), so the vast majority of the converted O⁻ atoms fall within ESA step 6. This energy signature provides a straightforward identification of ISN O in IBEX observations (Möbius et al. 2012). It should be noted that there is a small amount of Ne that contributes to the O distributions reported here. As detailed in Section 3, this Ne mixing is accounted for by using an effective mass for the combined O and Ne distribution.

The geometric factor calculation derived from data taken from calibration runs is detailed in Appendix A. The IBEX-Lo geometric factor for heavy neutral atoms, such as O, is controlled by the combination of surface conversion and sputtering to negative ions. The calculation takes into account the energy response curve of IBEX-Lo, how its peak energy varies with the incoming energy of the neutral atoms, and the absolute geometric factor as a function neutral atom energy.

There are five parameters that can be determined through analysis of IBEX interstellar O neutrals: the bulk flow speed, ${V}_{\mathrm{ISNO}\infty }$, the temperature, ${T}_{\mathrm{ISNO}\infty }$, the flow longitude, ${\lambda }_{\mathrm{ISNO}\infty }$, the flow latitude, ${\beta }_{\mathrm{ISNO}\infty }$, and the survival probability times the flow density, ${n}_{\mathrm{ISNO}\infty }$ × S_p. The survival probability S_p of O atoms from the termination shock to the point of observation at 1 au depends primarily on the ionization rate ${\beta }^{\mathrm{ion}}={\beta }_{1}^{\mathrm{ion}}({R}_{1}^{2}/{r}^{2})$, which is referenced at R₁ = 1 au by ${\beta }_{1}^{\mathrm{ion}}$ and scales as 1/r² where r is radial distance from the Sun. This scaling with radial distance neglects electron impact ionization, which makes a small contribution.

As detailed in Appendices C and D, in the fitting performed we solve for a scaling constant, A_k, in each interstellar season, k. The scaling constant is proportional to ${n}_{\mathrm{ISNO}\infty }$ × S_p. The survival probability is related to the ionization rate at 1 au,

$$\begin{eqnarray}&&{S}_{p}=\exp \left(\displaystyle \frac{{{mR}}_{1}^{2}{\beta }_{1}^{\mathrm{ion}}}{l}[\theta -{\theta }_{\infty }]\right),\end{eqnarray} \tag{ 1 }$$

where θ is the position angle at the time of observation and ${\theta }_{\infty }$ is the position angle when the neutral was far from the Sun (at $t=-\infty$), as defined by Lee et al. (2012). The angular momentum is given by l. The model integrates precisely along neutral trajectories, accurately taking into account the variation in the gravitational force, ionization rate, and survival probability with distance from the Sun. The model remains linearly proportional to the scaling constant A_k that varies with the ionization rate at 1 au and the interstellar density near the termination shock. We discuss how the scaling constant A_k is solved in Appendices C and D.

All parameters other than ${\beta }_{1}^{\mathrm{ion}}$ apply at or near the termination shock since we neglect charge-exchange processes that occur within the heliosheath. In particular, filtration of O due to charge-exchange causes the ionization of neutral O atoms predominantly in the outer heliosheath, constituting a loss process for neutrals that modifies the neutral distribution function. Therefore, it is important to take filtration into account when determining the interstellar parameters in the LISM beyond the heliosphere.

The procedure for finding interstellar parameters starts with a comparison between observed rate distributions as a function of spin-phase and corresponding modeled rate distributions. The rate distributions are accumulated over the good times for a given orbit. The modeled rate distributions are also determined at a series of points in time separated by a maximum of 0.8 days (this time interval was chosen as a compromise between efficiency and convergence; smaller time intervals gave almost precisely the same answer as the chosen 0.8 day interval). The modeled distributions are then compared to observed rate distributions to yield a χ² difference (Section 4) or a likelihood (Section 5). We vary interstellar parameters to minimize the χ² and maximize the likelihood, therefore finding the parameters that optimize the fits to observations.

Appendix B details our general approach to fitting and solving for uncertainties in fit parameters. Appendix C derives an analytical method for minimizing the χ² and thereby solving for background rates in each orbit and solving for a scaling factor, the interstellar O density at the termination shock times the ionization rate at a reference distance (1 au). This scaling factor is treated as an unknown with which all modeled rates scale linearly.

The other four parameters used in the fit are the bulk speed, temperature, longitude, and latitude of the interstellar O flow (${V}_{\mathrm{ISNO}\infty }$, ${T}_{\mathrm{ISNO}\infty }$, ${\lambda }_{\mathrm{ISNO}\infty }$, and ${\beta }_{\mathrm{ISNO}\infty }$). In searching these four parameter values to minimize the χ² or maximize likelihood, it is important to recognize that there is an underlying dependency between the parameters that influences the 4D dependence of the χ² and likelihood functions. These interdependencies between parameters should be understood in order to do an accurate and efficient search for the global χ² minimum.

The hyperbolic trajectory equation, along with the IBEX observation at perihelion, leads to a strict relationship between the speed at infinity ${V}_{\mathrm{ISN}\infty }$ and the flow angle ${\lambda }_{\mathrm{ISN}\infty }$ in ecliptic longitude. The relation is governed by the motion of neutrals in the Sun's gravity, which has one parameter, the ecliptic longitude λ_Peak where the interstellar bulk flow hits perihelion (referred to hereafter as the peak longitude). We assume a small angle between the ecliptic plane and the trajectory plane (Lee et al. 2012), which is appropriate given that the upstream direction of the interstellar flow is ∼5° above the ecliptic based on measurements of neutral He. The bulk flow velocity is defined as follows:

$$\begin{eqnarray}&&{V}_{\mathrm{ISN}\infty }=\sqrt{\displaystyle \frac{{{GM}}_{s}}{{R}_{1}}\left(\displaystyle \frac{-1}{\cos ({\lambda }_{\mathrm{ISN}\infty }+180^\circ -{\lambda }_{\mathrm{Peak}})}-1\right)}.\end{eqnarray} \tag{ 2 }$$

where R₁ denotes 1 au and M_s denotes the mass of the Sun. Equation (2) is strictly exact if the trajectory plane is identical to the ecliptic plane. Equation (2) applies for all species provided that they are acted upon only by gravity and ionization. (Radiation pressure can also be factored in; we have neglected it here for simplicity since we are primarily focused on interstellar O, which is not strongly affected by radiation pressure.)

One of the remarkable properties of the trajectory relation between interstellar speed and the peak longitude is that it naturally allows for a separation between primary and secondary components of ISN atoms. The charge-exchange process in the heliosheath renders secondary interstellar atoms much slower than the primary neutral atoms that come directly from the interstellar medium. In the case of interstellar O, primary atoms should have flow speed of ∼26 km s⁻¹, similar to the flow speed of primary He. However, the secondary component is slowed and heated relative to the primary component. The average flow speed of secondary O has not been estimated.

The flow of O is analogous in some ways to He. The primary neutral component of He has a flow speed of 25.4 ± 1.1 km s⁻¹ while the secondary component has a flow speed of ∼11.3 km s⁻¹. The large speed differential shifts the peak longitude by ∼26°. The Sun acts as a "gravitational lens" that separates the primary from secondary components. This same effect leads to the separation of the secondary O component from the primary component (Park et al. 2015). Further, the secondary O neutrals come from the hot plasma flow in the outer heliosheath that is highly asymmetric relative to the primary flow direction (see Kubiak et al. 2016). This leads to a significant observed offset of the apparent secondary flow at the heliospheric boundary in both flow longitude and flow latitude. Finally, the O secondary component is much better separated from the O primary flow than for He and H because O has a much lower thermal speed (due to its higher atomic mass, the thermal speed of O is ∼1/2 the thermal speed of He and ∼1/4 the thermal speed of H). The lower thermal speed of O creates a much more sharply peaked velocity distribution function than for the He and H species.

The angular width of the ISN flow distribution as a function of spin-angle ψ observed at IBEX is controlled solely by the thermal speed (expressed through the temperature) and the ISN bulk flow speed ${V}_{\mathrm{ISN}\infty }$ at infinity. Consequently, the width σ_ψ of the spin-angle distribution at the location of the bulk flow intercept at 1 au is defined by

$$\begin{eqnarray}&&{\sigma }_{\psi }^{2}=\displaystyle \frac{{k}_{B}{T}_{\mathrm{ISN}}}{{{mV}}_{\mathrm{ISN}\infty }^{2}}\displaystyle \frac{{(1+{w}_{\mathrm{ISN}\infty }^{2})}^{2}}{{w}_{\mathrm{ISN}\infty }^{2}(2+{w}_{\mathrm{ISN}\infty }^{2})},\end{eqnarray} \tag{ 3 }$$

where T_ISN is the temperature of the ISN and m is the mass. Here ${w}_{\mathrm{ISN}\infty }={V}_{\mathrm{ISN}\infty }/{V}_{E}$ is the dimensionless ISN flow speed normalized to the average speed of the Earth at 1 au, ${V}_{E}=\sqrt{{{GM}}_{s}/{R}_{1}}$ (Lee et al. 2012).

Note that the computation of the ISN temperature from the width of the observed flow distributions requires the mass of the species. The heavy ISN flow distribution is mostly O but also contains a sizable fraction of Ne (Bochsler et al. 2012; Park et al. 2014). Using the Ne/O ratio derived by Park et al. (2014) and factoring in the IBEX-Lo efficiencies for O and Ne from calibration, the observed Ne/O ratio has been translated into an effective mass of m_O = 16.85 ± 0.3 for the combined O and Ne distribution. This effective mass value was used for the determination of the O temperature, which represents a combined O and Ne temperature.

We begin the process of searching for interstellar parameters by choosing a guess for ${T}_{\mathrm{ISNO}\infty :0}$, ${\lambda }_{\mathrm{ISNO}\infty :0}$, and ${\beta }_{\mathrm{ISNO}\infty :0}$. We then vary ${V}_{\mathrm{ISNO}\infty }$ and find the minimum χ² in this dimension, ${V}_{\mathrm{ISNO}\infty :1}$. Our initial guess for ${T}_{\mathrm{ISNO}\infty :0}$, ${\lambda }_{\mathrm{ISNO}\infty :0}$, and ${\beta }_{\mathrm{ISNO}\infty :0}$ is based on the results for interstellar He (Schwadron et al. 2015). With the initial minimum for ${V}_{\mathrm{ISNO}\infty :1}$, and the guesses for ${\lambda }_{\mathrm{ISNO}\infty :0}$ and ${\beta }_{\mathrm{ISNO}\infty :0}$, we vary ${T}_{\mathrm{ISNO}\infty }$ to find the χ² minimum or maximum likelihood in this dimension, ${T}_{\mathrm{ISNO}\infty :1}$.

The next step is to find the χ² minimum or maximum likelihood while varying ${\lambda }_{\mathrm{ISNO}\infty }$. The relations previously specified (2) and (3) provide the means to consistently vary ${V}_{\mathrm{ISNO}\infty }$(${\lambda }_{\mathrm{ISNO}\infty }$) and ${T}_{\mathrm{ISNO}\infty }$(${\lambda }_{\mathrm{ISNO}\infty }$) as a function of the flow longitude. We specify the peak longitude ${\lambda }_{\mathrm{Peak}:1}$ based on the initial χ² minimum for bulk speed:

$$\begin{eqnarray}{\lambda }_{\mathrm{Peak}:1} & = & {\lambda }_{\mathrm{ISNO}\infty :1}+180^\circ -{\cos }^{-1}\\ & & \times \left(\displaystyle \frac{-1}{1+{R}_{1}{({V}_{\mathrm{ISNO}\infty :1})}^{2}/({{GM}}_{S})}\right).\end{eqnarray} \tag{ 4 }$$

With this value for the peak longitude, Equation (2) provides specification for determining ${V}_{\mathrm{ISNO}\infty }$(${\lambda }_{\mathrm{ISNO}\infty }$). Solutions are highly degenerate along the parameter tube. In other words, variation in longitude ${\lambda }_{\mathrm{ISNO}\infty }$ with covariation in ${V}_{\mathrm{ISNO}\infty }$(${\lambda }_{\mathrm{ISNO}\infty }$) and ${T}_{\mathrm{ISNO}\infty }$(${\lambda }_{\mathrm{ISNO}\infty }$) leads to only small changes in the χ² or likelihood. Therefore, the search for a global minimum in the χ² or a global maximum in likelihood most strongly depends on where in longitude along the parameter tube that solutions best match observations.

We use Equation (3) to provide a reference angular width ${\sigma }_{\psi :1}$ based on the initial determinations of χ² minima from ${V}_{\mathrm{ISNO}\infty :1}$ and ${T}_{\mathrm{ISNO}\infty :1}$. The temperature as a function of ${\lambda }_{\mathrm{ISNO}\infty }$ follows

$$\begin{eqnarray}{T}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISNO}\infty }) & = & \displaystyle \frac{{m}_{{\rm{O}}}{\sigma }_{\psi :1}^{2}}{{k}_{B}}{[{V}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISNO}\infty })]}^{4}\\ & & \times \displaystyle \frac{2{V}_{E}^{2}+{[{V}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISN}\infty })]}^{2}}{\{{V}_{E}^{2}+{[{V}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISN}\infty })]}^{2}\}{}^{2}}.\end{eqnarray} \tag{ 5 }$$

By varying the interstellar temperature and interstellar speed as indicated here we break the underlying degeneracy in the parameter dependence. This greatly increases the efficiency of finding the global χ² minimum or the maximum likelihood.

When we find a χ² minimum or a maximum likelihood in the search in interstellar longitude, ${\lambda }_{\mathrm{ISN}\infty :2}$, we update the estimates for the interstellar speed and temperature: ${V}_{\mathrm{ISNO}\infty :2}={V}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISN}\infty :2})$ and ${T}_{\mathrm{ISNO}\infty :2}={T}_{\mathrm{ISNO}\infty }({\lambda }_{\mathrm{ISN}\infty :2})$. Using these values for the longitude, ${\lambda }_{\mathrm{ISN}\infty :2}$, the speed, ${V}_{\mathrm{ISNO}\infty :2}$, and the temperature, ${T}_{\mathrm{ISNO}\infty :2}$, we perform a search in latitude to find the corresponding χ² minimum or maximum likelihood at ${\beta }_{\mathrm{ISNO}\infty :2}$.

As a final step in the analysis, we vary all parameters about the optimal fit. This serves both to guarantee that we have found an optimal fit, and allows us to map out the multidimensional shape of the χ² and likelihood functions of interstellar parameters. From the shape of χ² and likelihood functions, we derive curvatures in various parameter dimensions, from which fit and propagation uncertainties are determined (see Appendices B, E, and F).

The method described above represents an iterative process whereby initial guesses for the four interstellar parameters, the scale factor (which determines the interstellar density and the ionization rate at 1 au) and background rates are successfully improved. This scheme converges rapidly. In our case we found convergence to less than 5% of uncertainties in three iterations.

The behavior of the χ² function and the maximum likelihood function are also used to ascertain uncertainties in the derived parameters. The uncertainty in the scaling factor δA_k is detailed in Appendices C and D. We describe here the approach to finding the uncertainties in the other four interstellar parameters $({V}_{\mathrm{ISNO}\infty },{T}_{\mathrm{ISNO}\infty },{\lambda }_{\mathrm{ISNO}\infty },{\beta }_{\mathrm{ISNO}\infty })$. The implicit parameters are the background rates ${B}_{i=1..N}$ where N = 12 is the number of orbits, and the linear scaling factors ${A}_{k=1..Q}$ where Q = 2 is the number of interstellar seasons. These parameters are implicit because, as detailed in Appendices C and D, in every search through parameter space we solve directly for the background rates and scaling factors.

After specifying a guess for interstellar parameters, we then perform a search through parameter space by sweeping first the interstellar speed and then the interstellar temperature. For example, we may represent the one-dimensional (1D) sweep in speed x as ${\chi }^{2}={\chi }^{2}(x,{T}_{\mathrm{ISNO}\infty },{\lambda }_{\mathrm{ISNO}\infty },{\beta }_{\mathrm{ISNO}\infty })$. This represents a 1D cut in parameter space. And because the interstellar longitude is frozen, this 1D cut determines the width of the parameter tube for the interstellar speed, $\delta {V}_{\mathrm{ISNO}\infty :1}$, based on the curvature of χ² dependence on interstellar speed (see Appendix E):

$$\begin{eqnarray}&&{\chi }^{2}(x)={\chi }_{\mathrm{est}}^{2}+{C}_{2}{(x-{V}_{\mathrm{ISNO}\infty :1})}^{2},\end{eqnarray} \tag{ 6 }$$

where C₂ is the curvature of the fit. Similarly, the curvature of likelihood function is used to derive the uncertainty $\delta {V}_{\mathrm{ISNO}\infty :1}$ (see Appendix F).

We approach uncertainties using the method outlined in Appendices E and F, which is an extension of the approach discussed by Schwadron et al. (2013). There are two independent sources of uncertainty (see Appendix B): (1) the propagation of measurement uncertainties associated with the observed rates, and (2) the statistical uncertainty associated with the χ² fit. For specificity, we consider the derivation of the uncertainty of the bulk speed, although a similar approach is employed also for the other interstellar parameters. The propagation uncertainty is found by inverting the curvature matrix (one half times the Hessian matrix, Press et al. 1992). The propagation uncertainty is given by

$$\begin{eqnarray}&&{(\delta {V}_{\mathrm{ISNO}\infty :1}^{\mathrm{pr}})}^{2}={C}_{2}^{-1}.\end{eqnarray} \tag{ 7 }$$

The statistical uncertainty is

$$\begin{eqnarray}&&{(\delta {V}_{\mathrm{ISNO}\infty :1}^{\mathrm{st}})}^{2}={(\delta {V}_{\mathrm{ISNO}\infty :1}^{\mathrm{pr}})}^{2}{\chi }_{\mathrm{red}}^{2},\end{eqnarray} \tag{ 8 }$$

where ${\chi }_{\mathrm{red}}^{2}={\chi }_{\mathrm{est}}^{2}/\nu$ is the reduced χ² at the global minimum ${\chi }_{\mathrm{est}}^{2}$. The quantity ν is the number of degrees of freedom: $\nu =N-{N}_{o}-\gamma$ where N is the total number of data points and N_o is the total number of orbits. In addition to the N_o background rates, we have γ = 6 additional parameters: the four interstellar parameters and two linear scaling factors (one scaling factor for each of the two interstellar seasons in 2009 and 2010). Combining the two forms of uncertainty, we have a total uncertainty given by

$$\begin{eqnarray}&&\delta {V}_{\mathrm{ISNO}\infty :1}=\delta {V}_{\mathrm{ISNO}\infty :1}^{\mathrm{pr}}\sqrt{1+{\chi }_{\mathrm{red}}^{2}}.\end{eqnarray} \tag{ 9 }$$

Note that the speed uncertainty $\delta {V}_{\mathrm{ISNO}\infty :1}$ applies across the parameter tube and combines the propagation and statistical uncertainties. There is also a covariant uncertainty along the parameter tube (in longitude), $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$, that is the largest source of uncertainty for the interstellar speed of oxygen. The 1σ upper and lower limits of the speed, ${V}_{\mathrm{ISNO}\infty :1}\pm \delta {V}_{\mathrm{ISNO}\infty :1}$, determine the corresponding upper and lower limits of the peak longitude, ${\lambda }_{\mathrm{Peak}\pm \sigma }$, through Equation (4). Analogous solutions are found for uncertainties based on the curvature of the inverse likelihood function (Appendix F).

A similar procedure is used to define the 1σ limits on the temperature parameter tube. We determine the χ² dependence or likelihood dependence on temperature, varying only the temperature and leaving the other three parameters fixed. The statistical and propagation uncertainties apply across the parameter tube, $\delta {T}_{\mathrm{ISNO}\infty }$. Another source of uncertainty for the temperature arises from its dependence on interstellar speed in Equation (5).

The statistical and propagation uncertainties in flow longitude $\delta {\lambda }_{\mathrm{ISNO}\infty }$ and flow latitude $\delta {\beta }_{\mathrm{ISNO}\infty }$ are also determined. However, as already discussed, the variation in longitude is done while also varying speed and temperature along the parameter tubes. Therefore, the uncertainty in longitude is understood to apply along the parameter tube.

The data used here are from the 2009 and 2010 interstellar seasons. The 2009 season consists of IBEX orbits 14–19, which occurred during January 16–February 28. The 2010 season consists of orbits 61–67, which occurred during January 10–March 3. The data for orbit 62 was lost due to an anomaly that caused the reset of the spacecraft power system.

Table 2 shows the O ISN velocity vector and temperature results of the χ² analysis in 2009, in 2010, and for the combined data in 2009 and 2010. Figures 1–2 show the reduced χ² values separately for the flow latitude, longitude, speed, and temperature using data from each season, and from the combined 2009–2010 seasons in Figure 3. Figures 4 and 5 compare modeled and observed rates in each of the orbits studied for the optimal fitting parameter values (the third column in Table 2). The observational uncertainties shown and included in the analysis arise from Poisson counting statistics. Note that the χ² analysis automatically excludes bins that have no counts, since the uncertainty in these bins diverges. As such, Figures 4 and 5 show only the bins that contain 1 or more counts.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Results of the ${\tilde{\chi }}^{2}$ Minimization

Parameter	2009	2010	2009–10
${V}_{\mathrm{ISNO}\infty }$ (km s−1)	26.2	26.0	26.1
$\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$ a	0.7	0.4	0.3
$\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$ a	0.6	0.4	0.3
$\delta {V}_{\mathrm{ISNO}\infty }$ a	0.9	0.6	0.5
$\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ b	4.6	3.4	2.7
${T}_{\mathrm{ISNO}\infty }$ (K)	7950	6500	7160
$\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$ a	1300	1200	900
$\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$ a	1200	1100	700
$\delta {T}_{\mathrm{ISNO}\infty }$ a	1800	1600	1150
$\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ b	3900	2400	2100
${\lambda }_{\mathrm{ISNO}\infty }$ (°)	79.7	74.0	76.0
$\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$ c	4.5	3.3	2.6
$\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{st}}$ c	4.2	3.0	2.2
$\delta {\lambda }_{\mathrm{ISNO}\infty }$ c	6.2	4.4	3.4
${\beta }_{\mathrm{ISNO}\infty }$ (°)	−4.4	−4.3	−4.4
$\delta {\beta }_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$ (°)	0.4	0.3	0.3
$\delta {\beta }_{\mathrm{ISNO}\infty }^{\mathrm{st}}$ (°)	0.4	0.3	0.2
$\delta {\beta }_{\mathrm{ISNO}\infty }$ (°)	0.5	0.5	0.3

Notes.

^aThe quantities $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, $\delta {V}_{\mathrm{ISNO}\infty }$, $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, and $\delta {T}_{\mathrm{ISNO}\infty }$, represent propagation, statistical, and total uncertainties across the parameter tube for interstellar O speed and temperature, respectively. ^bThe quantities $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ and $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ represent the uncertainties along the parameter tube covariant with the uncertainty in flow longitude. ^cThe quantities $\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, and $\delta {\lambda }_{\mathrm{ISNO}\infty }$ represent propagation, statistical, and total uncertainties of interstellar longitude with interstellar speed and temperature varied along the parameter tube.

Download table as:  ASCIITypeset image

Results show a global χ² minimum as a function of all parameters. The reduced χ² is lower than 1 (e.g., close to 0.7) when combining both the 2009 and 2010 interstellar seasons. This suggests that the observational uncertainties used in the analysis may be overestimates. It is also possible that because of nonlinearities in our model, we have underestimated the number of free parameters. We have used the standard relation $\sigma =R/\sqrt{{N}_{c}}$ for the uncertainty of the measured rate R. Here, N_c is the number of counts in a given bin. This relation is accurate for standard uncertainties, but only asymptotically accurate for large count rates in the case of Poisson-distributed uncertainties. Had uncertainties been ∼16% smaller, our reduced χ² would have been very close to 1. The overestimation of uncertainties may be associated with our lowest count rate bins, which have rates near the background levels.

The final parameters derived from χ² fits are the scaling constants A_k for 2009 and 2010. Recall that A_k represents a density times a survival probability. For the 2009 period (k = 1), we find A₁ = (9.2 ± 2.3) × 10⁻⁶ cm⁻³ and for the 2010 period (k = 2), ${A}_{2}=(8.9\pm 2.4)\times {10}^{-6}$ cm⁻³. We relate these scaling constants to an interstellar density in Section 6.

Based on the results of the χ² analysis, we derive the parameter tube for the interstellar temperature (Figure 6, top) and speed (Figure 6, bottom) as a function of interstellar longitude. The derived temperature parameter tube for interstellar O and interstellar He are almost identical. At the flow longitude of 76°, the O temperature is ∼600 K less than the He temperature. This constitutes less than a 10% difference.

Zoom In Zoom Out Reset image size

Figure 6 reveals that the derived interstellar speeds are statistically distinct at the 1σ level. For example, at a flow longitude of 76°, the interstellar O speed would need to be ∼0.9 km s⁻¹ higher than the interstellar He speed.

The low counting statistics observed in some of the spin-phase bins requires that we take into account the departures from Gaussian probability distributions due to Poisson-distributed probabilities. Appendices F and D describe our approach to maximum likelihood. Results of the analysis show that the relations associated with the interstellar O parameter tube remain almost unchanged from 2009 to 2010. Table 3 shows the results for Peak longitude in 2009, 2010, and for the combined two (2009–2010) interstellar seasons. Since the peak longitude specifies the parameter tube, it is apparent that the these parameter tubes are quite stable for the two interstellar seasons.

Table 3.  Peak Longitude Derived from Maximum Likelihood Fits

Season	Peak Long. (°)
2009	1308
2010	1325
2009–2010	${131.6}_{-0.9}^{+0.8}$

Download table as:  ASCIITypeset image

Table 4 shows the O ISN velocity vector and temperature results of the maximum likelihood analysis for the combined data in 2009 and 2010. Figure 7 shows the fits for the reduced inverse likelihood $\tilde{l}=-2\mathrm{ln}{\rm{\Gamma }}/(N-M)$ where Γ is the joint probability associated with the data set, N = 120 is the number of data points, and $M=12+4+2=18$ (12 background rates, 4 interstellar parameters, and 2 model amplitudes) is the number of free parameters. Note that the interstellar O speed (lower right panel) is varied across the parameter tube with interstellar O temperature varied according to Equation (5) to maintain a fixed width of the spin-angle distribution. The interstellar temperature (lower left panel) and the interstellar inflow latitude (upper right panel) are varied independently with respect to other interstellar parameters. The interstellar flow longitude (top left panel) is varied with the interstellar speed and interstellar temperature varied along the parameter tube according to Equations (2) and (5), respectively.

Zoom In Zoom Out Reset image size

Table 4.  Maximum Likelihood Fits for 2009–2010 Seasons

Parameter	Opt. Fit.	Prop. Unc.a	Stat. Unc.a	Tot. Unc.a	Cov. Unc.b
	(${p}^{* }$)	(δppr )	(δpst )	(δp )	δpcov
${V}_{\mathrm{ISNO}\infty }$ (km s−1)	25.76	0.29	0.55	0.62	${}_{-4.8}^{+5.3}$
${T}_{\mathrm{ISNO}\infty }$ (K)	8450	910	1750	1970	${}_{-3500}^{+5000}$
${\lambda }_{\mathrm{ISNO}\infty }$ (°)	76.51	2.85c	5.45c	6.15c	NA
${\beta }_{\mathrm{ISNO}\infty }$ (°)	−4.47	0.26	0.50	0.56	NA

Notes.

^aThe quantities $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, $\delta {V}_{\mathrm{ISNO}\infty }$, $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, and $\delta {T}_{\mathrm{ISNO}\infty }$ represent propagation, statistical, and total uncertainties across the parameter tube for interstellar O speed and temperature, respectively. ^bThe quantities $\delta {V}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ and $\delta {T}_{\mathrm{ISNO}\infty }^{\mathrm{cov}}$ represent the uncertainties along the parameter tube covariant with the uncertainty in flow longitude. ^cThe quantities $\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{pr}}$, $\delta {\lambda }_{\mathrm{ISNO}\infty }^{\mathrm{st}}$, and $\delta {\lambda }_{\mathrm{ISNO}\infty }$ represent propagation, statistical, and total uncertainties of interstellar longitude with interstellar speed and temperature varied along the parameter tube.

Download table as:  ASCIITypeset image

Figures 8 and 9 compare model results and observations using the fit results from maximum likelihood. Since there is a finite probability for bins with 0 counts, these observed bin nulls are included in the fit.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on the results of the maximum likelihood analysis, we derive the parameter tube for the interstellar temperature (Figure 10, top) and speed (Figure 10, bottom) as a function of interstellar longitude. The derived temperature parameter tube for interstellar O and interstellar He are almost identical. At the flow longitude of 76°, the O temperature is ∼1000 K greater than the He temperature. This constitutes less than a ∼11% difference.

Zoom In Zoom Out Reset image size

Figure 10 reveals that the derived interstellar speeds are almost statistically distinct at the 1σ level. At the flow longitude of 76°, the interstellar O speed would need to be ∼0.95 km s⁻¹ higher than the interstellar He speed.

It is noteworthy that while the parameter tubes are close to one another in the 2009 and 2010 seasons, the inferred inflow longitudes are quite different. Using only 2009 data, we derive ${\lambda }_{\mathrm{ISNO}\infty }\sim 80\buildrel{\circ}\over{.} 3\pm 8\buildrel{\circ}\over{.} 6$, and using only 2010 data, we derive ${\lambda }_{\mathrm{ISNO}\infty }\sim 71\buildrel{\circ}\over{.} 4\pm 7\buildrel{\circ}\over{.} 8$. The large uncertainties and the large variation from season to season are indications of the significant degeneracy in interstellar parameters along the parameter tube.

The final parameters derived from maximum likelihood fits are the scaling constants A_k for 2009 and 2010. For the 2009 period (k = 1), we find A₁ = (1.1 ± 0.1) × 10⁻⁵ cm⁻3 and for the 2010 period (k = 2), ${A}_{2}=(1.0\pm 0.09)\times {10}^{-5}$ cm⁻³. We relate these scaling constants to an interstellar density in Section 6.

Table 5 shows results for the interstellar O density using a χ² and a maximum likelihood analysis. Since the scaling factor derived from fits represents an interstellar O density (near the termination shock) times the survival probability, we must have an estimate of the O ionization rate, taken from Bzowski et al. (2013a), in order to estimate the density.

Table 5.  Results for the Interstellar O Density

Parameter	2009	2010	2009–10
${\beta }_{1}^{\mathrm{ion}}$ a (10−7 s−1)	5.34 ± 0.53	5.05 ± 0.51	⋯
${n}_{\mathrm{ISNO}\infty }$(χ2) (10−5 cm−3)	${5.2}_{-1.6}^{+2.0}$	${4.6}_{-1.6}^{+1.9}$	${4.9}_{-1.1}^{+1.4}$
${n}_{\mathrm{ISNO}\infty }$(Max. Likelihood) (10−5 cm−3)	${6.3}_{-1.2}^{+1.5}$	${5.2}_{-0.9}^{+1.1}$	${5.8}_{-0.8}^{+0.9}$

Note.

^aIonization rates from Bzowski et al. (2013a).

Download table as:  ASCIITypeset image

The IBEX results can be directly compared to the abundances determined from both pickup ions and anomalous cosmic rays because they are sampling the same population of interstellar atoms (e.g., those where the filtration has already occurred). The interstellar density of He determined from pickup ions is estimated at 1.5 × 10⁻² cm⁻³ (Gloeckler 1996). The O/He abundance from pickup ions is estimated as $[{\rm{O}}/\mathrm{He}]={3.5}_{-1.4}^{+1.8}\times {10}^{-3}$ (Geiss et al. 1994). Taken together, these observations indicate a neutral density in the outer heliosphere (near the termination shock) for O of ${n}_{\mathrm{ISNO}\infty }={5.25}_{-2.1}^{+2.7}\times {10}^{-5}$ cm⁻³. From anomalous cosmic-ray data, Cummings et al. (2002) derives an O density in the outer heliosphere (near the termination shock) of ${n}_{\mathrm{ISNO}\infty }=(5.3\pm 0.8)\times {10}^{-5}$ cm⁻³. Both results are consistent, within uncertainties to the interstellar O density ${n}_{\mathrm{ISNO}\infty }={5.8}_{-0.8}^{+0.9}\times {10}^{-5}$ cm⁻³ found here.

The flow longitude derived for interstellar O of 760 ± 34 from χ² analysis and 765 ± 62 from a maximum likelihood fit is compatible with the derived flow longitude of interstellar He, 75.6 ± 1.4 (Schwadron et al. 2015). Drews et al. (2012) used STEREO data to determine the longitude of the O flow vector using the upwind crescent formed by O pickup ions. They found a flow longitude of 789 ± 31. This longitude is consistent with the IBEX value to within uncertainties. Note also that transport effects acting on pickup ions tend to increase the longitude observed for pickup ions relative to their progenitor neutral atoms. This effect would reduce the inferred flow longitude based on the pickup ion observations (Chalov & Fahr 1999), which likely brings the IBEX and STEREO observations into even closer agreement. Sokół et al. (2016) show that the latitudinal anisotropy of the rate of charge-exchange between ISN O and solar wind protons is due to the latitudinal structure of the solar wind. This anisotropy in latitude likely biases the longitude of the ISN O inflow derived from PUI observations (Drews et al. 2012) by an amount comparable to the difference between the PUI O inflow longitude and the most recent IBEX ISN He inflow longitude (McComas et al. 2015b).

We compare the results derived here for interstellar O (green data point) with interstellar He and other species in Figure 11. This figure is a variant of that developed by Schwadron et al. (2015c) showing the triangulation of the direction of the local interstellar magnetic field. Blue and cyan curves show the great circle projections of connecting interstellar He (blue for IBEX, and cyan for Ulysses; Witte 2004; Schwadron et al. 2015) with the interstellar H measured by SOHO/SWAN (Lallement et al. 2005, 2010). The dashed blue curves show the uncertainty limits of this projection. These curves define what is often referred to as the "neutral deflection plane" of the "B–V plane."

Zoom In Zoom Out Reset image size

The gray data point shows the direction of the ISN He "warm breeze" (Kubiak et al. 2014, 2016), the secondary neutral component created by charge-exchange between ISN He and the plasma in the heliosheath. The warm breeze is both substantially hotter and slower than the primary neutral He from the LISM. Kubiak et al. (2016) found a flow direction of the warm breeze, which is coplanar (the best-fit B–V plane including secondary He, primary He, and H is shown in gray in Figure 11) with the flow directions of ISN H and ISN He, and the direction to the center of IBEX Ribbon within uncertainties (Funsten et al. 2009b, 2013; Fuselier et al. 2009a; McComas et al. 2009a; Schwadron et al. 2009). This co-planarity lends support to the hypothesis that the warm breeze is the secondary population of ISN He and that the center of the Ribbon coincides with the direction of the local interstellar magnetic field.

One of the interesting features of the interstellar O direction found here is that the latitude is below that of both the primary interstellar He component and the secondary component (the warm breeze). Figure 12 provides a schematic that helps explain the latitude shifts in the secondary He warm breeze, and interstellar O, which is more strongly filtered by the heliosheath than primary neutral He. The draping of the local interstellar magnetic field causes compression and shifts the pressure maximum in observed ENAs (in the globally distributed flux, which excludes the ribbon) south of the upwind direction (Schwadron et al. 2014). This added compression south of the upwind direction thins the outer heliosheath, providing preferred access to primary neutrals from the interstellar medium due to weakened filtration. In contrast, secondaries created within the outer heliosheath should have a slight bias for being created at latitudes above the upwind direction. These effects would tend to shift the secondary upwind direction above (north of) the primary He upwind direction in latitude, but shift the interstellar O upwind direction below (south of) the primary He upwind direction in latitude, as observed (Figure 11).

Zoom In Zoom Out Reset image size

We have presented the first detailed analysis of IBEX ISN O observations, providing entirely new information on the flow properties of this important interstellar population. The key areas of development are summarized as follows:

• 1.  
  Deriving the local interstellar properties of oxygen is important for understanding nucleosynthesis. Because interstellar O undergoes significant amounts of charge-exchange, its ionization rate is important for backing out abundances from measurements made inside the heliosphere. One of the key parameters derived from our fitting analysis is the survival probability times the ISN O density, ${n}_{\mathrm{ISNO}\infty }{S}_{p}=(9.2\pm 2.3)\times {10}^{-6}$ cm⁻³ based on χ² fitting and (1.1 ± 0.1) × 10⁻⁵ cm⁻³ based on maximum likelihood fitting. For ionization rates (referenced to 1 au) from Bzowski et al. (2013a; see Table 5), our observations require a density of ${n}_{\mathrm{ISNO}\infty }={4.9}_{-1.1}^{+1.4}\times {10}^{-5}$ cm³ for χ² fitting and ${n}_{\mathrm{ISNO}\infty }={5.8}_{-0.8}^{+0.9}\times {10}^{-5}$ cm⁻³ based on maximum likelihood fitting. These estimates are remarkably close to previous estimates (e.g., Geiss et al. 1994; Gloeckler 1996; Cummings et al. 2002).
• 2.  
  We have derived the interstellar parameter tube relationships for primary interstellar O. The relationship between interstellar temperature and longitude overlaps with and is consistent with the corresponding parameter tube for primary interstellar He. In contrast, the parameter tube for interstellar O speed is systematically higher than for primary interstellar He. This difference is likely the result of enhanced filtration of O atoms through the heliosheath. The interaction of interstellar of O in the heliosheath is somewhat analogous to the interaction of primary interstellar H since both species have relatively large charge-exchange ionization rates. Charge-exchange interactions for both primary neutral H and O in the outer heliosheath represent a loss process, which has been quantified (e.g., Izmodenov et al. 2004). Such loss processes act preferentially on slower atoms due to the energy dependence of the charge-exchange cross section (Izmodenov et al. 2001; Katushkina & Izmodenov 2010). As a result, simulations show that filtration causes a ∼1–2 km s¹ average increase of the primary flow speed and a few hundred kelvin reduction in the effective temperature of the primary distribution function (Izmodenov et al. 2001). Note that we observe a ∼1 km s⁻¹ faster primary population of O compared to He. The signature (faster O flow) is characteristic of the filtration process.
• 3.  
  The upwind flow latitude of interstellar O is +44 ± 03 for χ² fitting and +447 ± 056 for maximum likelihood fitting, which is southward of the interstellar upwind flow latitude of interstellar He, +512 ± 027. The upwind flow latitude of secondary He (the warm breeze) is found by Kubiak et al. (2016) to be +1195 ± 067. The deviation between the warm breeze and primary He flow suggests that the asymmetric latitude structure of the heliosheath leads to the preferential addition of secondary neutral atoms north of the upwind direction, and therefore causes the secondary upwind direction to be northward of the upwind primary direction. This also suggests that the southward deflection of the upwind flow latitude of interstellar O may be the result of filtration in which the compressed heliosheath southward of the primary upwind direction provides preferential access to the heliosphere for primary oxygen atoms (see Figure 12).
• 4.  
  The temperature of primary interstellar O is found to be remarkably consistent with that of He. The optimal fit temperature of oxygen at a reference flow longitude of 76° is ∼600 K less (for χ² fitting) and ∼1000 K greater (for maximum likelihood fitting) than the corresponding helium temperature at this longitude. This difference in temperatures is extremely small, within ∼10%. We measure the states of neutral atoms, which interact with ions in the interstellar medium through collisions and charge-exchange. (Note that neutral O thermodynamics should be closely coupled to neutral H thermodynamics, which in turn is coupled to protons through charge-exchange and through ion-neutral dipole coupling; Spangler et al. 2011.) Collisions will tend to create Maxwell–Boltzmann-shaped velocity distributions with roughly equal temperatures for different species. However, mass-dependent heating, which is often found in even moderately turbulent plasmas such as the solar wind near 1 au (e.g., Hefti et al. 1998), would render ion velocity distributions closer to a state with equal thermal (or random kinetic) speeds. This would cause heavy ion distributions such as O to have much higher temperature than light ion distributions. The fact that ISN atom distributions show roughly equal temperatures for He and O indicates that collisions maintain equilibrium temperatures of ions and neutrals. High levels of turbulence in the interstellar medium would cause differentiation in the temperatures of O and He. Therefore, turbulence must be relatively weak in the local interstellar cloud near the solar system.

We have provided the first detailed analysis of ISN oxygen properties using IBEX data. We have shown that the interstellar oxygen parameters are similar to those of the He primary component. In particular, the best-fit O and He temperatures are within ∼10% at a reference flow longitude of 76°, which provides a strong indication that the local interstellar plasma near the Sun is relatively unaffected by turbulent heating. We also observe key differences between primary oxygen and helium. The parameter tube (relation between speed and longitude) for oxygen has interstellar speeds larger than the corresponding speeds in the parameter tube for He, and the flow latitude for oxygen is closer to the ecliptic plane. Both of these differences may be the result of enhanced filtration of interstellar oxygen due to its higher charge-exchange ionization rate compared to helium. Further, we derive the interstellar O neutral density that, within uncertainties, is consistent with previous estimates. Thus, we use IBEX data for the first time to probe the interstellar properties of oxygen.

We are grateful to the many individuals who have made the IBEX project possible. This work is supported by the Interstellar Boundary Explorer mission as a part of NASA's Explorer Program and partially by NASA SR&T Grant NNG06GD55G, the Swiss National Science Foundation, PRODEX, and the Polish National Science Center grant 2015-18-M-ST9-00036. Support for NAS also provided by a NSF grant (Sun-2-Ice, grant AGS1135432). Finally, we thank the International Space Science Institute for supporting discussions of the maximum likelihood approach used here as part of the research team: Radiation Interactions at Planetary Bodies (http://www.issibern.ch/teams/interactplanetbody/).

The geometric factor of IBEX-Lo for heavy neutral atoms, such as O, is controlled by the combination of surface conversion and sputtering to negative ions. As a consequence, the peak of the effective response in energy E_Peak is always noticeably below the energy of the incoming neutral atom E_In, with this energy reduction factor R = E_Peak/E_In decreasing with the energy E_In of the neutral atoms. In addition, the sensor detects these neutral atoms as negative ions over a broad range in energy around the peak energy E_Peak due to the large acceptance range ΔE/E = 0.7 of the ESA and the broad range over which sputtered ions are generated. Therefore, the geometric factor that needs to be applied to a given combination of incoming energy E_In and IBEX-Lo energy step n (E-Step n = 1 8) depends on the absolute efficiency as a function of E_In, the reduction factor R, and where along the energy response function, relative to E_Peak, the energy of energy step n E_n lies.

To obtain such a geometric factor calibration, the energy response curve of the sensor is needed, how its peak energy varies with the incoming energy of the neutral atoms, and the absolute geometric factor as a function neutral atom energy. It turns out that based on various calibration runs with an O beam during final calibration before flight, a combined normalized energy response curve $\epsilon (E/{E}_{\mathrm{Peak}})$ could be constructed, normalized to the peak efficiency and the peak energy for each beam energy. The best-fit energy response curve to the combined calibration data set yields

$$\begin{eqnarray}&&\epsilon ({E/E}_{\mathrm{Peak}})=\exp \left(-\displaystyle \frac{{(E/{E}_{\mathrm{Peak}}-1)}^{2}}{2{\delta }^{2}}\right),\end{eqnarray} \tag{ 10 }$$

where δ = 0.59. This functional dependence remains constant over the energy range of IBEX-Lo. The reduction factor R(E_In) as a function of beam energy E_In has been derived from a series of energy response measurements as:

$$\begin{eqnarray}&&R({E}_{\mathrm{In}})={R}_{0}-{R}_{1}\mathrm{log}({E}_{\mathrm{In}}),\end{eqnarray} \tag{ 11 }$$

where R₀ = 1.04 and R₁ = 0.134. ISN O arrives during the early spring ISN flow observation season in the IBEX frame at 1 au with a bulk energy of ≈540 eV. If secondary O is included in the analysis, the range of O energies in the IBEX frame falls between 540 and 435 eV, the energy of O atoms arriving on bounding parabolic trajectories starting with 0 eV outside the heliosphere. This energy range is between IBEX-Lo energy E-steps 5 and 6, which are the only steps of interest here. The signals in E-step 4 and the lower steps are completely controlled by ISN He (Möbius et al. 2009, 2012), and the O signal is too weak in E-step 7 to be of much use.

To determine the absolute geometric factor for E-steps 5 and 6, we used calibration results obtained with a neutral O beam E_In = 279 and 601 eV, the nominal incoming energies for these two steps, with the IBEX-Lo response taken at the same energy step and all the way down to energy step 1. Using Equations (10) and (11) the energy response curve was then fitted to the set of calibration results for E-step 5 and 6, adopting the same ±35% uncertainty for the absolute efficiency of each calibration point. To get the values for the two bounding energies of ISN O, the geometric factors at E_Peak for the two E-steps were linearly interpolated. Table 6 shows the resulting geometric factors for the ISN bulk flow energy (540 eV) and for ISN O on a parabolic trajectory (435 eV). Listed are the geometric factors at E_Peak for these input energies at the top, followed by the geometric factors for observations in E-step 5 and 6.

Table 6.  Geometric Factors for ISN O

Parameter	ISN O	ISN O (ISN Traj)
${E}_{1\mathrm{au}}$ (eV)	540	435
GF at EPeak (cm2 sr)	9.6 ± 3.4 10−5	8.7 ± 3.0 10−5
GF at E-Step 5 (cm2 sr)	8.8 ± 3.1 10−5	8.6 ± 3.0 10−5
GF at E-Step 6 (cm2 sr)	5.6 ± 2.0 10−5	2.2 ± 0.8 10−5

Download table as:  ASCIITypeset image

We examine the fitting of the M-parameter statistical model $V(x;{p}_{1},{p}_{2},...,{p}_{M})$ to the given data set $\{{x}_{i},{y}_{i}\pm {\sigma }_{{yi}}\}{}_{i=1}^{N}$, by minimizing the χ²

$$\begin{eqnarray}&&{\chi }^{2}({p}_{1},{p}_{2},...,{p}_{M})=\displaystyle \sum _{i=1}^{N}\,{\sigma }_{{yi}}^{-2}{[{y}_{i}-V({x}_{i};{p}_{1},{p}_{2},\ldots ,{p}_{M})]}^{2}.\end{eqnarray} \tag{ 12 }$$

Here y_i is the number of counts. The global minimum of the χ² give the parameter optimal values, $({p}_{1}^{* },\,{p}_{2}^{* },...,{p}_{M}^{* })$, by solving the system of equations,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\chi }^{2}({p}_{1},{p}_{2},\ldots ,{p}_{M})}{\partial {p}_{1}}=0,\\ &&\displaystyle \frac{\partial {\chi }^{2}({p}_{1},{p}_{2},\ldots ,{p}_{M})}{\partial {p}_{2}}=0,\\ &&\qquad \qquad \qquad ...,\\ &&\displaystyle \frac{\partial {\chi }^{2}({p}_{1},{p}_{2},\ldots ,{p}_{M})}{\partial {p}_{M}}=0.\end{eqnarray} \tag{ 13 }$$

After solving the system of these M equations, we derive the parameter optimal values as functions of data points,

$$\begin{eqnarray}{p}_{1}^{* } & = & {p}_{1}^{* }({y}_{1},{y}_{2},\ldots ,{y}_{N}),\\ {p}_{2}^{* } & = & {p}_{2}^{* }({y}_{1},{y}_{2},\ldots ,{y}_{N}),\\ & ..., & \\ {p}_{M}^{* } & = & {p}_{M}^{* }({y}_{1},{y}_{2},...,{y}_{N}).\end{eqnarray} \tag{ 14 }$$

For simplicity, we have indicated only dependence on data (y) values.

The statistical error (also called the curvature error) of the optimal value ${p}_{m}^{* }$ is given by

$$\begin{eqnarray}&&\delta {p}_{m,\mathrm{st}}^{* }=\sqrt{{\chi }_{\mathrm{red}}^{2}\cdot 2{H}_{{mm}}^{-1}},\end{eqnarray} \tag{ 15 }$$

for m = 1, 2, ..., M and ${\boldsymbol{H}}$ is the Hessian matrix of the χ² at the global minimum. The quantity ${H}_{{mm}}^{-1}$ is the mth diagonal element of the Hessian inverse matrix (Livadiotis 2007). The estimated χ² value is

$$\begin{eqnarray}&&{\chi }_{\mathrm{est}}^{2}={\chi }^{2}({p}_{1}^{* },\,{p}_{2}^{* },...,{p}_{M}^{* }),\end{eqnarray} \tag{ 16 }$$

and ${\chi }_{\mathrm{red}}^{2}={\chi }_{\mathrm{est}}^{2}/(N-M)$ is the reduced χ² value (the degrees of freedom are N − M).

We generalize the statistical error for use in maximum likelihood fitting techniques. Consider an arbitrary distribution function P_i = f(y_i, V_i). The joint probability distribution over the set $\{{y}_{i}\}{}_{i=1}^{N}$ is ${\rm{\Gamma }}={\prod }_{i=1}^{N}{P}_{i}$, or, by taking its logarithm,

$$\begin{eqnarray}&&\mathrm{ln}{\rm{\Gamma }}=\displaystyle \sum _{i}^{N}\,\mathrm{ln}{P}_{i}=\displaystyle \sum _{i=1}^{N}\,\mathrm{ln}f({y}_{y},{V}_{i}).\end{eqnarray} \tag{ 17 }$$

The likelihood is defined by this joint probability, but it is more often used with the inverse function and its logarithm,

$$\begin{eqnarray}&&{\ell }=-c\displaystyle \sum _{i=1}^{N}\,\mathrm{ln}[f({y}_{i},{V}_{i})],\end{eqnarray} \tag{ 18 }$$

where c is an arbitrary constant (typically taken to be c = 2). The inverse likelihood ℓ is a function of the fitting parameters, ℓ(p₁, p₂, ..., p_M) and needs to be minimized, corresponding to the minimization of the χ² in the case of normally distributed data. The Taylor expansion of ${\ell }({p}_{1},{p}_{2},...,{p}_{M})$ around its minimum value, ${{\ell }}^{* }({p}_{1}^{* },\,{p}_{2}^{* },\ldots ,\,{p}_{M}^{* })$, gives

$$\begin{eqnarray}{\ell }({p}_{1},{p}_{2},...,{p}_{M}) & = & {{\ell }}^{* }({p}_{1}^{* },\,{p}_{2}^{* },\ldots ,{p}_{M}^{* })\\ & & +\displaystyle \sum _{i=1}^{N}\,{\left(\displaystyle \frac{\partial {\ell }}{\partial {p}_{i}}\right)}_{{\boldsymbol{p}}={{\boldsymbol{p}}}^{* }}({p}_{i}-{p}_{i}^{* })\\ & & +\displaystyle \frac{1}{2}\displaystyle \sum _{i=1}^{N}\,\displaystyle \sum _{j=1}^{N}\,{\left(\displaystyle \frac{{\partial }^{2}{\ell }}{\partial {p}_{i}\partial {p}_{j}}\right)}_{{\boldsymbol{p}}={{\boldsymbol{p}}}^{* }}\\ & & \times ({p}_{i}-{p}_{i}^{* })({p}_{j}-{p}_{j}^{* }).\end{eqnarray} \tag{ 19 }$$

Since this expansion is performed about an extremum, the second term in the expansion drops out, resulting in the following:

$$\begin{eqnarray}{\ell }({p}_{1},{p}_{2},...,{p}_{M}) & = & {{\ell }}^{* }({p}_{1}^{* },\,{p}_{2}^{* },...,{p}_{M}^{* })\\ & & +\displaystyle \sum _{i=1}^{N}\,\displaystyle \sum _{j=1}^{N}\,{A}_{{ij}}({p}_{i}-{p}_{i}^{* })({p}_{j}-{p}_{j}^{* }),\end{eqnarray} \tag{ 20 }$$

where the curvature matrix ${\boldsymbol{A}}$ is half the Hessian,

$$\begin{eqnarray}&&{A}_{{ij}}=\displaystyle \frac{1}{2}{H}_{{ij}}=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\partial }^{2}{\ell }}{\partial {p}_{i}\partial {p}_{j}}\right)}_{{\boldsymbol{p}}={{\boldsymbol{p}}}^{* }}.\end{eqnarray} \tag{ 21 }$$

The statistical independent uncertainty in this case is given by

$$\begin{eqnarray}&&\delta {p}_{i,\mathrm{st}}\approx \sqrt{\displaystyle \frac{2{{\ell }}^{* }}{N-M}{H}_{{ii}}^{-1}},\end{eqnarray} \tag{ 22 }$$

There is another type of error that characterizes the uncertainty in parameter optimal values and involves the propagation of the measurement uncertainties $\{{\sigma }_{{yi}}\}{}_{i=1}^{N}$. This propagation uncertainty is given by

$$\begin{eqnarray}&&\delta {p}_{m,\mathrm{pr}}^{* }=\sqrt{{\sum }_{i=1}^{N}{\left(\displaystyle \frac{\partial {p}_{m}^{* }}{\partial y}\right)}^{2}{\sigma }_{{yi}}^{2}}\end{eqnarray} \tag{ 23 }$$

for m = 1, 2, ..., M.

The two types of uncertainty are fundamentally distinct: the curvature error is not strongly dependent on the measurement uncertainties $\{{\sigma }_{{yi}}\}{}_{i=1}^{N}$. In contrast, the propagation error is based directly on a weighted sum of the squares of measurement uncertainties. For example, if all measurement uncertainties are equal to σ, the propagation error becomes directly proportional to the common measurement uncertainty, $\delta {p}_{m,\mathrm{pr}}^{* }\,\propto \,\sigma$. Because both uncertainties are distinct, they must both be included in the final uncertainty estimation. A rough approximation for the total uncertainty (see, for example, Schwadron et al. 2013) is given by

$$\begin{eqnarray}&&\delta {p}^{* }=\sqrt{{(\delta {p}_{\mathrm{st}}^{* })}^{2}+{(\delta {p}_{\mathrm{pr}}^{* })}^{2}}.\end{eqnarray} \tag{ 24 }$$

In this section we discuss how minimization of χ² can be used to derive orbit-by-orbit background rates and the scaling factor that multiplies the spin-phase distribution. Observed rate distributions as a function of spin-phase are represented by ${\psi }_{{ij}}^{o}$ where the first index i corresponds to the orbit number, and the second index j refers to the spin-phase. These rate distributions are accumulated over the good times for a specific orbit. The model has dependencies on the instantaneous velocity and pointing of the spacecraft, and the instantaneous pointing of the IBEX-Lo instrument. As a result, the model simulations must also accumulate rates over the good time periods, and over the angles swept out over a given spin-phase bin. Therefore, the rate distribution modeled ${\psi }_{{ij}}^{m}$ simulates the actual instrument response precisely.

The modeled rate distributions scale linearly with a factor that depends on the density times the survival probability, ${n}_{\mathrm{ISNO}\infty }$ × S_p. Therefore, it is convenient to define an overall scaling factor ${A}_{k}=({n}_{\mathrm{ISNO}\infty }/{n}_{\mathrm{ref}})\times ({S}_{p}/{S}_{p-\mathrm{ref}})$ for the modeled rates. Here n_ref and ${S}_{p-\mathrm{ref}}$ are a reference density and a reference survival probability. The scaling factor A_k should be roughly fixed in each interstellar season, but could vary from season to season. The scaling factor therefore depends on the year of observation, indexed by k. The model distribution is expressed:

$$\begin{eqnarray}&&{\psi }_{{ij}(k)}^{m}={A}_{k}{\psi }_{{ij}}^{m^{\prime} }.\end{eqnarray} \tag{ 25 }$$

Only specific orbits i exist in a given year and we treat the index (k) as an implicit value. The quantity ${\psi }_{{ij}}^{m^{\prime} }$ represents simply the unnormalized modeled rate,

$$\begin{eqnarray}&&{\psi }_{{ij}}^{m^{\prime} }=\displaystyle \frac{{\psi }_{{ij}(k)}^{m}}{{n}_{\mathrm{ref}}{S}_{p-\mathrm{ref}}}.\end{eqnarray} \tag{ 26 }$$

This unnormalized model distribution incorporates the effects of survival probability, integration through the instrument response, and all the factors involved in translating fluxes from the outer heliosphere to IBEX at 1 au.

There is also a background rate, B_i, that typically varies from orbit to orbit. We treat this background rate as a constant rate in the spin-phase distribution in a given orbit, which is added to the rate derived from interstellar atoms. Therefore the modeled spin-phase distribution in a given orbit is defined by

$$\begin{eqnarray}&&{M}_{{ij}(k)}={A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i}\end{eqnarray} \tag{ 27 }$$

and the corresponding χ² is given by

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=\mathrm{orbits}}\,\displaystyle \sum _{j=\mathrm{spinphases}}\,{({A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i}-{\psi }_{{ij}}^{o})}^{2}/{\sigma }_{{ij}}^{2}.\end{eqnarray} \tag{ 28 }$$

where σ_ij is the uncertainty in the observed rate.

We solve for the background rate and scaling factor by minimizing the χ² explicitly. We solve $\partial {\chi }^{2}/\partial {B}_{i}=0$, which yields

$$\begin{eqnarray}&&{B}_{i}={\bar{O}}_{i}-{A}_{k}{\bar{M}}_{i}\end{eqnarray} \tag{ 29 }$$

where the effective observed rate ${\bar{O}}_{i}$ is

$$\begin{eqnarray}&&{\bar{O}}_{i}=\displaystyle \frac{{\sum }_{j=\mathrm{spinphases}}{\psi }_{{ij}}^{o}/{\sigma }_{{ij}}^{2}}{{\sum }_{j=\mathrm{spinphases}}1/{\sigma }_{{ij}}^{2}}\end{eqnarray} \tag{ 30 }$$

and the effective model rate ${\bar{M}}_{i}$ is

$$\begin{eqnarray}&&{\bar{M}}_{i}=\displaystyle \frac{{\sum }_{j=\mathrm{spinphases}}{\psi }_{{ij}}^{m^{\prime} }/{\sigma }_{{ij}}^{2}}{{\sum }_{j=\mathrm{spinphases}}1/{\sigma }_{{ij}}^{2}}.\end{eqnarray} \tag{ 31 }$$

In order to solve for the scaling factor A_k and its uncertainty $\delta {A}_{k}^{* }$, we express the χ² as a quadratic function of A_k with a minimum at ${A}_{k}^{* }$:

$$\begin{eqnarray}&&{\chi }^{2}={C}_{0}+{C}_{2}{({A}_{k}-{A}_{k}^{* })}^{2}.\end{eqnarray} \tag{ 32 }$$

By setting Equation (32) equal to Equation (28) we may solve for the coefficients C₀, ${A}_{k}^{* }$, and C₂:

$$\begin{eqnarray}&&{C}_{2}=\displaystyle \sum _{i=\mathrm{orbits}}\,\displaystyle \sum _{j=\mathrm{spinphases}}\,{({\psi }_{{ij}}^{m^{\prime} })}^{2}/{\sigma }_{{ij}}^{2}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{A}_{k}^{* }={C}_{2}^{-1}\displaystyle \sum _{i=\mathrm{orbits}}\,\displaystyle \sum _{j=\mathrm{spinphases}}\,({\psi }_{{ij}}^{o}-{B}_{i}){\psi }_{{ij}}^{m^{\prime} }/{\sigma }_{{ij}}^{2}\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{C}_{0}=\displaystyle \sum _{i=\mathrm{orbits}}\,\displaystyle \sum _{j=\mathrm{spinphases}}\,{({\psi }_{{ij}}^{o}-{B}_{i})}^{2}/{\sigma }_{{ij}}^{2}-{C}_{2}{({A}_{k}^{* })}^{2}.\end{eqnarray} \tag{ 35 }$$

This formulation is similar to that used by Schwadron et al. (2013), where it is shown that the statistical uncertainty of the χ² minimum from a quadratic form is

$$\begin{eqnarray}&&\delta {A}_{k,\mathrm{st}}^{* }=\sqrt{{C}_{0}/(\nu {C}_{2})},\end{eqnarray} \tag{ 36 }$$

where ν is the number of degrees of freedom.

The propagation uncertainty $\delta {A}_{k,\mathrm{pr}}^{* }$ follows from the traditional treatment of error propagation

$$\begin{eqnarray}\delta {A}_{k,\mathrm{pr}}^{* } & = & \sqrt{{\displaystyle \sum }_{i=\mathrm{orbits}}\,{\displaystyle \sum }_{j=\mathrm{spinphases}}{\left(\displaystyle \frac{\partial {A}_{k}^{* }}{\partial {\psi }_{{ij}}^{o}}\right)}^{2}{\sigma }_{{ij}}^{2}}\\ & = & 1/\sqrt{{C}_{2}}.\end{eqnarray} \tag{ 37 }$$

We note further that C₀ is the minimum of the χ² and that the reduced χ² (denoted ${\chi }_{\mathrm{red}}^{2}$) at the global minimum is given by,

$$\begin{eqnarray}&&{\chi }_{\mathrm{red}}^{2}={C}_{0}/\nu .\end{eqnarray} \tag{ 38 }$$

Therefore, the statistical uncertainty is directly related to the propagation uncertainty,

$$\begin{eqnarray}&&{(\delta {A}_{k,\mathrm{st}}^{* })}^{2}={(\delta {A}_{k,\mathrm{pr}}^{* })}^{2}{\chi }_{\mathrm{red}}^{2},\end{eqnarray} \tag{ 39 }$$

which results in a total uncertainty given by

$$\begin{eqnarray}&&{(\delta {A}_{k}^{* })}^{2}=(1+{\chi }_{\mathrm{red}}^{2}){(\delta {A}_{k,\mathrm{pr}}^{* })}^{2}.\end{eqnarray} \tag{ 40 }$$

This result can be readily interpreted. The minimum reduced χ² amplifies the propagation uncertainty. A poor fit results in a large value of ${\chi }_{\mathrm{red}}^{2}$ and thereby increases the total uncertainty. The propagation uncertainty represents the minimum possible uncertainty in the presence of a particularly good fit. A statistically likely fit with ${\chi }_{\mathrm{red}}^{2}\sim 1$ results in a total uncertainty that is $\sqrt{2}$ larger than the propagation uncertainty.

In this case, there are two sets of variables to determine from explicit χ² minimization: the background rate B_i determined on an orbit-by-orbit basis; and the scaling factor A_k determined for each interstellar season in a given year. The task is to find a χ² minimum as a function of both the background rates and the scaling factor. We find this two-dimensional (2D) χ² minimum through successive iterations: we find the background rates, then use those background rates in the solution for the amplitude, and then solve again for the background rates, and so forth. This scheme converges quite quickly (typically within three steps) to the 2D χ² minimum.

We use maximum likelihood to calculate background rates and a model scaling factor in a manner similar to the derivation in Appendix C. Specifically, the model is specified as

$$\begin{eqnarray}&&{M}_{{ij}(k)}={A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i},\end{eqnarray} \tag{ 41 }$$

where A_k is the scaling factor and B_i is the background rate. In this case, we must specify the expected counts ${m}_{{ij}(k)}$,

$$\begin{eqnarray}&&{m}_{{ij}(k)}=({A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i}){\rm{\Delta }}{t}_{{ij}(k)}\end{eqnarray} \tag{ 42 }$$

where ${\rm{\Delta }}{t}_{{ij}(k)}$ is the exposure time. Using this expectation value, the probability of observing counts ${o}_{{ij}(k)}$ is

$$\begin{eqnarray}&&{P}_{{ij}(k)}=\displaystyle \frac{{m}_{{ij}(k)}^{{o}_{{ij}(k)}}\exp (-{m}_{{ij}(k)})}{{o}_{{ij}(k)}!}.\end{eqnarray} \tag{ 43 }$$

The inverse likelihood is expressed as

$$\begin{eqnarray}{{\ell }}_{k} & = & -2\displaystyle \sum _{{ij}(k)}\,\mathrm{ln}({P}_{{ij}(k)})\\ & = & -2\displaystyle \sum _{{ij}(k)}\,[{o}_{{ij}(k)}\mathrm{ln}{m}_{{ij}(k)}-{m}_{{ij}(k)}-\mathrm{ln}({o}_{{ij}(k)}!)],\end{eqnarray} \tag{ 44 }$$

where the summation extends over only those orbits associated with a given season k.

The minimization of ℓ_k with respect to the background rate B_i, $\partial {{\ell }}_{k}/\partial {B}_{i}=0$, leads to the relation

$$\begin{eqnarray}&&\displaystyle \sum _{j}\,\displaystyle \frac{{o}_{{ij}(k)}}{{A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i}}=\displaystyle \sum _{j}\,{\rm{\Delta }}{t}_{{ij}(k)}.\end{eqnarray} \tag{ 45 }$$

The constant B_i is varied using a search algorithm to find the optimal value. We minimize the inverse likelihood with respect to the amplitude, A_k, $\partial {{\ell }}_{k}/\partial {A}_{k}=0$, leading to

$$\begin{eqnarray}&&\displaystyle \sum _{{ij}}\,\displaystyle \frac{{o}_{{ij}(k)}{\psi }_{{ij}}^{m^{\prime} }{\rm{\Delta }}{t}_{{ij}(k)}}{{A}_{k}{\psi }_{{ij}}^{m^{\prime} }+{B}_{i}}=\displaystyle \sum _{{ij}}\,{\psi }_{{ij}}^{m^{\prime} }{\rm{\Delta }}{t}_{{ij}(k)}.\end{eqnarray} \tag{ 46 }$$

where the summation extends over only those orbits associated with a given season k. A search algorithm is used to find the optimal amplitude.

As discussed in Appendix B, there are two forms of uncertainty included in our calculations: the statistical and propagation uncertainty. In a previous section (Appendix B), because we are able to solve for the model scaling factor analytically, it is relatively straightforward to directly solve for both the statistical and propagation uncertainties. The difficulty arises however in nonlinear forward models where there is no closed form analytic solution to determine best χ² fits for model parameters. In this case, we must find alternative solutions to derive the best-fit parameters and their uncertainties. This appendix develops a relatively straightforward and robust technique for deriving these uncertainties.

The χ² analysis used in this paper is five dimensional. It is useful here to consider a nonlinear forward model with one parameter, p. The ${\chi }^{2}={\chi }^{2}(p)$ is minimized to find a best-fit at p = ${p}^{* }$. We consider a series of N measurements $\{{O}_{i}\}{}_{i=1}^{N}$ with uncertainties $\{{\sigma }_{i}\}{}_{i=1}^{N}$ and corresponding simulated values $\{{M}_{i}\}{}_{i=1}^{N}$. The χ² is defined:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{N}\,{\left(\displaystyle \frac{{M}_{i}-{O}_{i}}{{\sigma }_{i}}\right)}^{2}.\end{eqnarray} \tag{ 47 }$$

Because the χ² has a minimum, its dependence on p can be approximated as a quadratic:

$$\begin{eqnarray}&&{\chi }^{2}={C}_{0}+{C}_{2}{(p-{p}^{* })}^{2}\end{eqnarray} \tag{ 48 }$$

The statistical uncertainty is then

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{st}}^{* })}^{2}=\displaystyle \frac{{C}_{0}}{\nu {C}_{2}},\end{eqnarray} \tag{ 49 }$$

where ν is the number of degrees of freedom: $\nu =N-M$ and the number of model parameters is M = 1.

The propagation uncertainty is calculated as follows:

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{pr}}^{* })}^{2}=\displaystyle \sum _{i=1}^{N}\,{\left(\displaystyle \frac{\partial {p}^{* }}{\partial {O}_{i}}\right)}^{2}{\sigma }_{i}^{2}.\end{eqnarray} \tag{ 50 }$$

The difficulty in this case is that there is no closed-form solution for ${p}^{* }({O}_{i})$, which makes it cumbersome to solve Equation (50).

One option for solution is to approach the problem numerically. We go through each measurement and vary the measured value through O_i ± σ_i. As we perform this variation, we solve for the χ² and find how variation in O_i causes variation in the χ² minimum. This allows us to solve for $\partial {p}^{* }/\partial {O}_{i}$. We then move to the next data point $i+1$, vary the data point ${O}_{i+1}\pm {\sigma }_{i+1}$, find the change in ${p}^{* }$, and solve for $\partial {p}^{* }/\partial {O}_{i+1}$.

There are two problems with this solution. First, if the number of measurements N is large, which is typical, then the calculation becomes extremely costly to compute. The second problem is that finding the χ² minimum often involves interpolation. Therefore, finding the effects of extremely small changes to χ² and the resulting small changes in ${p}^{* }$ is likely to be difficult to compute accurately.

There is another way to approximate the changes to ${p}^{* }$, that is far more straightforward to implement and interpret. We consider the solution for ${p}^{* }$ as a weighted average of individual terms p_i, which represent the model parameter based on a given observation O_i. The average for ${p}^{* }$ is weighted by the inverse of parameter variances ${\sigma }_{{pi}}^{-2}$ (to be detailed as follows):

$$\begin{eqnarray}&&{p}^{* }={\left(\displaystyle \sum _{i=1}^{N}{\sigma }_{{pi}}^{-2}\right)}^{-1}\displaystyle \sum _{i=1}^{N}\,{\sigma }_{{pi}}^{-2}{p}_{i}.\end{eqnarray} \tag{ 51 }$$

In formulating this solution, we use the model to estimate how individual terms in the summation, p_i, change based on variations in the observed data O_i. Specifically, we use partial derivative $\partial {M}_{i}/\partial p$ to estimate how a change in a given observation ΔO_i leads to a change in the model parameter Δp_i:

$$\begin{eqnarray}&&{\rm{\Delta }}{p}_{i}={\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{-1}{\rm{\Delta }}{O}_{i}.\end{eqnarray} \tag{ 52 }$$

In order to turn this estimate of Δp_i into an estimate of p_i itself we must simply include a constant ${\bar{O}}_{i}$, such that

$$\begin{eqnarray}&&{p}_{i}\approx {\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{-1}({O}_{i}+{\bar{O}}_{i}).\end{eqnarray} \tag{ 53 }$$

Since we are concerned only with small changes in p_i resulting from changes in O_i, it is important only that the constant ${\bar{O}}_{i}$ is fixed while O_i is varied.

Similar reasoning is applied to find the variance of the parameter ${\sigma }_{{pi}}^{2}$ based on the variance observed ${\sigma }_{i}^{2}$ using $\partial {M}_{i}/\partial p$ to convert observed variances into parameter variances:

$$\begin{eqnarray}&&{\sigma }_{{pi}}^{2}={\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{-2}{\sigma }_{i}^{2}.\end{eqnarray} \tag{ 54 }$$

By substituting Equations (53) and (54) into (51), we recover the following solution for ${p}^{* }$

$$\begin{eqnarray}{p}^{* } & \approx & {\left(\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{-2}{\left[\displaystyle \frac{\partial {M}_{i}}{\partial p}\right]}^{2}\right)}^{-1}\\ & & \times \displaystyle \sum _{i=1}^{N}\,{\sigma }_{i}^{-2}\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)({O}_{i}+{\bar{O}}_{i}).\end{eqnarray} \tag{ 55 }$$

Equation (55) is now substituted into (50) to find the variance in ${p}^{* }$:

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{pr}}^{* })}^{2}={\left(\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{-2}{\left[\displaystyle \frac{\partial {M}_{i}}{\partial p}\right]}^{2}\right)}^{-1}.\end{eqnarray} \tag{ 56 }$$

This solution for the propagation uncertainty may be used directly. There is a further simplification when we consider the original form of the χ² minimum. We treat the model solutions as a second-order expansion about the optimal solution:

$$\begin{eqnarray}&&{M}_{i}\approx {M}_{i}^{* }+{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* }(p-{p}^{* })+\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{\partial }^{2}{M}_{i}}{\partial {p}^{2}}\right)}^{* }{(p-{p}^{* })}^{2}.\end{eqnarray} \tag{ 57 }$$

Substitution of this second-order expansion into the χ², Equation (47), results in the following when we retain terms up to the second order in $(p-{p}^{* })$:

$$\begin{eqnarray}{\chi }^{2} & = & \displaystyle \sum _{i=1}^{N}\displaystyle \frac{{({M}_{i}^{* }-{O}_{i})}^{2}}{{\sigma }_{i}^{2}}+2(p-{p}^{* })\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* }\displaystyle \frac{({M}_{i}^{* }-{O}_{i})}{{\sigma }_{i}^{2}}\\ & & +{(p-{p}^{* })}^{2}\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{-2}\left[{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* 2}\right.\\ & & +\left.{\left(\displaystyle \frac{{\partial }^{2}{M}_{i}}{\partial {p}^{2}}\right)}^{* }({M}_{i}^{* }-{O}_{i})\right].\end{eqnarray} \tag{ 58 }$$

Because the solution exists at the minimum of the χ², the term that is the first order in $(p-{p}^{* })$ must vanish:

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{N}\,{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* }\displaystyle \frac{({M}_{i}^{* }-{O}_{i})}{{\sigma }_{i}^{2}}=0.\end{eqnarray} \tag{ 59 }$$

Notice also that we may now relate the constants in Equation (48) to the terms in the expansion:

$$\begin{eqnarray}&&{C}_{0}=\displaystyle \sum _{i=1}^{N}\,\displaystyle \frac{{({M}_{i}^{* }-{O}_{i})}^{2}}{{\sigma }_{i}^{2}},\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{C}_{2}\approx \displaystyle \sum _{i=1}^{N}\,{\sigma }_{i}^{-2}{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* 2},\end{eqnarray} \tag{ 61 }$$

where we have neglected the second-order derivative ${({\partial }^{2}{M}_{i}/\partial {p}^{2})}^{* }$ from C₂. This approximation, which is commonly used, simplifies the analysis and, in practice, does not substantially change the amplitude of the second order term. This association, when applied to Equation (56) implies that

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{pr}}^{* })}^{2}\approx {C}_{2}^{-1}.\end{eqnarray} \tag{ 62 }$$

As in Appendix B, this then implies that

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{st}}^{* })}^{2}={\chi }_{\mathrm{red}}^{2}{(\delta {p}_{\mathrm{pr}}^{* })}^{2}\end{eqnarray} \tag{ 63 }$$

where ${\chi }_{\mathrm{red}}^{2}={\chi }_{\mathrm{est}}^{2}/(N-M)$ is the minimum value of the reduced χ². The total uncertainty is

$$\begin{eqnarray}&&{(\delta {p}^{* })}^{2}=(1+{\chi }_{\mathrm{red}}^{2}){(\delta {p}_{\mathrm{pr}}^{* })}^{2}.\end{eqnarray} \tag{ 64 }$$

As was found previously in Appendix B, the statistical uncertainty is an amplification of the propagation uncertainty where the term ${\chi }_{\mathrm{red}}^{2}$ serves as an amplification factor.

The propagation uncertainty is related to the covariance. For example, in Press et al. (1992), the curvature matrix (one half times the Hessian matrix) is found to be

$$\begin{eqnarray}&&{\alpha }_{{kl}}=\displaystyle \sum _{i=1}^{N}\,\displaystyle \frac{1}{{\sigma }_{i}^{2}}{\left(\displaystyle \frac{\partial {M}_{i}}{\partial {p}_{k}}\displaystyle \frac{\partial {M}_{i}}{\partial {p}_{l}}\right)}^{* },\end{eqnarray} \tag{ 65 }$$

where p_k is a given parameter in a multidimensional fit. The covariance matrix ${\boldsymbol{C}}$ is the inverse of the matrix [α]:

$$\begin{eqnarray}&&{\boldsymbol{C}}={[\alpha ]}^{-1}\end{eqnarray} \tag{ 66 }$$

and the diagonal elements of the covariance matrix constitute the variances of the parameters involved in the fit, which we show are in fact the propagation uncertainties,

$$\begin{eqnarray}&&{(\delta {p}_{k,\mathrm{pr}}^{* })}^{2}={C}_{{kk}},\end{eqnarray} \tag{ 67 }$$

When applied to a single-dimensional fit, we find that the curvature matrix reduces to the following constant

$$\begin{eqnarray}&&\alpha =\displaystyle \sum _{i=1}^{N}\,\displaystyle \frac{1}{{\sigma }_{i}^{2}}{\left(\displaystyle \frac{\partial {M}_{i}}{\partial p}\right)}^{* 2},\end{eqnarray} \tag{ 68 }$$

which is identically equal to C₂. And the inverse of α is identically equal to the variance of the propagation uncertainty, ${(\delta {p}_{\mathrm{pr}})}^{2}=1/\alpha$, as found previously, Equation (62). Substitution of Equation (62) into (63) therefore yields

$$\begin{eqnarray}{(\delta {p}_{\mathrm{st}}^{* })}^{2} & = & {\chi }_{\mathrm{red}}^{2}/\alpha \\ & = & {\chi }_{\mathrm{red}}^{2}\cdot 2{H}^{-1}\end{eqnarray} \tag{ 69 }$$

where H is the Hessian. Equation (69) is nothing less than a 1D form of Equation (15). The derivation comes full-circle, revealing the consistency of the derived terms for the statistical and propagation uncertainties.

The discussion in Section 3 applies to the derivation of the parameter tube for interstellar parameters, which exists in 4D parameter space. This constitutes a multi-parameter fit using χ² minimization. The parameter tube is itself an expression of the underlying relationship between interstellar parameters. And because the interstellar speed and interstellar temperature are varied along the parameter tube when finding the χ² minimum in the interstellar longitude, the covariance between parameters is handled explicitly. As a result, it is appropriate to treat uncertainties using the 1D equations for the propagation and statistical uncertainties: Equations (62) and (63). It must be appreciated however that the uncertainty in interstellar speed and interstellar temperature define the width of the parameter tube in their respective dimensions. The uncertainty in the interstellar longitude applies along the parameter tube and provides the most significant covariant uncertainty in determining the interstellar speed and temperature.

It is also apparent that utilization of the covariance alone determines only the propagation uncertainty. The statistical uncertainty is an additional and independent contribution that inflates the total uncertainty.

The main difference between χ² minimization and maximum likelihood is the approach to metric used to relate observations and model predictions. We consider bin i of the number of counts observed o_i or modeled m_i. Note that there is an important distinction between the counts within bins used here and the rates used in Appendices B and C. The exposure time Δt_i is used to relate rates to observed or modeled counts:

$$\begin{eqnarray}&&{o}_{i}={O}_{i}{\rm{\Delta }}{t}_{i}\end{eqnarray} \tag{ 70 }$$

$$\begin{eqnarray}&&{m}_{i}={M}_{i}{\rm{\Delta }}{t}_{i}.\end{eqnarray} \tag{ 71 }$$

We treat the observed events within each bin o_i as a random Poisson distribution with expectation value m_i,

$$\begin{eqnarray}&&{P}_{i}({m}_{i})=\displaystyle \frac{{m}_{i}^{{o}_{i}}\exp (-{m}_{i})}{{o}_{i}!}.\end{eqnarray} \tag{ 72 }$$

The χ² is related to the joint probability associated with a Gaussian distribution, P_Gi, by ${\chi }^{2}=-2{\sum }_{i}\mathrm{ln}({P}_{{Gi}})$. The corresponding quantity in the case of Poisson distribution is

$$\begin{eqnarray}&&{\ell }=-2\displaystyle \sum _{i}\,\mathrm{ln}({P}_{i})=-2\displaystyle \sum _{i}\,({o}_{i}\mathrm{ln}({m}_{i})-{m}_{i}-\mathrm{ln}[{o}_{i}!]).\end{eqnarray} \tag{ 73 }$$

The maximum likelihood is found by minimizing ℓ with respect to a parameter p:

$$\begin{eqnarray}\displaystyle \frac{\partial {\ell }}{\partial p} & = & 0\\ & = & 2\displaystyle \sum _{i}\,\left(\displaystyle \frac{{m}_{i}-{o}_{i}}{{m}_{i}}\right)\displaystyle \frac{\partial {m}_{i}}{\partial p}.\end{eqnarray} \tag{ 74 }$$

The counterpart for solving for an extremum in the minimization of χ² is

$$\begin{eqnarray}\displaystyle \frac{\partial {\chi }^{2}}{\partial p} & = & 0\\ & = & 2\displaystyle \sum _{i}\,\left(\displaystyle \frac{{m}_{i}-{o}_{i}}{{\sigma }_{i}^{2}{\rm{\Delta }}{t}_{i}^{2}}\right)\displaystyle \frac{\partial {m}_{i}}{\partial p}.\end{eqnarray} \tag{ 75 }$$

The variance can be derived precisely based on the curvature of $\mathrm{ln}({P}_{i})$:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\sigma }_{i}^{2}}=\displaystyle \frac{{\partial }^{2}(-\mathrm{ln}{P}_{i})}{\partial {M}_{i}^{2}}.\end{eqnarray} \tag{ 76 }$$

It should be noted in this relation that the inverse variance is related to the curvature with respect to variation in the model, as opposed to the curvature with respect to variation in observation. The reason for basing the uncertainty on model variation is that the model is a continuous variable, and, as a result, the curvature can be defined accurately even when the observed counts are low. Applying the relation in Equation (76) to the Poisson distribution, we find the following definition for the inverse variance:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\sigma }_{i}^{2}}=\displaystyle \frac{{o}_{i}}{{M}_{i}^{2}}.\end{eqnarray} \tag{ 77 }$$

We apply the same approach as that in Appendix C to derive the propagation uncertainty,

$$\begin{eqnarray}{(\delta {p}_{\mathrm{pr}}^{* })}^{2} & = & {\left(\displaystyle \sum _{i=1}^{N}{\sigma }_{i}^{-2}{\left[\displaystyle \frac{\partial {M}_{i}}{\partial p}\right]}^{2}\right)}^{-1}\\ & = & {\left(\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{o}_{i}}{{m}_{i}^{2}}{\left[\displaystyle \frac{\partial {m}_{i}}{\partial p}\right]}^{2}\right)}^{-1},\end{eqnarray} \tag{ 78 }$$

where we have used Equation (77) to relate the variance to expected counts in a bin. We expand the ℓ to the second order in $(p-{p}^{* })$:

$$\begin{eqnarray}{\ell } & = & {{\ell }}^{* }+2(p-{p}^{* })\displaystyle \sum _{i}{\left(\displaystyle \frac{\partial {m}_{i}}{\partial p}\right)}^{* }\\ & & +{(p-{p}^{* })}^{2}\displaystyle \sum _{i}{\left(\displaystyle \frac{{\partial }^{2}{m}_{i}}{\partial {p}^{2}}\right)}^{* }{\left(\displaystyle \frac{{m}_{i}-{o}_{i}}{{m}_{i}}\right)}^{* }\\ & & +{(p-{p}^{* })}^{2}\displaystyle \sum _{i}{\left(\displaystyle \frac{\partial {m}_{i}}{\partial p}\right)}^{* 2}{\left(\displaystyle \frac{{o}_{i}}{{m}_{i}^{2}}\right)}^{* },\end{eqnarray} \tag{ 79 }$$

where ${{\ell }}^{* }={\ell }({p}^{* })$. We neglect terms involving second derivatives of m_i and the first-order term vanishes since we are at an extremum. With these simplifications, we find the following quadratic form for ℓ:

$$\begin{eqnarray}&&{\ell }={C}_{0}+{C}_{2}{(p-{p}^{* })}^{2}\end{eqnarray} \tag{ 80 }$$

$$\begin{eqnarray}&&{C}_{0}={{\ell }}^{* }\end{eqnarray} \tag{ 81 }$$

$$\begin{eqnarray}{C}_{2} & = & \displaystyle \sum _{i}{\left(\displaystyle \frac{\partial {m}_{i}}{\partial p}\right)}^{* 2}{\left(\displaystyle \frac{{o}_{i}}{{m}_{i}^{2}}\right)}^{* }\\ & = & \displaystyle \sum _{i}{\sigma }_{i}^{-2}{\left(\displaystyle \frac{\partial {m}_{i}}{\partial p}\right)}^{* 2}.\end{eqnarray} \tag{ 82 }$$

Therefore, we find that the propagation and statistical uncertainties are

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{pr}}^{* })}^{2}={C}_{2}^{-1}\end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray}&&{(\delta {p}_{\mathrm{st}}^{* })}^{2}=\displaystyle \frac{{{\ell }}^{* }}{(N-M){C}_{2}}.\end{eqnarray} \tag{ 84 }$$