Distance to the Energetic Neutral Hydrogen Source from the Heliotail

In this study we solve for the distance to the ENA source, specifically, that in the fast SW propagating through the heliotail. Our method relies on IBEX observations of high energy ENAs from both the slow and fast SW in the IHS, as we describe below. The distance to the ENA source also depends on the latitudinal structure of the SW. Figure 3 illustrates the simplified SW structure in a meridional cut through the heliosphere. At solar minimum the fast SW is emitted from PCHs centered near the poles and extending down to midlatitudes (McComas et al. 1998; Figure 3(a)). As the SW flows past the TS, it is carried down the heliotail (Figure 3(b)).

Zoom In Zoom Out Reset image size

In Figure 3(b) we show a simplified geometry of the fast and slow SW as it propagates down the heliotail, assuming that the SW flow trajectory is directed downwind immediately after crossing the TS. We discuss the implications of this assumption in Section 4. The system of equations for this geometry is

$$\begin{eqnarray}\begin{array}{rcl}\sin \alpha & = & \displaystyle \frac{{L}_{1}}{{R}_{1}},\\ \sin \beta & = & \displaystyle \frac{{L}_{1}}{B+{R}_{2}}=\displaystyle \frac{{L}_{1}+{L}_{2}}{A+B+{R}_{2}}.\end{array}\end{eqnarray} \tag{ 1 }$$

Angle α defines the elevation angle from the equatorial plane of the fast–slow SW boundary emitted from the Sun during solar minimum. Angle β defines the elevation angle at which IBEX observes the transition from low to high energy ENAs emitted from the heliotail. We aim to solve for the sum of distances A + B, where A is the LOS path length in the fast SW (up to the cooling length surface), and B is the LOS path length in the slow SW.

First, we determine angle α from Ulysses SWOOPS observations of the SW. IBEX began taking measurements in 2009, near the end of solar cycle 23. Considering that it takes at least 2 yr for the delay between measurements of the SW and ENAs from the heliosphere (e.g., McComas et al. 2017, 2018; Zirnstein et al. 2017a; Reisenfeld et al. 2019), the timing of Ulysses' third, fast polar scan when it was ∼1.4–2.4 au from the Sun approximately corresponds to IBEX observations from 2009 and later. According to observations of the north and south PCH areas (Karna et al. 2014), the PCHs extended to their minimum latitudes around the time of Ulysses observations in ∼2007. Therefore, Ulysses observations during its third polar scan are a reasonable measurement of the minimum latitude that the large-scale, PCH fast SW was emitted from the Sun.

Figure 4 shows SW speed as a function of latitude for Ulysses' third polar scan in 2007. We determine angle α by finding the highest latitudes at which Ulysses crosses into the slow SW of the streamer belt—these crossings are seen as abrupt changes in SW speed. Starting from the poles, we compute a running standard deviation in SW speed using the three nearest bins in latitude and find when the deviation first increases above 50 km s⁻¹. This latitude is found in both the northern and southern hemispheres, yielding the angular width of the slow SW (vertical dashed blue lines). By halving the angular width of the slow SW from Ulysses observations, we find the elevation, or half, angle α = 368 ± 044. The uncertainty of α is calculated as the square root sum of the squares of the uncertainties of each latitude found by the running standard deviation. Note that α is the minimum latitude that the fast SW is emitted from the PCHs in the IBEX epoch, and thus does not change over the times we observe.

Zoom In Zoom Out Reset image size

We note that there are potential systematic uncertainties in α. Specifically, the latitude at which Ulysses crosses into the slow SW streamer belt changes as a function of solar longitude as the Sun rotates with respect to Ulysses position. This is visible in Figure 4. However, we choose angle α at the highest latitude that Ulysses crossed the boundary, since we expect fast–slow SW stream interactions within the streamer belt to wear down the SW structure at large distances from the Sun and thus decrease the SW speed at low- to midlatitudes. The abrupt changes in SW speed observed by Ulysses at low- to midlatitudes is likely from corotating interaction regions (CIRs; e.g., Borovikov et al. 2012). Borovikov et al. (2012) demonstrated that as CIRs wear down with distance from the Sun, the low latitude slow SW may extend to even higher latitudes. We note that this behavior, if true, would imply our value for α would be underestimated. However, as explained later in Section 3.4, we provide a minimum distance to the ENA source in the heliotail, and an underestimation of α also implies an underestimation of the ENA source distance, consistent with our results.

