A Three-dimensional Map of the Heliosphere from IBEX

For this study, we utilize the first 11 yr of IBEX-Hi ENA flux measurements from the five energy passbands centered on 0.71, 1.11, 1.74, 2.73, and 4.29 keV (these are also referred to ESA steps 2–6). The IBEX spacecraft maintains a nearly Sun-pointed spin axis, repointing in ecliptic longitude back toward the Sun every ∼7.5 days during the period between launch and an orbit-stabilizing maneuver in 2011 and then every ∼4.5 days thereafter (McComas et al. 2011, 2012). As IBEX spins, it samples a fixed longitudinal swath of the sky over all latitudes. Each repointing shifts the longitude by 4°–5° to sample the next adjacent swath. Over the course of 6 months, the spin axis rotates through 180° to produce a complete all-sky map. Because of the ∼30 km s⁻¹ motion of the Earth (and therefore IBEX) around the Sun, half of the data in a 6 month map is collected over the half spin when the IBEX spacecraft is moving toward the incoming observed ENAs (Ram data), and half is taken over the half spin when it is moving away from them (anti-Ram data). In the following 6 months, another all-sky map is collected, but now the portion of the sky that was previously sampled in the Ram direction is sampled in the anti-Ram direction, and vice versa.

A consequence of the alternate Ram/anti-Ram sampling is that a given region appears considerably brighter when observed in the Ram direction than when observed in the anti-Ram direction, complicating direct comparison. Therefore, time variation studies often make use of yearly maps, where the Ram portions of two adjacent 6 month maps are combined, resulting in "Ram-only" maps. Similarly, yearly "anti-Ram-only" maps can be constructed. Such maps allow for the most directly comparable measurements of time variation of the whole sky (McComas et al. 2012, 2014, 2017, 2020).

Here we take a different approach in order to maximize the time resolution for applying the sounding method. Specifically, we use 6 month map sets where the ENA fluxes have been Compton–Getting (CG) corrected by transforming the data from the spacecraft frame to the solar inertial frame (McComas et al. 2010). The direct CG correction transforms both the flux and the measurement energies from the spacecraft frame to a solar inertial frame. Because the transformed energies depend on the look angle of the IBEX-Hi sensor relative to the spacecraft velocity vector, the measurement energies in a CG-corrected map vary across the map. We prefer to work with maps having the same measurement energy for all of the pixels; thus, for a given passband, a "monoenergetic" map is created, consisting of fluxes shifted from the direct CG-transformed values to values at a common solar inertial frame energy. For convenience, the energies assigned to the monoenergetic CG maps correspond to the original spacecraft-frame central energy values for each of the passbands. The ENA fluxes are also adjusted for the survival probability of ENAs traveling 100 au from beyond the TS to 1 au (Bzowski 2008) so that they reflect the expected ENA flux levels in the HS. All data used in this study are taken from the latest and most complete release of the IBEX-Hi data (McComas et al. 2020), which incorporated all of these corrections.

As mentioned in the Introduction, applying the sounding method to regions of the sky containing the ENA ribbon is complicated by the fact that IBEX is observing ENAs from two distinct sources, as the GDF originates in the HS and the ribbon emissions are likely coming from the outer HS, beyond the HP. Thus, to apply the sounding method for mapping out the boundaries of the heliosphere, one must either ignore the regions of the sky containing the ribbon altogether or separate and remove the ribbon flux. Because so much of the sky is covered by the ribbon, we have decided to take the ribbon separation and removal approach. A ribbon separation technique introduced by Schwadron et al. (2011) Requires summing multiple maps together to get good enough statistics to accurately identify ribbon versus GDF flux at the level of individual 6° × 6° pixels. Because we wish to carry out ribbon separation on individual 6 month maps in order to preserve as much time resolution as possible, we have devised a different technique.

We summarize the technique here and describe the procedure in detail in the Appendix. We transform the ecliptic-frame ENA maps described above into a "ribbon-centered" frame, where the ribbon center is located at the "north pole" of the maps and the ribbon itself appears as a constant latitude feature centered on a latitude 75° from the pole. In this manner, it is easy to take longitudinal cuts across the ribbon and fit functions approximating the shape of the ribbon and the GDF in order to partition the ribbon flux from GDF. To improve statistics, prior to fitting, we average the flux in longitude every 24° (4 pixels), resulting in 15 longitude strips, each with 30 6° wide bins spanning pole-to-pole.