To check the validity of our value for α, we also analyze SW speeds reconstructed from interplanetary scintillation observations of the SW (Tokumaru et al. 2010, 2012) following the methodology of Sokół et al. (2013, 2020). We calculate the average SW latitudinal structure from 2006.5 to 2009.5, which overlaps Ulysses observations and the minimum of solar cycle 23. We find that the latitude at which the interplanetary scintillations-derived SW speed equals a value typically used to identify the fast SW boundary (650 km s⁻¹) is also 368. This reinforces the results derived from Ulysses observations.

Next, we define β as the elevation angle where the ENA spectral index is balanced by ENAs from the slow and fast SW in the IHS (ENAs produced along distances B and A, respectively). To first order this implies that the LOS-integrated fluxes produced along A are equal to those produced along B, i.e., A = B (note, however, there is a small correction to this, discussed below). One might assume this means the ENA spectral index at angle β is exactly halfway between that from the slow versus fast SW. While this is a reasonable approximation, in general, this is not true if the ENA fluxes from the slow and fast SW are scaled vastly differently. To solve this correctly, we first calculate the ENA flux along β at energy E as the weighted mean of ENA fluxes from the slow and fast SW (e.g., Appendix C in Zirnstein et al. 2016b),

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\beta }\left(E\right) & = & \sum _{i=1}^{N}\left(\displaystyle \frac{{J}_{{\rm{s}},i}\left(E\right)}{{\sigma }_{J{\rm{s}},i}^{2}(E)}+\displaystyle \frac{{J}_{{\rm{f}},i}\left(E\right)}{{\sigma }_{J{\rm{f}},i}^{2}(E)}\right)/\sum _{i=1}^{N}\left(\displaystyle \frac{1}{{\sigma }_{J{\rm{s}},i}^{2}(E)}+\displaystyle \frac{1}{{\sigma }_{J{\rm{f}},i}^{2}(E)}\right),\\ {\sigma }_{J\beta }\left(E\right) & = & 1/\sqrt{\sum _{i=1}^{N}\left(\displaystyle \frac{1}{{\sigma }_{J{\rm{s}},i}^{2}(E)}+\displaystyle \frac{1}{{\sigma }_{J{\rm{f}},i}^{2}(E)}\right)}.\end{array}\,\end{eqnarray} \tag{ 2 }$$

We calculate J_β and its uncertainty σ_Jβ at 1.74, 2.73, and 4.29 keV using IBEX data in directions where the minimum (fast SW, J_f) and maximum (slow SW, J_s) spectral indices are found in 2009–2013 (representing a time period reflecting solar minimum conditions), and then compute the spectral index γ_β over J_β from 1.74 to 4.29 keV. This represents the spectral index value at which we need to compute the width of the ENA spectral index distribution shown in Figure 5(a). Note that Equation (2) is averaged over multiple directions, i = 1, ..., N. The directions are chosen as the N lowest and N highest spectral index values in Figure 5(a). For the results presented in this study, we choose N  = 3, i.e., the 3 lowest and 3 highest points from Figure 5(a), yielding γ_β = 2.07 ± 0.04. However, we tested our results for N  = 1, 2, 3, 4, and 5. Our choice of N = 3 does not significantly affect the calculation of angle β (results vary by <8% between N  = 1, ..., 5).

Zoom In Zoom Out Reset image size

To calculate the width of the spectral index distribution in Figure 5(a), we find the intersections of γ_β with the time-averaged data on both sides of the spectral peak and determine the latitudes of these intersections (vertical dashed green lines; Figure 5(a)). Halving the angular width of this region yields β. The uncertainty in β is calculated by fitting two separate lines to the two IBEX data points nearest the intersection, including the extremes of their uncertainties, i.e., fit lines to [(θ_i, γ_i + σ_γi), (θ_i+1, γ_i+1 − σ_γi+1)] and [(θ_i, γ_i − σ_γi), (θ_i+1, γ_i+1 + σ_γi+1)], where θ is the latitude, σ_γ is uncertainty in spectral index γ, and γ_i < γ_β < γ_i + 1. We solve these two linear equations for when the spectral index γ equals the half maximum value and then difference the θ-intercepts. Halving these differences yields the estimated uncertainty of β. We find β from 2009 to 2013 is 253 ± 210. Note that the uncertainty in γ_β is much smaller than the difference in γ_i and γ_i+1.

For the latter years of IBEX observations (2014–2017), the latitudinal structure of the ENA spectral index is different, reflecting conditions closer to solar maximum (Zirnstein et al. 2017a, 2020). As shown in Figure 5, the spectral index at high latitudes in 2009–2013 is <1.5 and increases to ∼2.5 at low latitudes. However, in 2014–2017, the spectral index at high latitudes is more variable and, on average, is ∼1.8, reflecting a mixture of slow and fast SW during solar maximum. We argue that γ_β found in 2009–2013 is also applicable to 2014–2017, since it should represent the average spectral index of ENAs produced from slow and fast SW, and thus reflects the latitude at which the spectrum is balanced by ENAs produced along distances A and B. If we take γ_β found in 2009–2013 and apply it to the data taken in 2014–2017, we find that β is 172 ± 152. In Section 4 we discuss in more detail the origin of the spectral index observed in 2014–2017.

We note that the derivation of angle β from Equation (1) includes the effects of the differences in PUI energies of the slow versus fast SW plasma in the IHS. Differences in the ENA spectra at low versus high latitudes is direct evidence that the average PUI energy (or temperature) is different in the slow and fast SW, and thus is inherently taken into account in our analysis.

We assume that the distances to the TS R₁(α) ≅ R₂(β). We make this assumption based on the TS geometry results presented by McComas et al. (2019). These authors derived the 3D TS geometry based on Voyager 1 and 2 observations of the TS crossings on the upwind side and magnetic disconnections from the acceleration source of anomalous cosmic rays on the downwind side of the TS. Based on their results, the distance to the TS at points R₁ and R₂ is approximately 160 au for both cases. We note that this distance is an estimate of the time-averaged distance to the TS. It is shown by Equation (1) that if the distance to the TS at the fast SW latitude is larger (R₂ > R₁), then the cooling depth distance will decrease. On the other hand, if R₂ < R₁, then the cooling depth distance will increase. If the dynamic pressure of the SW is largely independent of latitude, as derived from Ulysses observations (e.g., McComas et al. 2008), then the distance to the TS at different latitudes will be similar, as was assumed in our analysis. However, changes over time due to solar transients or during the transition from solar minimum to maximum may change this assumption (e.g., Pogorelov et al. 2017), but only if the changes occur over large enough timescales to affect the LOS-integrated ENA flux measurements. Therefore, for the time-averaged SW near solar minimum, we believe the assumption of R₁ = R₂ is a reasonable approximation. Uncertainties in the distance to the TS over time due to changes over the solar cycle are discussed in Section 4.

With these assumptions, we solve Equation (1) for B, yielding

$$\begin{eqnarray}&&B={R}_{\mathrm{TS}}\left(\displaystyle \frac{\sin \alpha }{\sin \beta }-1\right),\end{eqnarray} \tag{ 3 }$$

where R_TS = R₁ = R₂.

In order to calculate the distance to the ENA source, i.e., cooling length (Schwadron et al. 2014; Zirnstein et al. 2016a) C = A + B, next we need to solve for A. As stated earlier, one possible solution is that the ENA spectral index γ_β shown in Figure 5 is produced when the LOS-integrated fluxes from the fast and slow SW are created along equal distances, i.e., A = B. However, this assumption ignores the fact that the energetic parent protons in the IHS experience an extinction effect from losses by charge exchange. This effect is important to consider since our goal is aimed at deriving the cooling length C which itself depends on the extinction effect by charge exchange. Therefore, one should set