For each of these strips, we identify a set of latitude bins on either side of the ribbon (but not including it) and fit a second-order polynomial to capture the shape of the GDF beneath the ribbon. The flux associated with this fit is subtracted from the total flux in the ribbon bins, and a Gaussian representing the cross-sectional shape of the ribbon is fit to the remaining flux. For each bin in the ribbon, we ratio the Gaussian fit to the sum of the Gaussian and second-order polynomial fits and use this ratio to partition the total observed flux in the bin between the ribbon and the GDF. For a given latitude, the same ratio is applied to each of the four initial 6° × 6° ribbon-centered pixels composing a bin. What results are two maps, one of the GDF and one of the ribbon flux. The separated maps are then transformed back into the ecliptic frame. Figure 1 gives an example of the separation technique applied to the 2013B maps. See the Appendix for further discussion.

Zoom In Zoom Out Reset image size

The GDF maps in Figure 1 show that our technique successfully removes the ribbon. The pressure enhancement associated with the heliospheric nose appears as a nearly circular structure centered ∼20° south of the upwind direction, owing to the asymmetric pressure of the external interstellar magnetic field (ISMF) draped around the heliosphere (McComas & Schwadron 2014). The enhancement associated with flux coming from the tailward direction appears at the edges of the maps (McComas et al. 2013). For this study, we are not concerned with the ribbon maps, although one could apply the sounding method to these as well to find the distance to the ribbon, which will be the focus of a future study. We acknowledge that applying our ribbon separation technique to 24° wide strips and then using the result to partition fluxes in individual 6° × 6° pixels leads to an under- or oversubtraction of the ribbon flux in a given pixel. However, as we shall see in the next section, the sounding method will be applied to groups of summed 6° × 6° pixels, or "macropixels" (again, to improve statistics); thus, the partitioning error present at the native pixel scale will average out.

Before applying the sounding method to the ENA data, we divide the sky into 56 regions, or macropixels, of roughly the same solid angle of roughly a quarter steradian each. We select macropixels that are large enough to reduce the statistical noise within a pixel to a level where statistics do not significantly influence the time correlation but small enough to capture a reasonable spatial resolution. Figure 2 shows an overlay of the boundaries of the 56 macropixels on a representative ENA map for context. The sky is divided into seven latitude bands of macropixels, which we will respectively refer to as the equatorial pixels (12 pixels between 15°N and 15°S), the north and south midlatitude pixels (12 pixels each spanning 15°N–45°N and 15°S–45°S), the north and south high-latitude pixels (9 pixels each spanning 45°N–80°N and 45°S–80°S), and 2 polar pixels (each extending from 80°N/S to the poles). In longitude, the macropixels are oriented such that the longitude of the heliospheric upwind direction (−105°, 5°) is aligned with the center of a macropixel, which we refer to as the nose pixel. The Voyager 1 and 2 locations are also shown and are located near the centers of two macropixels, which we refer to as the V1 and V2 pixels, respectively.

Zoom In Zoom Out Reset image size