$$\begin{eqnarray}\begin{array}{rcl}{{Ae}}^{-\displaystyle \int {\eta }_{f}\displaystyle \frac{{{dL}}_{f}}{{u}_{f}}} & = & {{Be}}^{-\displaystyle \int {\eta }_{s}\displaystyle \frac{{{dL}}_{s}}{{u}_{s}}},\\ {{Ae}}^{-\left\langle {\eta }_{f}\displaystyle \frac{{L}_{f}}{{u}_{f}}\right\rangle } & \cong & {{Be}}^{-\left\langle {\eta }_{s}\displaystyle \frac{{L}_{s}}{{u}_{s}}\right\rangle },\end{array}\end{eqnarray} \tag{ 4 }$$

where η = n_Hσ(v)v is the charge-exchange rate, n_H is the neutral H density, σ is the energy-dependent, charge-exchange cross section (Lindsay & Stebbings 2005), v is the speed of ENAs, L_f and L_s are the average streamline distances from the TS to the IBEX LOS along A (in fast "f" SW) or B (in slow "s" SW), respectively, and u is the plasma flow speed for the fast or slow SW. The exponential terms represent the probability of survival by charge exchange. We make the approximation that η_f,s and u_f,s are constant over the average streamline distances L_f,s.

To solve Equation (4), first we note that the rates of charge exchange depend on the dominant energy of ENAs from the slow and fast SW affecting the measured spectral index. Since we are calculating the spectral index over energy bins with central energies from 1.74 to 4.29 keV, we can reasonably assume that the measured ∼1.7 keV ENAs are dominated by ENAs created in the cooler, slow SW, and the measured ∼4.3 keV ENAs are dominated by ENAs created in the hotter, fast SW. Even though the slow SW also creates lower energy ENAs (<1.7 keV), these likely do not affect the spectral index calculated over 1.74–4.29 keV as much as ∼1.7 keV ENAs (see IBEX-Hi ESA response functions in Funsten et al. 2009). Therefore, the ratio of the charge-exchange rates in the fast and slow SW η_f/η_s = [σ(v)v]_{4.3 keV}/[σ(v)v]_{1.7 keV} ≅ 1.2.

The typical speed of the slow SW upstream of the TS is between ∼350 and 500 km s⁻¹ and the typical speed of the fast SW is between ∼700 and 800 km s⁻¹ (Figure 4), yielding a ratio range of 1.4 < u_f/u_s < 2.3, and on average u_f/u_s ≅ 1.85. The ratio u_f/u_s will change in the IHS depending on the TS compression ratio. While there are no in situ observations of the TS in the direction of the heliotail, steady-state MHD simulations of the heliosphere suggest that the compression ratios at the upwind and downwind TS are similar (e.g., Zirnstein et al. 2018b; Heerikhuisen et al. 2019). The biggest differences between the upwind and downwind TS are their relative distances from the Sun due to the stagnating flow in the upwind direction and outflow in the downwind direction.

If we assume that the TS compression ratio for the slow SW in the downwind direction is similar to that observed by Voyager 2 (∼2.5; Richardson et al. 2008), and the compression ratio for the fast SW must be greater due to its higher Mach number (e.g., Fitzpatrick 2014), then the ratio of fast to slow SW speed in the IHS becomes smaller. We estimate this ratio by solving the ideal MHD jump conditions at a perpendicular shock. Using Voyager 2 observations, but scaling the proton density down to account for the longer distance to the downwind TS (84–160 au; McComas et al. 2019), we are able to acquire a compression ratio of ∼2.5 after adding a PUI population with 25% of the proton density and temperature of 4.7 × 10⁶ K. By substituting an SW speed 1.85 times faster than that observed by Voyager 2, the compression ratio increases to ∼3.4. With this higher compression for the fast SW, the ratio of fast-to-slow SW in the IHS is most likely 1 < u_f/u_s < 1.7, and on average u_f/u_s ∼ 1.4. Note that the compression ratio is in general determined from the ratio of densities across the TS, and only from the radial velocity component if the shock normal is radial. Since the TS is likely not Sun-centered nor strictly spherical, our analysis above has uncertainties depending on the actual orientation of the shock normal. However, given (1) the large uncertainties in our final results, (2) our derivation of the minimum distance to the ENA source, and (3) that the TS is more likely quasi-perpendicular than not at low-to-midlatitudes, this is likely a reasonable assumption for the purposes of our study.

Finally, the distances L_f and L_s are the average streamline lengths that the plasma takes to travel from the TS through the fast and slow SW, respectively, intersecting the LOS along A + B. Since A is located farther from the TS than B (see Figure 3(b)), it is obvious that L_f must always be greater than L_s. ENAs created along distance B originate from plasma parcels that travel different streamline distances from the TS to different points along B. One might assume that the average streamline distance L_s must connect to B at its halfway distance, such that X  = B/2 in Figure 3(c). However, this is not necessarily true. A generalized solution for L_s can be calculated by assuming that the nearly triangular area between L_s, X, and the TS, Area₁ ≅ (1/2)X (X + τ₁)tanβ, equals half the entire area between min(L_f), B, and the TS, Area₂ ≅ (1/2)B(B + τ₁+ τ₂)tanβ. By solving this condition, a set of equations can be derived to relate L_s and the minimum of L_f,

$$\begin{eqnarray}\begin{array}{rcl}\min ({L}_{{\rm{f}}}) & = & \displaystyle \frac{B+{\tau }_{1}+{\tau }_{2}}{\cos \beta },\\ {L}_{{\rm{s}}} & = & \displaystyle \frac{\sqrt{{\tau }_{1}^{2}+2B\left(B+{\tau }_{1}+{\tau }_{2}\right)}}{2\cos \beta }.\end{array}\end{eqnarray} \tag{ 5 }$$

Solving Equations (5) for the ratio of min(L_f)/L_s yields

$$\begin{eqnarray}&&\displaystyle \frac{\min ({L}_{{\rm{f}}})}{{L}_{{\rm{s}}}}=\displaystyle \frac{2\left(B+{\tau }_{1}+{\tau }_{2}\right)}{\sqrt{{\tau }_{1}^{2}+2B\left(B+{\tau }_{1}+{\tau }_{2}\right)}}\gt \sqrt{2}.\end{eqnarray} \tag{ 6 }$$

Here we show that ${\eta }_{f}\displaystyle \frac{{L}_{f}}{{u}_{f}}\geqslant {\eta }_{s}\displaystyle \frac{{L}_{s}}{{u}_{s}}$ and thus A ≥ B. Knowing the ratios η_f/η_s = 1.2 and u_f/u_s < 1.7, we find that for A ≥ B to be true, L_f/L_s > 1.7/1.2 ∼ $\sqrt{2}$, which has been proven by Equation (6) for the minimum of L_f when we limit A = 0. Thus, as A increases, L_f must increase and L_s remains constant, and the condition L_f/L_s > $\sqrt{2}$ is easily satisfied.

Therefore, instead of solving Equation (4) for A and estimating values for η, L, and u, we can reasonably assume that A ≥ B, yielding C = A + B ≥ 2B as a minimum for the cooling length of ENAs. Thus, C is given as

$$\begin{eqnarray}&&C\ \geqslant 2{R}_{\mathrm{TS}}\left(\displaystyle \frac{\sin \alpha }{\sin \beta }-1\right).\end{eqnarray} \tag{ 7 }$$

The uncertainty in C, σ_C, is calculated by the propagation of uncertainties for α (σ_α) and β (σ_β), yielding

$$\begin{eqnarray}&&{\sigma }_{C}=\displaystyle \frac{2{R}_{\mathrm{TS}}}{\sin \beta }\sqrt{{\cos }^{2}\left(\alpha \right){\sigma }_{\alpha }^{2}+{\left(\displaystyle \frac{\sin \alpha \cos \beta }{\sin \beta }\right)}^{2}{\sigma }_{\beta }^{2}}.\end{eqnarray} \tag{ 8 }$$

The solutions for Equations (7) and (8) are shown in Figure 6, where we set R_TS as a free parameter. We also show the predicted distance to the TS in this particular direction from McComas et al. (2019), i.e., R_TS = 160 au. With this distance, we find that C ≥ 129 ± 35 au in 2009–2013 and C ≥ 329 ± 56 au in 2014–2017. The total distances of the ENA source from the Sun are ≥289 ± 35 au and ≥489 ± 56 au, respectively. The increase in time of the distance to the ENA source is a potential indication of the motion of the ENA source down the heliotail.

Zoom In Zoom Out Reset image size