For the top four energies (1.1–4.29 keV), the statistical uncertainties in the macropixels are in the 1%–4% range (1σ), except at the poles, where the uncertainties are ∼0.2%, because the exposure times are nearly continuous. At 0.71 keV, the mean statistical uncertainty is somewhat higher, at 6%, due to the lower sensitivity of the IBEX-Hi instrument at this energy. For a few macropixels where heliospheric viewing is limited by foreground contamination (e.g., the Earth's magnetosphere and magnetosheath), the statistical variation is considerably larger, as discussed in Section 3.5. Note that the macropixel averaging does not remove systematic variations that can impact the time correlation. There are still variations in the data from imperfect subtraction of background and limitations of the accuracy of the CG correction method as noted by McComas et al. (2018a). The cases where statistical and systematic noise do have an impact on the analysis will be noted.

Given an SW disturbance propagating outward at the SW speed v_SW at a particular longitude b and latitude λ, the delay time it takes for the effect of that disturbance to be observed in the return ENA signal is estimated by Equation (3) of Zirnstein et al. (2018),

$$\begin{eqnarray}\,\begin{array}{rcl}\left\langle {t}_{\mathrm{tb}}{\left(b,\lambda \right)}_{i}\right\rangle & = & \displaystyle \frac{{d}_{\mathrm{TS}}(b,\lambda )}{{v}_{\mathrm{sw}}(b)}+\displaystyle \frac{3}{2}\displaystyle \frac{{l}_{\mathrm{HS}}(b,\lambda )}{{v}_{\mathrm{ms}}}\\ & & +\displaystyle \frac{{d}_{\mathrm{TS}}(b,\lambda )+{l}_{\mathrm{HS}}(b,\lambda )/2}{{v}_{\mathrm{ENA}}({E}_{i})},\end{array}\,\end{eqnarray} \tag{ 1 }$$

where $\left\langle {t}_{\mathrm{tb}}{\left(b,\lambda \right)}_{i}\right\rangle$ is the average trace-back time for ENAs originating in the HS and detected in the ith energy passband of IBEX-Hi, d_TS and l_HS are the (unknown) TS distance and HS thickness, v_ms is the average magnetosonic speed in the HS, and v_ENA(E) is the speed of an ENA observed by IBEX at a solar rest-frame energy E.

The first term on the right-hand side is the time it takes for an SW parcel to travel from the Sun to the TS. The second term is the time it takes a pressure pulse generated at the TS traveling at the magnetosonic speed to propagate to the HP and halfway back. This takes into account the finding of McComas et al. (2018a) and Zirnstein et al. (2018) that only after a pressure pulse rebounds from the HP and propagates back to the TS does the plasma pressure fully adjust to a change in the SW dynamic pressure. The reason we use the time it takes for the pulse to only propagate halfway rather than all the way back to the TS is because we are calculating the average trace-back time for ENAs formed between the HP and the TS, and we assume that the midpoint of the HS is the average formation point for ENAs. The last term represents the time it takes ENAs formed at the midpoint to travel back to IBEX.

In this work, the period of observation is a full solar cycle; so, to arrive at realistic values for t_tb, we need to account for changes in the SW speed as a function of latitude over the solar cycle; this is especially important for high latitudes, which have much higher speeds on account of the large circumpolar coronal holes around solar minimum, compared to slower and more variable speeds around solar maximum (e.g., McComas et al. 1998). We therefore calculate time-dependent SW speeds using the method of Sokół et al. (2015, 2020) in which SW speeds over time and versus latitude are derived from interplanetary scintillation (IPS) observations (Tokumaru et al. 2012). Details of the derivation of the SW speed from multistation IPS measurements, which are used in the method by Sokół et al. (2015, 2020), are given in Tokumaru et al. (2010). Note that, as in Reisenfeld et al. (2016), we reduce the outgoing SW speed by 10% from these observed values to account for the deceleration of the SW due to mass loading by pickup ions. For v_ms, we use a value of 314 km s⁻¹, the HS magnetosonic speed derived from Voyager observations by Rankin et al. (2019).

The action of the SW dynamic pressure can be related to the ENA flux observed by IBEX-Hi via the plasma pressure integrated over the LOS through the HS (see Schwadon et al. 2011):

$$\begin{eqnarray}&&{P}_{s}\bullet L=\displaystyle \frac{2\pi {m}^{2}}{3{n}_{{\rm{H}}}}{\int }_{{E}_{\min }}^{{E}_{\max }}\displaystyle \frac{{dE}}{E}\displaystyle \frac{{j}_{\mathrm{ENA}}\left(E\right)}{\sigma \left(E\right)}{\left(\left|v\right|\right)}^{3}.\end{eqnarray} \tag{ 2 }$$

The limits of integration span the energy range of IBEX-Hi CG-corrected observations of the ENA flux j_ENA(E), where E is the ENA energy in the Sun's rest frame, and v is the corresponding velocity. Thus, Equation (2) does not account for the total HS pressure, as there are significant contributions from outside the IBEX-Hi energy range. However, we expect this energy range (0.5–6 keV) to be the range most sensitive to variability in the outgoing SW and its associated pickup ions. Since we are interested in tracking variations and less interested in the absolute pressure, the fact that we are not measuring the total HS pressure is not of concern. The other quantities on the right-hand side of Equation (2) are m, the proton mass; n_H, the interstellar neutral hydrogen density, taken to be 0.127 cm⁻³ (Swaczyna et al. 2020); and σ (E), the charge-exchange cross section (Lindsay & Stebbings 2005).

The quantity P_S denotes the "stationary" pressure, or the internal pressure of the HS plasma if it were at rest in the Sun's rest frame. The distance L is the thickness of the primary ENA-forming region, interpreted here as the thickness of the HS. The product P_S •L is referred to as the "LOS-integrated pressure," and it is determined entirely by observations. As pointed out by Schwadron et al. (2011), the nonstationary component of the LOS-integrated pressure scales more or less linearly with the stationary component; thus, for the purpose of correlating changes in the HS with the 1 au SW dynamic pressure, it is sufficient to work with P_S •L.

In practice, since ENAs are measured in discrete energy passbands, the integral in Equation (2) must be calculated in a piecewise manner. We assume that j_ENA(E_i ) is a measurement of the flux precisely at the central energy E_i of the ith passband, and that the flux between adjacent central energies follows a power law with spectral index γ_i . Then, Equation (2) can be evaluated by the sum of a set of LOS-integrated partial pressures:

$$\begin{eqnarray}\,\begin{array}{rcl}{\rm{\Delta }}{\left({P}_{s}\bullet L\right)}_{i} & \cong & \displaystyle \frac{4\pi \sqrt{2m}}{3{n}_{{\rm{H}}}}{\displaystyle \int }_{{E}_{i}}^{{E}_{i+1}}\displaystyle \frac{\sqrt{E}}{\sigma \left(E\right)}\\ & & \times \,{j}_{\mathrm{ENA}}\left({E}_{i}\right){\left(\displaystyle \frac{E}{{E}_{i}}\right)}^{-{\gamma }_{i}}{dE}.\end{array}\,\end{eqnarray} \tag{ 3 }$$

Note that ENAs detected in the different passbands of IBEX-Hi have different trace-back times. This means that to carry out the time correlation, the energy integral in Equation (2) must be performed by summing over the five partial pressures given by Equation (3), each with a different time offset.

In order to determine the distance to the ENA-forming region, we time-correlate HS plasma pressures calculated using Equations (1)–(3) with 1 au SW dynamic pressure observations. We will now walk through the steps in this procedure by reference to Figure 4, which shows the sounding method applied to a macropixel centered at the ecliptic longitude/latitude of (−75°, 30°). Figure 4 (top panel) shows partial and total LOS-integrated plasma pressures in the HS versus the trace-back date, that is, the IBEX observation date minus the trace-back time: t_IBEX − t_tb. The total pressure given by Equation (3) (black curve) is calculated from the partial pressures (colored curves) for the date range when the trace-back dates for all five IBEX-Hi energy passbands overlap. The total pressure curve is scaled arbitrarily to keep it on the same plot range as the partial pressures.

Zoom In Zoom Out Reset image size

Also shown in Figure 4 (top panel) is the 1 au SW dynamic pressure calculated from the OMNI-2 SW archive (red line; King & Papitashvili 2005), plotted at the actual times of observation. This curve is smoothed by a 1 yr running average to remove annual oscillations due to the ±7° yearly excursion of the Earth in heliographic latitude. The smoothing is also warranted because the propagation time of a pressure pulse through the HS to the HP and back is of order a year; thus, the ENA flux collected along the LOS will be effectively smoothed on this timescale.

The "best" trace-back time is dictated by evaluating the correlation between the total HS plasma pressure and the 1 au SW dynamic pressure for a range of HP distances. Examining Equation (1), we see that the sounding method does not uniquely determine both d_TS and l_HS; thus, we must use an alternate means of fixing one of these. Since d_TS is perhaps better constrained by other methods, we fix d_TS and vary l_HS in the cross-correlation calculation (Figure 4, bottom panel). For the example presented here, we set d_TS = 100 au (the rationale for this will be discussed in Section 3.4) and then vary l_HS between zero and 140 au in 5 au increments. For each value of l_HS, the cross-correlation coefficient is evaluated. Taking the peak cross-correlation coefficient, we find that a best-fit value for l_HS is 35 au (blue diamond in bottom panel), leading to an HP distance of 135 au.

The correlation in this example is of excellent quality, with a high maximum correlation coefficient of 0.94. The results are not always this unambiguous because there is sometimes considerable systematic noise in the ENA data, as discussed in Section 2.3, which can complicate the correlation analysis. Fortunately, the Sun has given us some help. Inspecting the correlation between the total LOS pressure and the SW dynamic pressure, the dominant feature driving the correlation is the rapid rise in the SW dynamic pressure that occurs in 2014. This pressure enhancement later appears in the IBEX-Hi ENA observations starting in late 2016. It first shows up at 4.3 keV, appearing as a strong enhancement ∼30° south of the interstellar upwind direction, and then expands outward in each successive sky map (McComas et al. 2018b, 2019a). The enhancement appears at lower energies later than at the top energies, which is a consequence of the longer travel times for lower-energy ENAs from the HS to IBEX. By the end of 2019, the ENA enhancement has spread to encompass roughly two-thirds of the sky at 4.3 keV (McComas et al. 2020).

In the regions of the sky where the ENA enhancement associated with the 2014 dynamic pressure enhancement is present, it turns out that the cross-correlation results are much more stable pixel-to-pixel if the ENA enhancement is given more weight. This is done in practice by carrying out two correlations, one for the full time series (black line in Figure 4, bottom panel) and one for only the last 5 yr of the time series (green line), where the pressure enhancement dominates. We then take the arithmetic mean of the correlation coefficients (orange line) and assign the best trace-back time to the maximum of the mean correlation. In this manner, we give additional weight to the pressure enhancement but still use information from the full time series to set the HP distance. Note that in this case, it does not matter which correlation is used, since both correlation methods peak at the same choice of HS thickness. As we will see in Section 3.5, using the weighted cross-correlation makes a difference. In regions of the sky that the pressure enhancement has not reached, we simply use the full time series correlation.

We note that the boundaries of the TS and HP change with time to maintain the pressure balance between the SW and interstellar medium (ISM) total pressures. Since the SW pressure is dominated by the dynamic pressure of the outflowing supersonic wind, which evolves throughout the solar cycle, this can cause the boundaries to move. The question therefore arises as to the validity of carrying out a time correlation across 11 yr to determine the distance to the HP if the HP is in constant motion. Time-dependent studies that use realistic inputs for the variability of the SW over the solar cycle predict that the TS distance varies by perhaps ∼±10 au and the HP by even less, ∼±5 au (Izmodenov et al. 2005; Washimi et al. 2011; Izmodenov & Alexashov 2020). As we discuss in Section 3.5, we estimate the uncertainty in the distance to the HP based on the sounding method to be at best ±10 au; thus, the expected temporal motion of the HP is well within the method's uncertainties. Furthermore, as we shall see in Section 4, the derived distance to the HP is fairly insensitive to the choice of TS model (see next section), even though there is considerable variation in TS distances between the models, a variation that is larger than the ±10 au variation of the TS distance over the solar cycle. Therefore, we do not expect the time variation of the TS and HP distances to significantly affect our results.

As mentioned in the previous section, the pressure correlation technique does not uniquely determine both the TS distance and the HS thickness. We therefore have assumed different shapes for the TS and then evaluated the distance to the HP using TS distances derived from these shapes. Another reason for using different TS models is that this will tell us how sensitive the resulting HP shapes are to the choice of TS shape. We carry out the procedure with three models.

• 1.  
  A uniform Sun-centered 100 au radius TS model. This is chosen to demonstrate that the resulting HP distance variations are driven not by any inherent asymmetry in the TS model but rather by the ENA data themselves.
• 2.  
  The spherical TS shape derived by the McComas et al. (2019b) analysis of ACRs measured by the two Voyager spacecraft as they traversed the IHS (the Voyager TS model). As described in the Introduction, this analysis provides a best-fit spherical TS with a radius of 117 au and a center offset from the Sun. Although the true TS shape is not likely to be purely spherical, it does represent the simplest shape that can be fit to four points: the two Voyager TS crossing points and the derived locations of the TS boundary in the tail. This TS model is used because it represents the most data-driven approach to finding the shape of the HP that can be used at this time and also yields an asymmetric TS shape qualitatively similar to that expected from previous IBEX data analyses (e.g., Schwadron et al. 2014; Reisenfeld et al. 2016; Zirnstein et al. 2017a; McComas et al. 2020).
• 3.  
  A TS shape derived from the MHD-plasma/kinetic-neutral model of Zirnstein and Heerikhuisen (the Z-H TS model). The Z-H model solves ideal MHD equations for the SW and interstellar plasmas, treated as a single fluid, and Boltzmann's equation for neutral H atoms, treated kinetically (e.g., Pogorelov et al. 2008; Heerikhuisen & Pogorelov 2010). In their model, SW data from the OMNI database are used to represent SW conditions up to 37° from the solar equator, and Ulysses measurements are used for SW conditions for latitudes above 37°, thus approximately reproducing solar-minimum conditions in solar cycle 23 and yielding a nearly latitude-independent SW dynamic pressure. The local ISM boundary conditions (i.e., ISMF, plasma and neutral density, flow velocity, temperature) are derived from a combination of multi-spacecraft measurements (for more details, see Zirnstein et al. 2021). The rationale for the use of such a model is that this gives us a TS shape driven by detailed heliospheric physics, including kappa-distributed protons in the HS, the draping of the ISMF around the HP consistent with IBEX ribbon measurements, and reproduction of the average TS distances observed by Voyagers 1 and 2 (∼89 au).

The time correlation is carried out for each of the 56 macropixels shown in Figure 2 and repeated for each of the three TS models using the method described in Section 3.3. The best-fit HP distance is taken as that which gives the maximum cross-correlation coefficient. There is a great deal of variation in the appearance of the ENA time series, depending on the region of the sky, and a corresponding variety in the quality of the correlations between the ENA and SW time series. To assess the quality of the correlations, we categorize the correlation behavior of the different sky pixels into four categories. Here we describe each category and look at representative examples. These are all taken from the runs with the 100 au TS model, but what we describe applies generally, as the correlations behave nearly the same for all of our TS models.

Category 1. The first group is the largest category and consists of 21 pixels (38% of the total) having correlations that look very similar to the (−75°, 30°) pixel shown in Figure 4. As shown by the color-coding in Figure 2, these pixels are located in the upwind hemisphere, spanning from the equator to high latitudes (but not the polar pixels). All but two of these pixels are within 75° longitude of the upwind direction. The ENA time series is characterized by a well-defined pressure enhancement, and it closely follows the SW time series across the entire time span. For all of these, the best-fit correlation coefficient is greater than 0.75. Contained in this category are the V1 and V2 pixels, centered at (−105°, 30°) and (−75°, −30°), respectively (Figures 5(a) and (b)). Also in this category is the nose pixel (−105°, 0°). Because of the tightly constrained correlations of the category 1 pixels, a significant drop in the quality of the correlation is readily visible by making quite small time shifts. We therefore set the uncertainty in the best-fit time correlation at the resolution of the ENA time series, or ±3 months (half the time between sky maps). This corresponds to changing the assumed HP distance by ±10 au. We therefore set the uncertainty in the HP distance for category 1 pixels at ±10 au.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fact that the ENA time series consistently track the SW dynamic pressure even up to high latitudes provides strong supporting evidence for the idea that the SW dynamic pressure is a latitude invariant. This is an important assumption for this study, because without it, the ability to apply the sounding method would be limited to only the ecliptic. This assumption also factors into the reconstruction of the SW density profiles as a function of heliolatitude from the SW speeds derived from IPS observation data (Sokół et al. 2015, 2020). Our work indirectly relies on these SW density profiles, as they are used to calculate ENA survival probabilities applied to the IBEX data (Bzowski 2008; McComas et al. 2017, 2020, Sokół et al. 2020).

Category 2. The next category consists of pixels where the correlation is anchored by a well-defined pressure enhancement as in category 1, but the ENA time series does not necessarily cleanly track the SW time series. This category is almost as large as category 1, consisting of 20 pixels (35%). It includes the high-latitude pixels not in category 1, including the polar pixels (see Figure 2). At equatorial/low latitudes, it includes pixels on the flanks of the heliosphere as far as 105° from the upwind direction. Correlation coefficients range from as low as 0.23 up to 0.79. Even though the peak correlations are not as high as for the category 1 pixels, the pressure enhancement still puts a reasonably tight constraint on the HP distance. To quantify this, we look at how much the HP distance must change to reduce the correlation coefficient by 0.1. This varies somewhat from pixel to pixel but is in the range of ±10–20 au; thus, we put the typical uncertainty on category 2 pixels at ±15 au.

Because of their significance in previous work (Reisenfeld et al. 2012, 2016), correlations for the north and south polar pixels are shown in Figures 5(c) and (d), respectively. The statistical uncertainties for these two pixels, at 0.2%, are exceptionally low because the poles are observed continuously throughout the IBEX mission. Note now that the different correlation cases (green and black lines) look quite different from each other, in contrast to Figures 5(a) and (b). The mean of the coefficients (orange line) provides a more stable determination of the HS thickness. Although the correlations with the SW time series are not as tight as they are for category 1 pixels, they are still quite good, with coefficients of 0.65 for the north pole and 0.73 for the south.

The pixels that fall into category 2 are consistently further from the nose than the category 1 pixels. A possible explanation for the drop in the quality of the correlation for category 2 pixels is that because they are more distant from the upwind direction, the HS plasma at these locations is not only responding to pressure signatures radially transmitted across the TS but also to pressure structures convecting tailward within the HS; thus, pressure changes are due to an admixture of these two drivers. Fortunately, the 2014 pressure enhancement is large enough to dominate the internal processing, giving us a strong feature to correlate on.

Category 3. This category consists of pixels in which the 2014 pressure enhancement has yet to show up in the ENA time series. Not surprisingly, this consists exclusively of downwind pixels, numbering 10 in total (18%). They are all equatorial or low latitude and at least 105° in longitude from the upwind direction. The best-fit correlation coefficients are high; all but two are above 0.80, with the lowest at 0.68 and the highest at 0.93. This is in part due to the fact that there is not a lot of structure to correlate on, as the ENA time series for these pixels tends to have a flat downward slope that matches the general trend in the SW time series between 2003 and 2010. That is, the cross-correlation coefficient between two similarly sloped straight lines is intrinsically high.

Figure 5(e) shows the correlation for a typical pixel, centered at (45°, −30°) (located on the downwind, starboard side), having an "HP distance" of 290 au. This pixel has a peak correlation coefficient of 0.930, but we can see that the correlation is good for a very broad range of HS thicknesses for the reason above; thus, it is not clear how well constrained the distance is. We place "HP distance" in quotes because we now have to modify our interpretation of this distance when applying it to the downwind portion of the heliosphere. If the HP truly does truncate the heliosphere approximately at this distance in the tail, then this interpretation is fine. However, as we noted at the beginning of Section 3, due to the finite lifetime of the HS plasma parent to the observed ENAs, this must really be interpreted as a lower limit on the distance to the HP. In principle, the heliotail could extend hundreds of au downwind and we would never know it. Strictly speaking, the distance we have derived is simply the distance to the far edge of the ENA-forming region for ENAs in the IBEX-Hi energy range.

Figure 5(f) shows the directly downwind (antinose) macropixel at (75°, 0°), which has a peak cross-correlation coefficient of 0.77 occurring at a distance of 380 au. This value probably has a large uncertainty, but that said, given the overall shape of the ENA time series, its phasing with the SW time series strongly suggests that the ENA source must be of order 350+ au downwind. These broad correlation ranges are typical for category 3 pixels, and we place the uncertainty on the distance to the edge of the ENA formation region at ±30 au.

Category 4. The final category consists of a small number of macropixels where the ENA time series is particularly noisy, and for this reason, the cross-correlation analysis has a very high uncertainty. For these pixels, the cross-correlation is flat across a wide range of HS thickness values, leaving the choice of HS thickness ambiguous. Five pixels (9%) fall in this category, located in a part of the sky where IBEX's view of the heliosphere is partially obscured by the Earth's magnetosphere and magnetosheath; thus, there is limited usable data. This band, between ecliptic longitudes 120° and 150° on the starboard flank and 0° and 30° on the port flank, contains large amounts of statistical variability, which is further amplified by the CG correction. As a result, much of the structure in the ENA time series is noise, leading to false correlations. In these cases, we are forced to manually adjust the choice of HP distance based on patterns in the ENA time series from adjacent longitudes or the arrival of the pressure enhancement in the last few high-energy points. Figure 5(g), for the (135°, 0°) pixel, clearly shows the high level of statistical fluctuation in the ENA partial pressures, as well as a very flat correlation coefficient curve across a large range of HS thicknesses. The blue circle in the lower panel shows the calculated peak correlation at l_HS = 195 au, and the blue diamond shows the manually selected HS distance, l_HS = 100 au.

A list of the distances to the TS and HP for each of the 56 pixels and three different TS models, as well as the correlation coefficients for the given choice of HS thickness, is available in the online-only material as an Excel spreadsheet. The next section describes the shape of the heliosphere based on these distance determinations.

Although the shape of the HP is nearly independent of the choice of TS model, this does not mean that there are no important model-dependent differences to the heliosphere structure, particularly with regard to the thickness of the HS.

We first consider the 100 au TS model. As stated in Section 3.4, this model was chosen to demonstrate that asymmetries in the derived HP were not artifacts of the choice of TS model. In a sense, the 100 au TS acts as a "control" for the other models. Since it is not a particularly realistic TS model, we should not place much weight on conclusions drawn from it about the shape of the HS. For instance, in the region south of the nose where the minimum HP distance is 115 au, the HS appears to be about 15 au thick, which would seem unreasonably thin in light of pretty much any theoretical model of the HS. We know that since this region includes the pixel containing the Voyager 2 TS crossing, the TS distance in this region is more like 80–90 au, so the 100 au TS leads to an underestimate of the HP thickness there.

By design, the Voyager-based TS model (McComas et al. 2019b) accurately depicts the TS distance at the Voyager TS crossings, which means it likely represents the TS distance fairly well in regions near the nose. However, the four points defining this spherical model are all within 30° of the ecliptic, so the TS is not nearly as well constrained at high latitudes. Indeed, at high northern latitudes, the model predicts a TS distance of ∼150 au, which is likely an overestimate, as this leads to an HS thickness of 30–40 au, which for this model is as thin as it is around the nose, where we expect the HS to be at its thinnest owing to the nose being where the ISM flow and field exert maximum pressure. This result is also inconsistent with numerical models (e.g., Pogorelov et al. 2013; Izmodenov & Alexashov 2020). To be clear, McComas et al. (2019a, 2019b) did not claim that the TS is strictly spherical; rather, they fit an offset spherical shape on purely empirical grounds, as this is the most complex shape constrained by a four-point fit, which is the number of observed points.

Turning to the Z-H TS model, the shape of the TS is largely independent of latitude because the distance to the TS is determined by SW dynamic pressure, which is approximately latitude-independent, even though the SW speed varies significantly with latitude. Even during solar maximum, when the fast SW disappears, the Z-H TS model predicts a similar TS shape to solar minimum, so long as the SW dynamic pressure does not significantly change over time. Moreover, this model is constrained to reproduce the averaged TS crossing distances from Voyagers 1 and 2 (∼89 au), but it would need to incorporate time dependence to explain the observed ∼10 au asymmetry. This model provides a representation of the TS shape based on our best understanding of heliospheric physics. (Note that between the various heliospheric modeling research groups, although there is a wide range of predicted shapes for the HP, there is far less variation in the modeled TS itself.) Figures 8 and 9 show that the resulting HS is thinnest near the region of maximum compression by the LISMF and thicker over the north pole than the south and flares symmetrically about the port and starboard flanks of the TS.

Our global depiction of the heliosphere intersects with direct observation in two regions of the sky, in the pixels containing the Voyager 1 and 2 TS and HP crossings. This provides a means of validating the sounding method, at least for the upwind hemisphere. Table 1 gives the distances to these boundaries and the thickness of the HS for the three versions of the TS model and compares them to the distances determined from the Voyager 1 and 2 crossing distances. Note that the values for the Voyager model TS distances do not exactly match the actual Voyager distances because these pixels are averages over a 30° × 30° area of the sky.

The last column shows that the sounding method accurately determines the distance to the HP for the V2 pixel to well within the uncertainty of the method (±10 au). The HS thickness determined using the Voyager and Z-H TS models of 30 au is also quite close to the observed HS thickness of 35 au. It is no surprise that the 100 au model underestimates the HS thickness, since the TS distance is an overestimate. For the V1 pixel, the HP distances are consistently higher than the observed Voyager 1 HP crossing distance but are still at or near the uncertainty of the method. One must also keep in mind that the Voyager crossings only indicate the boundary locations at one point in time, and that the TS and HP boundaries have quite likely moved relative to one another by order of 5–10 au between their respective crossings by Voyager.

On the whole, the comparison to the Voyager distances validates the sounding method for determining accurate distances to the HP. Coupled with what we expect are realistic models of the TS, it is possible to determine a reasonably accurate measure of the thickness of the HS as well.