Oblique and rippled heliosphere structures from the Interstellar Boundary Explorer

Past analysis has shown that the heliosphere structure can be deduced from correlations between long-scale solar wind pressure evolution and energetic neutral atom emissions. However, this required spatial and temporal averaging that smoothed out small or dynamic features of the heliosphere. In late 2014, the solar wind dynamic pressure increased by roughly 50% over a period of 6 months, causing a time and directional-dependent rise in around 2–6 keV energetic neutral atom fluxes from the heliosphere observed by the Interstellar Boundary Explorer. Here, we use the 2014 pressure enhancement to provide a simultaneous derivation of the three-dimensional heliospheric termination shock (HTS) and heliopause (HP) distances at high resolution from Interstellar Boundary Explorer measurements. The analysis reveals rippled HTS and HP surfaces that are oblique with respect to the local interstellar medium upwind direction, with significant asymmetries in the heliosphere structure compared to steady-state heliosphere models. We estimate that the heliosphere boundaries contain roughly ten astronomical unit-sized spatial variations, with slightly larger variations on the HTS surface than the HP and a large-scale, southwards-directed obliquity of the surfaces in the meridional plane. Comparisons of the derived HTS and HP distances with Voyager observations indicate substantial differences in the heliosphere boundaries in the northern versus southern hemispheres and their motion over time.

Past analysis has shown that the heliosphere structure can be deduced from correlations between long-scale solar wind pressure evolution and energetic neutral atom emissions. However, this required spatial and temporal averaging that smoothed out small or dynamic features of the heliosphere. In late 2014, the solar wind dynamic pressure increased by roughly 50% over a period of 6 months, causing a time and directional-dependent rise in around 2-6 keV energetic neutral atom fluxes from the heliosphere observed by the Interstellar Boundary Explorer. Here, we use the 2014 pressure enhancement to provide a simultaneous derivation of the three-dimensional heliospheric termination shock (HTS) and heliopause (HP) distances at high resolution from Interstellar Boundary Explorer measurements. The analysis reveals rippled HTS and HP surfaces that are oblique with respect to the local interstellar medium upwind direction, with significant asymmetries in the heliosphere structure compared to steady-state heliosphere models. We estimate that the heliosphere boundaries contain roughly ten astronomical unit-sized s pa ti al variations, with slightly larger variations on the HTS s ur fa ce t h a n t he H P and a large-scale, southwards-directed obliquity of the surfaces in the meridional plane. Comparisons of the derived HTS and HP distances with Voyager observations indicate substantial differences in the heliosphere boundaries in the northern versus southern hemispheres and their motion over time.
The heliosphere surrounding our solar system is formed by the interaction between the solar wind (SW) and the partially ionized, local interstellar medium (LISM) 1 . The interstellar plasma, consisting mostly of H and He, is slowed at the bow wave upstream of the heliosphere 2,3 and diverted around the heliopause (HP) [4][5][6][7] . Interstellar neutral atoms, however, can cross the HP and enter the heliosphere [8][9][10] . Low energy interstellar neutrals are detected directly by the Interstellar Boundary Explorer (IBEX) 11 near Earth 12,13 and Ulysses GAS 14,15 , but they also may undergo charge exchange collisions inside the heliosphere. The ionization of interstellar neutrals in the supersonic SW and inner heliosheath (IHS) produces energetic pickup ions (PUIs) that dominate the plasma pressure. Through another charge exchange collision, PUIs create energetic neutral atoms (ENAs) at energies much greater than the interstellar neutrals. IBEX measures ENA fluxes at energies up to roughly 6 keV from all directions of the sky and has accumulated more than a solar cycle of ENA observations since 2009 (ref. 16 ).
It has become clear over the past decade that the heliosphere can respond globally to large-scale changes in the SW dynamic pressure. Voyager observations within the IHS have shown large variability in the magnetic field, thermal ion properties and transients propagating across the IHS and into the LISM [17][18][19] , as well as the dynamic characteristics of the heliosphere boundaries 5,6,20,21 , but only along their Article https://doi.org/10.1038/s41550-022-01798-6 With the only two in situ measurements of the heliosphere boundaries from the Voyager spacecraft 4,6,21 , as well as determinations of the heliospheric termination shock (HTS) structure on the flanks from the Voyagers' magnetic disconnection events in the IHS 27 , the heliospheric community has realized the importance of using ENA imaging to detect and correlate changes in the SW with ENA emissions across the sky. Reisenfeld et al. (ref. 28 ) recently demonstrated how IBEX measurements can be used to map the three-dimensional heliospheric structure on large scales using a combination of SW observations at 1 au and global simulations that inform us of the behaviour of ENA emissions in the IHS. With the intention of studying the time-averaged shape of the HP, the methods used by Reisenfeld et al. allowed for the respective trajectories. IBEX, with its ability to map the entire sky every 6 months, has revealed both gradual, long-term changes in ENA fluxes 16,22,23 and abrupt, short-term variability linked to changes in the SW dynamic pressure 24,25 . A large increase in SW dynamic pressure observed by ACE and Wind in late 2014 at 1 au ( Fig. 1) was reflected in enhanced ENA emissions measured by IBEX beginning in late 2016. Increased ENA fluxes were first seen roughly 30° below the nose of the heliosphere 24 (that is, the LISM upwind flow direction), followed by enhancements over larger regions of the sky later in time 16,25 . The spatially dependent response of heliospheric ENAs to the SW pressure change were shown to be caused by the asymmetric structure of the heliosphere 26 .  16 ). a, Ram maps correspond to times when observations are made in the spacecraft ram frame. b, Antiram (A-ram) maps are observations made in the antiram frame. Pixels in corresponding years between rows a and b are offset by 6 months. IBEX observations are corrected for the Compton-Getting effect when transforming from the spacecraft frame to the solar inertial frame and corrected for ENA losses between 100 and 1 au. c, SW dynamic pressure observed by ACE and Wind at 1 au in the ecliptic plane (black), smoothed over two CRs. Approximate time delays between SW and ENA observations are illustrated by the coloured grey bars. The running linear slope fit to SW pressure over ±3 CRs is shown in red. d, IPS observations cover the three CRs nearest to the peak change in SW dynamic pressure (CR 2,154-2,156, large black dots), centred on 2014.75. e, IPS-derived SW speeds as a function of heliolatitude during CR 2,154-2,156. Speeds are shifted uniformly to match OMNI at low latitudes during each CR (OMNI, coloured points). We use SW observations weight averaged over this period (black curve) to analyse IBEX observations. The grey contour represents the propagated standard deviation of the average. f, A large coronal hole in the southern hemisphere, visible in SDO/AIA observations as the dark colour spot, resulted in fast SW at mid-latitudes in CR 2,156 (image courtesy of NASA/SDO and the AIA science team).
Article https://doi.org/10.1038/s41550-022-01798-6 estimation of the HP boundary at large scales over nearly the entire sky. However, the boundaries of the HTS and HP are expected to move on the order of 10 au over a solar cycle 29,30 , or perhaps more if Rayleigh-Taylor and Kelvin-Helmholtz instabilities of the HP surface are strong and prevalent 31,32 . These variations are within the uncertainties of the time-averaged model demonstrated by Reisenfeld  A substantial advancement of the current study compared to previous analyses is the determination of the HTS shape directly from IBEX observations, without assuming the HTS shape a priori from preexisting models. We use observations of a single, global increase in SW dynamic pressure in late 2014 and two separate temporal features observed in IBEX ENA fluxes as they respond to this SW pressure event between 2016 and 2019 to derive the shape of the HTS and HP over a roughly 2 year time span. This methodology allows us to derive the HTS and HP surfaces at higher resolution than previous analyses, but it can only be applied to directions in the sky where ENA emissions respond strongly to the solar event, that is, where the IHS is closest to the Sun. Thus, our analysis is confined to the half of the sky centred on the direction where ENAs first responded to the global pressure event.
We use IBEX observations of roughly 1.4-6 keV ENA fluxes from 2014 to 2019 in our analysis. As IBEX orbits around Earth, it spins along a Sun-pointed axis allowing it to map the sky every 6 months. Data are collected in the spacecraft frame of reference as the Earth orbits the Sun, both in its 'ram' reference frame where it is moving towards the ENA source and the 'antiram' reference frame where it is moving away from the ENA source 33 . The data are transformed into the solar inertial frame by correcting for the Compton-Getting effect 34,35 and are corrected for ENA losses within roughly 100 au of the Sun 16 . Because the IBEX ribbon overlaps a substantial part of the globally distributed flux (GDF) near the upwind hemisphere 12 , most studies of the GDF require removal of the ribbon feature by using a combination of subtraction, masking, interpolation and reconstruction 28,[36][37][38] . However, our analysis does not require removal of the ribbon because the GDF signal that is rapidly changing in response to the SW pressure increase is much stronger than the slowly varying ribbon flux behind it at the electrostatic analyser (ESA) energy steps examined here. This is primarily due to the longer line of sight (LOS) thickness of the ribbon source region outside the HP 39,40 .
The spatially dependent response of heliospheric ENAs to the SW pressure increase indicates the asymmetric structure of the heliosphere boundaries, as was demonstrated by a global magnetohydrodynamic (MHD) simulation 26 . The timing of the ENA response is correlated with the time for magnetosonic wave propagation from the HTS to the HP and approximately halfway back (that is, near the middle of the ENA source region), which was interpreted as a measure of the time it takes the IHS to respond to global changes in SW pressure. This correlation was used to estimate the distance to the HP across the entire sky, averaged over a solar cycle 28 . On further analysis of global MHD simulations, however, we find that the response of ENAs in the IHS to a strong, global pressure change such as that which occurred in late 2014 (roughly 50% increase in dynamic pressure) is driven by both the magnetosonic wave speed and the flow advection speed, and the point in time where the pressure wave reflected from the HP interacts with the higher-pressure advecting flow somewhere in the IHS. This process is demonstrated in Fig. 2, which shows an illustration of the advection of the high-pressure plasma and travelling magnetosonic wave in the IHS. First, the high-pressure wave front released from the Sun in late 2014 has travelled halfway to the HTS by early 2015 (Fig. 2a) and reaches the nose-ward HTS in mid-2015 (Fig. 2b). After reaching the HTS, a pressure wave travelling at the fast magnetosonic speed is released and locally heats the thermal plasma as it travels through the IHS but does not yet noticeably increase ENA emissions. The pressure wave reaches the HP in late 2015/early 2016, after which a reflected wave travels back towards the HTS. By mid-2016, the advecting flow and reflected wave have met near the middle of the IHS and crossed each other: their interaction results in an increase in ENA emissions by adiabatic heating of the advecting plasma, which is observed roughly 6 months later at 1 au as increased intensities of roughly 4 keV ENAs. This process is demonstrated by simulations shown in Fig. 2. Assuming that the simulation can be used to predict how the heliosphere and ENA emissions qualitatively behave to a rapid increase in SW dynamic pressure, we can derive the distance to the HTS and HP directly from IBEX observations.

Results
Using the fact that the SW pressure increase observed in late 2014 was probably a global event 16 , we first identify the time in IBEX observations that ENAs began responding to the SW pressure change for each pixel in the sky. We limit our analysis to pixels within 90° of the approximate direction from which the ENAs first responded (255°.7, −27°), beyond which the ENA response has either not yet been observed or is too weak to identify. Figure 3 shows examples of ENA fluxes in different directions of the sky for ESA 4-6 and the substantial rise in ENA flux in response to the SW pressure increase. The response of ENA fluxes from the heliosheath is strongest for energies roughly 3-6 keV and becomes weaker at lower energies, although still visible down at energies in ESA 4 (1.4-2.5 keV). (c,f,i). A cubic spline is fit to IBEX data (blue curve with propagated uncertainty contour), and a running slope over ±0.5 years is calculated at each point of the spline (red dashed curve: note that the slope is normalized to the y axis range). We first determine a time range surrounding the time of maximum slope in ENA flux. This range is bounded by the times when the slope is 25% of the peak (red vertical dashed lines). The mean ENA response time is defined as the point of maximum slope within this range (red vertical solid line). However, if multiple peaks in the slope exist within this range (for example, in a), the mean response time is defined as the middle of the range. The initial ENA response time is found at the point of zero slope in ENA flux (orange vertical solid line) that occurs before the mean response time. Note that the uncertainties in the South Pole are smaller primarily due to the higher exposure time per pixel at the poles.
Article https://doi.org/10.1038/s41550-022-01798-6 We use cubic spline interpolation to interpolate ENA fluxes between IBEX data points with 0.01 year resolution and calculate the local linear slope of the spline over a running ±0.5 year window. The point of maximum slope signifies the time when the response of ENA fluxes to the SW pressure is changing the most rapidly. This point in time is approximately when the line-of-sight integrated ENA emissions has reached roughly 50% of its maximum, hereafter called the 'mean ENA response time'.

ENA response times
The moment in time before ENA fluxes first begin to rise is used to identify the inner boundary of the ENA emission region, that is, the HTS. As demonstrated in the simulation (Fig. 2), after the pressure front crosses the HTS and begins propagating through the IHS, the HTS moves outwards slightly and the line-of-sight integrated ENA flux decreases slightly right before any large increase in ENA emissions occurs. Therefore, we identify the time when the ENA flux is at a minimum just before the sharp rise as the 'initial ENA response time'. We identify the initial and mean ENA response times for all available pixels in the sky, which represent approximately 34-37% of the full sky area depending on the ESA, as shown in Fig. 4. We note, however, that there is a potential issue in using this time as the location of the HTS. The time at which we identify the higher-plasma pressure has reached the HTS is probably after the HTS has already begun moving outwards before observing any significant increase in ENA intensity. Therefore, we must interpret this as the maximum distance to the HTS compared to its state before modification by the higher-pressure plasma. The uncertainties in the initial ENA response time are significantly larger than the mean response time due to the variability in ENA flux that occurs before the heliosphere responds to the global SW pressure event.
After identifying the initial and mean ENA response times for ESA 4-6, shown in Fig. 4, we calculate the distance to the HTS, mean ENA source and HP using the steps described in the Methods and summarized here: first, we calculate the distance to the HTS by integrating (1) the time that the SW travels from the Sun to distance r and (2) the time that ENAs within each ESA travel from distance r back to Earth, until the total time passed equals the initial ENA response time observed by IBEX. We use measurements of SW speed, density, temperature and magnetic field from the OMNI database (in-ecliptic SW) 41  The initial (a,e,i) and mean (b,f,j) response times and their corresponding uncertainties (initial response uncertainties in panels c,g,k and mean response uncertainties in panels d,h,l) are shown for each pixel in the sky accepted in the analysis, separately for ESA 4 (a-d), ESA 5 (e-h) and ESA 6 (i-l). Uncertainties for the initial and mean ENA response times are calculated by propagating the uncertainty of each IBEX data point. An additional uncertainty in the initial response time is included by calculating the effect that fluctuations in ENA flux before the initial response time may have on the result. See Methods for more details.
Article https://doi.org/10.1038/s41550-022-01798-6 multi-fluid equations for conservation of mass, momentum and pressure of the SW protons, alphas, H + PUIs and He + PUIs, with source terms for interstellar H and He neutrals ionized by charge exchange and photoionization [43][44][45] , separately deriving distances to the HTS, r HTS , for each pixel. The calculation is weighted by the product of the instrument response function 46 and ENA source spectrum derived from the GDF 38 for ESA 4-6.

Heliosphere boundary distances
The distances to the mean ENA source, r ENA , and HP, r HP , are calculated simultaneously using the difference in time between the initial ENA response time, t HTS , and the mean ENA response time, t ENA , which is the time it takes for the advecting high-pressure plasma to cross the reflected wave over enough distance that the line-of-sight integrated ENA flux reaches half the final pressure state, as illustrated in Fig. 2. All flow advection and wave speeds are calculated from the coupled multi-fluid transport equations, which are advected across the HTS using the single-fluid shock adiabatic equation. Figure 5 shows the results of our analysis, after culling pixels where the ENA fluxes did not show clear behaviour related to the SW event, the uncertainties were too high or there were data gaps. Uncertainties in the distances are calculated by propagating the uncertainties of multiple variables through the analysis. Detailed descriptions of the distance derivation, data culling and uncertainty propagation procedures are provided in the Methods section.  Fig. 4 and weight averaged over IBEX energy passbands 4-6: to heliosphere termination shock (a), to mean ENA source (b) and to HP (c). d-f, Propagated uncertainties of the raw distance results: to HTS (d), to mean ENA source (e) and to HP (f). g-i, A 10° statistical smoothing is performed to fill in gaps for illustrative purposes: to HTS (g), to mean ENA source (h) and to HP (i). A minimum of three nearby pixels within 10° (measured from the pixel centres) is required for filling; otherwise, the data is not modified and existing data are not changed. Similar to previous analyses 28,47 , the closest positions of the heliosphere boundaries are not centred on the LISM upwind direction but rather a few tens of degrees below the nose of the heliosphere (Fig. 5g-i). Cross sections in the ecliptic and meridian planes (Fig. 6) show the boundary surfaces are highly oblique in the meridian plane with respect to the nose, tilted roughly 30° southwards from the ecliptic plane. There is evidence for spatial variations of the boundary distances over angular scales roughly 10° and larger, as is evident in the 'wavy' structure of the surfaces between neighbouring pixels. We estimate the significance of these variations by performing a minimized fit of a quadratic polynomial to the surfaces in the ecliptic plane (both xy and radius-longitude (r-ϕ)) and meridian plane (both xz and radiuslatitude (r-θ)) and find a statistically significant standard deviation ranging between roughly 3 and 10 au in the ecliptic plane and roughly 5-16 au in the meridian plane (ranges indicate ±1-sigma uncertainty ranges). While it appears that there are larger variations in the HTS surface within these planes, our analysis cannot determine them to be statistically significant compared to variations in the HP surface. Spatial variations in the heliosphere boundary distances on the order of 10 au are probably too large to be caused purely by differences in measurement time of ≲1 years; therefore, they may be signatures of persistent ripples or fluctuations along the heliosphere boundary surfaces. Table 1 shows examples of distances to the heliosphere boundaries in several directions of the sky, with comparisons to Voyager observations. Voyager 1 and 2 measurements of the HTS distance from the Sun, although separated in time by roughly 3 years, show an asymmetry of roughly 10 au. Distances derived from IBEX observations taken approximately 10 years later show a larger asymmetry of roughly 25 au, but with a large uncertainty of roughly 17 au. IBEX measurements of the distance to the HTS in the Voyager 1 and 2 directions are separated in time by roughly 0.5 years, but it is unlikely that the asymmetry reported here can be explained by motions over less than 1 year. The observed asymmetry may potentially be linked to (1) north-south asymmetries in the SW mass flux 48 , where SOHO/ SWAN observations of back-scattered Lyman-α radiation suggest the existence of higher SW mass flux and/or dynamic pressure in the northern hemisphere in 2014 compared to the southern hemisphere, which might create an asymmetric heliosphere shape, in contrast to the SW mass flux observed in late 2003 relevant to Voyager 1's HTS crossing 48 , or (2) the pressure exerted by the interstellar magnetic field on the southern hemisphere of the heliosphere [49][50][51] . Global, three-dimensional models of the SW-LISM interaction with dynamic SW boundary conditions have suggested that substantial distortions of the HTS surface might occur over the course of a solar cycle 30,32,51-54 , but the large asymmetries reported here, if statistically significant, have yet to be reproduced by any model.

Discussion
The distances to the HP in the Voyager directions as observed by IBEX are intriguing and potentially controversial. The analysis suggests that the distance to the HP in the Voyager 1 direction is r HP = 131 ± 9 au as observed in 2016.6. This result, while appearing farther than 122 au where Voyager 1 crossed the HP, is still consistent with the fact that Voyager 1 crossed the HP in late 2012 and remained outside ever since 4,5 . In 2016.6, Voyager 1 was 136 au from the Sun, and therefore slightly outside the HP derived from IBEX observations. This suggests an increase in the distance to the HP in the few years after Voyager 1 crossed into interstellar space. We also note that Reisenfeld et al. (ref. 28 ) derived a similar distance to the HP near Voyager 1, although using temporal correlations over a solar cycle, indicating that their result was largely driven by the 2014 SW pressure event.
The distance to the HP in the Voyager 2 direction derived from our analysis is r HP = 103 ± 8 au as observed in 2015.9. At this time, Voyager 2 was 109 au from the Sun and it did not cross the HP until late 2018 at a distance of 119 au (refs. 6,7,55 ). Our results are consistent with Voyager 2 measurements within 1-sigma uncertainty. However, if the HP was as close as roughly 111 au from the Sun in 2015.9, then the HP must have then moved outwards after 2015.9 before Voyager 2 crossed it in 2018. 85. Dynamic heliosphere simulations qualitatively show this outward-moving behaviour of the HP 51,56 , although we must point out that nearly all models have difficulty reproducing Voyager measurements quantitatively. We note, however, that while we have attempted to include all known uncertainties in our analysis, such as the SW speed uncertainty (Methods), potential unquantified variables may contribute to these results.
IBEX has operated successfully and made numerous discoveries over the past 13 years. Using IBEX observations, this study provides high-resolution maps of the heliosphere's HTS and HP surfaces and their spatial variations (Fig. 7). While it is expected that IBEX will continue operating and taking measurements for the near future, a new NASA mission planned for launch in 2025, called the Interstellar Mapping and Acceleration Probe (IMAP) 57 , will improve on IBEX's capabilities by measuring ENA fluxes over a larger energy range with greater accuracy and temporal resolution. IMAP is equipped with three neutral atom imagers, IMAP-Lo, IMAP-Hi and IMAP-Ultra, which will measure neutral atom fluxes from 0.005 to 1 keV, 0.4 to 16 keV and 3 to 300 keV, respectively. With their greater sensitivity, the IMAP ENA imagers will be able to produce full sky maps every 6 months and partial sky maps every 3 months, allowing us to quantify variability in the outer heliosphere at twice the cadence of IBEX. Moreover, IMAP will orbit around L1 and thus not be affected by Earth's magnetosphere. Finally, uncertainties in the SW speed as a function of latitude (for example, discrepancies between  Fig. 5a-c. l IHS is the thickness of the IHS calculated as the difference between r HP and r HTS . Uncertainties (σ r ) are derived from uncertainty maps in Fig. 5d-

Data selection and initial processing
We analyse IBEX-Hi observations of ENAs measured within ESA energy passbands 4-6 (with a full-width at half-maximum from 1.4-2.5, 2.0-3.8 and 3.1 to 6.0 keV, respectively) starting from 2014 as part of data release 16 (ref. 16 ). ENAs measured in ESA passbands 4-6 have the highest signal to noise ratio due to the high rate of transmission through the instrument compared to lower energy passbands 46 . These ENA fluxes also show the quickest and strongest responses to the SW pressure increase, which was less noticeable at energies ≲2 keV (refs. 24,25 ). While ESA 4 fluxes show the weakest response to the SW pressure event (Fig. 3), our analysis accounts for this by yielding higher uncertainties in identifying the timing of the event (Fig. 4). ENA fluxes are observed in the spacecraft 'ram' frame (as IBEX is moving towards its look direction) and 'antiram' frame (as IBEX is moving away from its look direction), covering the sky every 6 months. We use data transformed into the solar inertial frame and corrected for ENA losses between 1 and 100 au, which removes effects of losses due to ionization close to the Sun. Each pixel in the sky has a unique time of observation as Earth orbits around the Sun each year, where IBEX starts taking observations for each ram map in the beginning of each year near 180° longitude and each antiram map starts near 0° longitude. IBEX fills in the sky over time with increasing longitude over the course of 6 months. Before analysis, we apply a smoothing to each pixel by calculating the statistically weighted average of all pixels within 9° and applying the average to the centre pixel. This process smooths fluctuations between closely neighbouring pixels that may be a by-product of imperfect background subtraction that is performed independently for each IBEX orbital swath 33,58,59 . The smoothing also improves the capability of our analysis on deriving time delays from IBEX observations. Note, however, that because the pixels near the poles have smaller solid angle areas, spatial smoothing will inherently combine more pixels that are observed at substantially different times throughout the year. Therefore, we limit the spatial average to pixels within 9° that have measurement times within 0.25 years of the centre pixel. Because IBEX constructs all sky maps over a period of 6 months, the front half of the sky for ram measurements is constructed over the first half of the year and the back half of the sky is constructed over the second half. Because of this, data on either side of ecliptic longitude roughly 180° in ram maps are separated by 1 year in time, and data on either side of ecliptic longitude roughly 0° in antiram maps are separate by 1 year in time. Therefore, smoothing is not applied across ecliptic longitudes 180° and 0° for ram and antiram data, respectively.
Next, we apply an initial culling of the data before our analysis. First, we remove all pixels more than 90° from (255°, −27°), which is the approximate location in the sky when heliospheric ENAs first responded to the late-2014 SW pressure increase 24 . We only analyse pixels in this half of the sky because, for most observations outside this region, there has not yet been a substantial response in ENA flux to the SW pressure increase and therefore making the derivation of heliospheric distances not currently possible. Second, we remove any pixels where there are data gaps at any point in 2014-2019. After this culling, 877 of pixels in the sky remain (or roughly 48.1% of the area of the sky) out of a possible total of 1,800. We note that certain sections of the sky may have culled pixels next to unculled pixels. For example, there is a patch of culled pixels near Voyager 2 (Fig. 5a-c) that indicate potential issues in the data near that region of the sky. Therefore, extra care should be taken when interpreting results from these regions.

Calculation of initial ENA and mean ENA responses
ENA fluxes from the outer heliosphere respond to the large increase in SW dynamic pressure a few years after in-ecliptic spacecraft first observed the SW pressure increase in late 2014. The ENA response is identified by an increase in ENA flux occurring over roughly 1-2 years. Since heliospheric ENAs cannot originate closer than the HTS, the initial rise in 3-6 keV ENA fluxes is used to identify the time at which ENAs first reacted to the SW pressure increase as it crossed the HTS. As the ENA flux continues to rise over time, the rate of increase maximizes and then gradually stops increasing. We identify the middle of this time period as the mean of the ENA source region in the IHS, as described below.
The 'initial ENA response' (hereafter referred to as t HTS ) and 'mean ENA response' times (t ENA ) are identified first by performing a cubic spline interpolation of the ENA flux in each non-culled pixel after 2014 with a high temporal resolution (see examples in Fig. 3). The uncertainty of the spline interpolation is calculated by propagating the data uncertainties (Calculation and propagation of uncertainties section). In the next step, the local linear slope of the spline interpolation is calculated by fitting a line to the spline and uncertainties using least-squares minimization over a ±0.5 year window. A 1-year-wide window is chosen since it represents the time over which IBEX makes at least three observations. This yields the local slope of the interpolated ENA fluxes, as shown by the red dashed curves in Fig. 3. We then find the times at which the slope reaches 25% of the peak slope in this time period (red dashed vertical lines in Fig. 3) and find the time at which the local slope is maximum, which we determine to be the 'mean ENA response time' (red solid vertical lines in Fig. 3). If there are multiple, large peaks in the local slope within this range, then the middle of the range is chosen as the mean response time (for example, Fig. 3a). We note that our choice of using a cubic spline to fit to IBEX data is arbitrary, and points of local maxima or minima in the slopes may shift if a different functional form was used. A higher temporal cadence of measurements from IMAP 57 may be necessary to better constrain the appropriate fitting function. The 'initial ENA response' time t HTS is determined to be the point of local minimum in ENA flux before the mean ENA response time, which we argue is an indication of the time at which the SW pressure increase had reached the HTS and began propagating through the IHS. However, there is a degree of uncertainty as to whether the time at which the local minimum in ENA flux is truly the location of the HTS. The first reason for this is the suggestion from simulations that, as the SW pressure increase reaches the HTS, the HTS first begins to move away from the Sun as the plasma with increased pressure begins propagating through the IHS. As shown in Fig. 3 of McComas et al. (ref. 24 ), a simulation of the response of ENAs from the IHS to the SW dynamic pressure increase first resulted in a slight decrease in ENA flux before a rise in ENA flux began. This decrease appears to be a response to the outwards motion of the HTS due to the increase in SW pressure, which initially decreases the LOS-integrated ENA flux. The outwards motion of the HTS before the rise of ENA flux at 1 au is observed represents a potential uncertainty in the location of the HTS. The second reason for potential uncertainty is the existence of strong fluctuations in ENA flux observed before the rapid rise, which is evident in some pixels near the nose of the heliosphere (Fig. 3a,d,e,g). These fluctuations may adversely affect our ability in finding the 'true' minimum in ENA flux before the rapid rise occurs. A description for how we include these uncertainties in our analysis is given in the section Calculation and propagation of uncertainties.
After performing our analysis, a culling is applied to the results using several criteria. Results are removed if (1) the final mean ENA response time is determined to be after 2019 where we do not have enough data to confidently determine that the ENA fluxes have stopped increasing (5, 1 and 2 pixels for ESA 4, 5 and 6); (2) there are three or more peaks in the slope with heights >50% of the peak slope, making it difficult to identify the actual mean ENA response time (4, 11 and 30 pixels for ESA 4, 5 and 6); (3) the propagated uncertainty of the mean ENA response time is >1 year (76, 26 and 8 pixels for ESA 4, 5 and 6); and (4) finally, some pixels are manually removed due to complexities in the observations making it difficult to determine the ENA response times, for example, multiple peaks are visible similar to criterion no. 2, or there is no clear step function-like rise of the ENA flux (119, 111 and 145 pixels for ESA 4, 5 and 6). After final culling, 673, 728 and 692 pixels in the sky remains for ESA 4, 5 and 6, respectively.

Calculation of distances to the HTS, mean ENA source and HP
After deriving the initial ENA response time (t HTS ) and mean ENA response time (t ENA ) for each accepted pixel in the sky and each ESA 4-6, we calculate the distances to the HTS, mean ENA source region and HP for each pixel. First, the distance to the HTS, r HTS , is calculated by integrating the time for SW propagation from r 0 = 1 au to distance r, plus the time for ENA propagation from r back to r 0 , until it yields the observed initial ENA response time t HTS , such that where u SW is the SW speed and v ENA is the ENA speed. The SW speed is solved as a function of distance from the Sun using spherically symmetric, steady-state fluid transport equations for mass, momentum, magnetic field (B) and pressure (p) of the SW proton ('SWH + '), alpha ('SWHe ++ '), H + PUI ('PUIH + ') and He + PUI ('PUIHe + ') mixture with photoionization and charge exchange source terms, given as 43-45,60-62 1 r d dr (rBu) = 0, (6) where subscript i represents different ion species (that is, mass and pressure terms in equations (3) and (5)). The mass source terms are given as where n represents different neutral species. The momentum source terms are where the first summation over ions i is for charge exchange between SW protons (i = 1) or H + PUIs (i = 2) and neutral H, and the second summation over neutrals n is for photoionization of interstellar H (n = 1) and He (n = 2). The pressure source terms are  Since the probability of neutralization of SW alphas, which is dominated by double charge exchange with interstellar He, results in <1% loss in mass over roughly 100 au (ref. 63 ), their contribution to the momentum and pressure of the plasma mixture is negligible and thus their source terms are ignored.
The relative speeds for charge exchange in the mass, momentum and pressure source terms are given as 61 and the charge exchange cross section for H-H + is 64 Note that the gain/loss of He + by charge exchange for He-He + and H-He + are substantially smaller than photoionization of He and can be ignored.
The total mass density is ρ = m H n SWH + + m He n SWHe ++ + m H n PUIH + + m He n PUIHe + , where we assume all ions co-move at bulk SW speed u, γ = 5/3 is the adiabatic index, p i is the thermal pressure of each ion species, and ν H = 1.44 × 10 −7 s −1 and ν He = 1.14 × 10 −7 s −1 are the neutral H and He photoionization rates at 1 au, respectively, during Carrington rotations (CRs) 2,154-2,156 (ref. 65 ) that are varied with latitude following Bzowski et al. (ref. 66 ). Initial conditions for the SW near the ecliptic are derived from OMNI in-ecliptic observations averaged over the time of the SW pressure increase (CR 2,154-2,156). PUI H + and He + densities are initially zero at r 0 . The SW speed and density at higher latitudes θ are extracted from IPS observations during CR 2,154-2,156 (ref. 42 ). IPS-based SW speeds are first derived from electron density fluctuations along lines of sight near the Sun by defining a power law relationship between those density fluctuations Δn and SW speed u, such that Δn ∝ u a . The power law slope a is approximated by comparing with in-ecliptic SW measurements from the OMNI database and Ulysses observations at high latitudes. Between 1985 and 2008, a value of a = −0.5 was found to derive speeds that best matched OMNI and Ulysses measurements. After 2008, however, a larger, positive slope value of a = 1.0 was required. Tokumaru et al. (ref. 42 ) concluded that the reason for this difference is probably changes in the relationship between the density fluctuations and SW speed with different solar cycles (see their study for more details). However, while their derived SW speeds matched better to OMNI using a power law slope of a = 1.0, it was still overestimating OMNI SW measurements, particularly in 2014. Therefore, we shift the published IPS-derived SW speeds down by 85, 61 and 70 km s −1 for CR 2,154, 2,155 and 2,156, respectively, to better match OMNI. Considering that the reason for this shift is not well understood, we include a 15% relative uncertainty of the SW speeds in our analysis. After shifting the SW speeds, the plasma density as a function of latitude is calculated assuming constant dynamic pressure with latitude, that is, [ρu 2 ] θ = [ρu 2 ] θ=0 , based on analyses of Ulysses observations 67 . We note that the assumption of latitudinal invariance of SW dynamic pressure does not significantly affect our results. The most important factor that affects the timing for SW propagation to the HTS is the SW speed measured at 1 au, while the SW density acts as a less effective, higher-order factor in determining the mass-loading of the SW from 1 au to the HTS. We include an uncertainty for SW density in our analysis, as described further below, but it does not contribute significantly to the uncertainties of the distance results.
We note that IPS-derived SW speeds show an abrupt increase in the southern hemisphere in CR 2,156 compared to CR 2,154 and 2,155. This behaviour appears to be caused by an emission of fast SW from a large coronal hole at mid-latitudes in the southern hemisphere in late 2014 as seen in the Solar Dynamics Observatory (SDO) Atmospheric Imaging Assembly (AIA)/Helioseismic and Magnetic Imager (HMI) observations (Fig. 1f) (https://sdo.gsfc.nasa.gov/data/aiahmi/). This coronal hole is persistent over multiple CRs to early 2015, indicating that the fast SW speeds in the southern hemisphere in CR 2,156 are important to include in our analysis. Because fast SW streams will interact with slow streams preceding them in time, and since our model cannot simulate the fast-slow SW speed interaction, we give larger weighting to SW speeds in CR 2,156 when calculating the weighted average in Fig. 1e (25% for CR 2,154, 25% for CR 2,155 and 50% for CR 2,156). The weighted standard deviation of the average, shown in grey in Fig. 1e, indicates that a relative uncertainty of 15% applied to SW speeds at all latitudes is sufficient to capture the potential uncertainties in our model. If we were to use SW speeds from CR 2,156 only, our HTS and HP distances would move slightly outwards at roughly 0 to −45° latitudes, but the time at which we begin the SW propagation at 1 au would be roughly 0.05 years after that used in our analysis. Ultimately, this would move our heliosphere boundaries only by a few au, and thus not enough to make a statistically significant difference.
The temperature of SW protons, if just solved using equations (2)-(6), yield values below 1,000 K at the HTS. However, Voyager 2 observations clearly show that the SW proton temperature does not decrease adiabatically with distance from the Sun and slightly increases with distance beyond roughly 20 au from the Sun 68 . The reason for this non-adiabatic heating has been studied in detail in the past and is probably due to turbulent heating by waves excited by interstellar PUI injection 45,[69][70][71][72][73][74] . While it is beyond the scope of our analysis to include a turbulent heating source term for SW protons in equations (2)-(6), we can put a lower limit on the SW proton temperature that is roughly consistent with Voyager observations. Thus, when solving the transport of SW proton pressure, we force their temperature to always be ≥10 4 K. This assumption does not significantly affect our results, however, since the interstellar H and He PUIs dominate the internal pressure of SW in the outer heliosphere. We note that New Horizons' SWAP observations show the H + PUI temperature is roughly 4 × 10 6 K at 30 au from the Sun in late 2014 (ref. 75 ), which is close to our model prediction of 4.0 × 10 6 K at the same distance. This does not necessarily suggest our model is consistent with SWAP at other times or distances from the Sun, because SWAP observations show PUIs experience non-adiabatic heating from a physical process that is not yet fully understood.
The total thermal pressure p = p e + 4 ∑ i=1 p i includes the pressure of electrons and all ion components. We assume quasi-neutrality is maintained throughout the system. The temperature of electrons in the outer heliosphere is not well known, but there is reason to believe they contain non-negligible suprathermal distributions. Electrons may be substantially heated at interplanetary shocks, maintaining high internal energies compared to the thermal SW protons [76][77][78] . Therefore, we assume that electron temperatures are ten times higher than the SW protons. The interstellar neutral H density, n H , is extracted from a global, three-dimensional steady-state simulation of the heliosphere based on the methodology in Zirnstein et al. (ref. 79 ). The simulation boundary conditions at 1 au are similar to the previous work, but the interstellar neutral H density was increased such that the H density near the upwind HTS is consistent with recent measurements from New Horizons' SWAP 80 . Interstellar neutral H distribution is assumed to be Maxwellian moving at a bulk speed of u H = 22 km s −1 and their inflow direction is (252°.2, 9°) 81,82 .
The interstellar neutral He density, n He (r, θ), is calculated analytically for a cold gas 8,83 , such that Article https://doi.org/10.1038/s41550-022-01798-6 where n He,∞ = 0.015 cm −3 is the interstellar neutral He density far from the Sun and λ He = 0.5 au is the size of the He density depletion region due to ionization 84,85 , G is the gravitational constant, M is the solar mass, u He,∞ = 25.4 km s −1 is the interstellar neutral He speed with inflow direction (255°.7, 5°.1) 86 , μ = 0 is the gravity compensation factor due to solar radiation pressure, θ is the angle of vector r from the neutral He upwind direction and θ j = θ if p j > 0 and θ j = 2π − θ if p j < 0 (see ref. 83 for more details). By solving equations (2)-(6), the bulk SW speed u (r) is calculated for each pixel direction in the sky as a function of distance r from the Sun in equation (1), therefore allowing us to derive r HTS . Next, the distance from the HTS to the mean ENA source and HP is calculated. First, we estimate the HTS compression ratio by solving the shock adiabatic equation for a perpendicular shock, given as where R is the (unique) shock compression ratio, β u is the upstream plasma beta and M u is the upstream plasma Mach number. We calculate the upstream SW bulk flow speed, u u , effective thermal pressure, p u , for all electron+ion components and effective mass density, ρ u , and magnetic field B u , from the solution of equations (2)- (6). We note that the effective specific heat ratio γ pressure term γp u need not assume that the index γ = 5/3 for all ion species, since it is possible that interstellar PUIs behave non-adiabatically due to their unique occupation of phase space and behaviour in the SW 75 . To allow for this possibility, we assumed that γ SW = γ PUI = γ = 5/3 for all particles but introduce a relative uncertainty that accounts for the possibility that γ might range between 1.33 and 2.0 (Calculation and propagation of uncertainties section). Thus, the effective specific ratio upstream of the HTS is The final step before calculating the mean ENA source and HP distances involves calculating the fast magnetosonic wave speed in the IHS, which we assume is the dominant wave speed in the IHS. The effective pressure term downstream of the HTS, p d , is readily calculated from the Rankine-Hugoniot jump conditions for a perpendicular shock as a function of the upstream plasma properties 87 . We also include a contribution of pressure from anomalous cosmic rays (ACRs) that may be on the order of 30% of the total pressure 18,88,89 . Thus, the total effective plasma pressure downstream of the HTS is modified to be p d,tot = p d / (1 − f ACR ), where f ACR = 0.3. The fast magnetosonic wave speed is then calculated as where B d = B u R. The IHS plasma flow speed immediately downstream of the HTS is derived using the shock compression ratio from equation (13), such that u d,0 = u u /R, as well as the fast magnetosonic wave speed using the downstream plasma pressure. The downstream flow speed and wave speed are used to simultaneously calculate the radial distance from the HTS through the IHS at which the mean ENA source and HP are located. This is done by performing an iterative, binary search for the optimal solution for the position of the HP that, using the previously derived IHS flow and wave speeds, yields the correct time delay observed by IBEX. The search behaves as: (1) Define an initial search range of Δr i HP,min < Δr i HP < Δr i HP,max (where i represents the step iteration), assuming Δr i HP,min = 2 au and Δr i HP,max = 70 au. The distance to the HP from the HTS is assumed to be in the middle of the range, that is, Δr i HP = (Δr i HP,min + Δr i HP,max ) /2.
(2) Calculate time Δt 1 it takes for the forward-propagating wave at speed u fw (r) = u w + u a (r) to travel from r HTS to r HP , where u w is the fast magnetosonic wave speed and u a (r) is the advecting flow speed (see details below).  Steps 2-6 are iteratively repeated until the optimal choice for Δr HP (and thus Δr ENA ) is found with an accuracy of <0.5 au.
We note that this process computes a radial distance from the HTS with an estimation for IHS plasma flow deflection away from the radial vector. There are no direct observations of IHS plasma flow deflection except for measurements from the Voyager spacecraft over two directions in the sky, which may not be applicable over all directions in our analysis. From global, steady-state simulations, we expect that the plasma flow is slowed near the IHS stagnation point and deflected away from it, although the existence and location of a stagnation point depends on asymmetries induced by the interstellar magnetic field, time-dependent solar cycle effects and corotating interaction regions, and instabilities developing near the HP 32,90-95 .
We first approximate the amount of flow deflection with help from the global heliosphere simulation 79 . We calculate the average flow deflection angle in the IHS plasma as a function of direction in the sky from the simulation, weighted by the 4.3 keV ENA source in the IHS (see ref. 79 , section 2.2). From the simulation, we find minimal deflection (roughly 0°) near the simulated stagnation point located near ecliptic (267°, −4°), and the deflection angle increases to a maximum of roughly 45° at an angle of roughly 40° away from the stagnation point, nearly symmetric in longitude and latitude. However, the true IHS stagnation point is probably roughly 30° below the nose as determined from IBEX and Voyager observations 47 . Therefore, we use the information from the simulation but modify it to better match these and other observations. We define a function such that the flow deflection angle is zero at (255°, −27°), increases proportional to √φ (where φ is the angular separation of the pixel from the stagnation point) and maximizes as 45° at φ = 40°. This indicates a decrease of the radial plasma flow speed u a (r) by a factor of cos 45 ∘ = 0.7. Next, we incorporate information from Voyager observations. While Voyager 1 observations indicate slowing Article https://doi.org/10.1038/s41550-022-01798-6 can be as large as roughly 50% halfway through the IHS 96 , Voyager 2 observations show less slowing (roughly 25%) 97 and recent analyses suggest radial plasma flow velocities derived from Voyager 1 energetic particle measurements may be inaccurate 97,98 . Moreover, these observations at Voyager 1 and 2 are probably coupled to time-dependent, solar cycle effects that are nearly impossible to predict for our analysis 95 . Thus, our analysis can only include a rough approximation of this effect. We approximate the IHS plasma flow speed as a function of distance r from the HTS according to where u d,0 is the initial downstream flow speed, the second term introduces slowing where Γ = 0.5, and the final term requires the flow to reach 0 at the HP for any arbitrary value of Γ. For the nominal value of Γ = 0.5, equation (16) requires that the radial flow speed decreases to 75% halfway through the IHS, similar to Voyager 2 measurements, and drops more quickly closer to the HP. The distances r HTS , Δr ENA and Δr HP are solved as a function of ENA speed, v ENA , and must be integrated over each IBEX ESA energy passband. Because the IBEX-Hi ESA passbands cover a relatively wide range of ENA energies with a full-width at half-maximum of roughly 60% (ref. 33 ), these results must be repeated for a range of ENA energies over ESA passbands 4-6 and weighted by the instrument energy-dependent response functions and ENA flux spectra. Therefore, we solve for r HTS , Δr ENA and Δr HP over a range of ENA speeds and average the results as The weights W HTS and W ENA are calculated as a function of the IBEX ESA energy response function R(v) 46,99 and the observed GDF ENA spectral indices η HTS and η ENA measured by IBEX at times t HTS and t ENA , respectively. Because the observations of t HTS and t ENA are made at different times, the ENA spectral index is different for the HTS, mean ENA source and HP distance results (note that we use the same weight for Δr ENA and Δr HP because the same observation time is used to derive them). We estimate η HTS and η ENA as a function of longitude, latitude and time using results from Swaczyna et al. (ref. 38 ), who performed a spherical harmonic decomposition of the IBEX GDF observations after separating out the ribbon and provided full sky maps of the GDF at all IBEX-Hi energies. We use Compton-Getting and survival-probability corrected GDF results derived by their analysis and compute the ENA spectral indices between ESA 3-5 and ESA 4-6 as a function of longitude, latitude and time to estimate η HTS and η ENA in our results. Spectral indices between ESA 3-5 are used for the distance calculations for ESA 4, and the spectral indices between ESA 4-6 are used for the distance calculations for ESA 5 and 6. The derived spectral indices are interpolated in time at t HTS and t ENA for each pixel in the sky and used in equations (17)- (19).
The results for ⟨r HTS ⟩, ⟨Δr ENA ⟩ and ⟨Δr HP ⟩ are obtained for each ESA 4-6 after integrating equations (17)- (19). Then, we combine the results over energy by averaging the distances with weights determined by the propagated variances. Their corresponding uncertainties are calculated by the propagation of uncertainties of multiple variables used in the analysis (next section).

Calculation and propagation of uncertainties
Our analysis includes multiple sources of uncertainties and propa gates the uncertainties when calculating the distances to the HTS, mean ENA source and HP. The parameters with uncertainties are listed below: (1) IBEX ENA fluxes, J ENA . We propagate the statistical uncertainties of the IBEX ENA fluxes through the analysis. The relative uncertainties are typically a few percent and therefore do not contribute significantly to most of the results.  Fig. 3a,e,g, is not accounted for in the simulation and therefore is used as an estimate of uncertainty in our analysis. This fluctuations before the pressure increase may be realistic effects of the outer heliosphere, or potentially due to imperfect Compton-Getting corrections or background subtractions within ESA 5 and 6 (refs. 33,100 ). Regardless of the origin, we attempt to account for this uncertainty by (1) calculating the standard deviation in IBEX ENA flux over a 1-year period before t HTS , that is, s J , (2) adding s J to the flux at the initial response time, that is J (t HTS ) + s J , (3) finding the point in time after t HTS when the observed flux J = J (t HTS ) + s J = J (t * ) and (4) calculating the difference t * − t HTS . This difference is added in quadrature to the other propagated uncertainties. This uncertainty represents the largest uncertainty in many pixels of the sky. (3) Mean ENA response time, t ENA . The uncertainty of the mean ENA response time is propagated through the calculation of the heliospheric boundary distances. This uncertainty is primarily composed of the propagated ENA flux uncertainties. While variability in ENA flux exists before the large response of ENAs to the SW pressure event may affect the initial ENA response time significantly, the same may occur with the distance to the HP, but changes to the HP occur more slowly over time and it is expected that wave reflection from the HP happens before the boundary moves outwards by any noticeable distance. The uncertainties in the mean ENA response time are small due to the smooth gradient in ENA flux and are largely due to uncertainties in the IBEX data.  80 ) recently updated the interstellar neutral H density within the outer heliosphere, yielding 0.127 ± 0.015 cm −3 at the upwind HTS. This is roughly 40% higher compared to previous work, which is obtained from the first outer heliosphere measurements of interstellar H + PUIs by New Horizons SWAP that provided a better estimation of the parent neutral H density. The uncertainty of n H obtained in that study includes an estimated uncertainty of the charge exchange cross section of the order of 10%. Therefore, to avoid double counting of the uncertainty of the cross section, we only include the combined uncertainty of interstellar neutral H density from the Swaczyna et al. analysis. (7) SW proton temperature 'floor'. We noted that the SW proton temperature solved using equations (2)-(6) in the supersonic SW yield unrealistically low values without including turbulent heating source terms. To account for this, we force the SW proton temperature to be at least 10 4 K at each step in the solution. Clearly, there is significant uncertainty in this approach; therefore, we assume a relative uncertainty of 100% of this parameter and propagate it through the analysis. (8) SW electron temperature, T e . The temperature of electrons in the supersonic SW is assumed to be ten times higher than the SW protons, which is an assumption based largely on extensive theoretical calculations [76][77][78]102 . Therefore, we assume a relative uncertainty of 100% (that is, ranging between 0 and 20 times the SW proton temperature) and propagate it through the analysis.
(9) Specific heat ratio of PUIs, γ. Due to the non-thermal distribution of PUIs and their preferential nature of heating at shocks in the outer heliosphere, it is not clear what the specific heat ratio of PUIs is near the HTS. For simplicity, we assume γ = 5/3. New Horizons' SWAP observations of non-adiabatic PUI heating in the outer heliosphere show that the 'cooling index' of PUIs, α, which is related to the specific heat ratio as α = 1/ (γ − 1), is roughly 2.1 halfway to the HTS and may increase to roughly 2.9 at the HTS 75 . Because this is the only direct evidence of the specific heat capacity of PUIs, we assume a relative uncertainty of 20% for γ, such that within 1-sigma of uncertainty, γ may be between 1.33 and 2.0 (or α varies between 3 and 1, respectively). (10) HTS compression ratio, R. The kinetic nature of particle heating and acceleration at the HTS probably means that our use of the single-fluid, ideal shock adiabatic equation to derive the HTS compression has some level of uncertainty. While the compression ratio derived using equation (13) yields values that appear consistent with measurements from Voyager 2 (ref. 103 ) and predictions from particle-in-cell simulations 104 , we include a 1-sigma relative uncertainty of 10% for the HTS compression ratio in our analysis. (11) IHS plasma flow speed, u d . Downstream of the HTS, the plasma flows through the IHS and is deflected away from the radial vector and slowed by compression or deflection near the HP. Considering the substantial differences in Voyager 1 and 2 observations (or differences in interpretations of the data), and what little is known about the global IHS plasma flow, we introduce an uncertainty of the flow slowing factor Γ = 0.5 in an amount of σ Γ = 1/√2 − Γ, such that the radial plasma speed halfway through the IHS may be between 75 and 50% slower than its initial speed at the HTS. (12) ACR pressure contribution in the IHS, p ACR . In our analysis, we assume that ACRs contain 30% of the total effective thermal pressure of the IHS plasma, based on a recent analysis of Voyager and IBEX observations 18 . We assume a relative uncertainty of 33% of this parameter and propagate it through the analysis. (13) ENA source region thickness, Δl ENA . To determine the optimal position of the HP, we must assume a distance between the backwards-propagating reflected wave and forwards-propagating advecting flow that coincides with the mean ENA response time.
Based on simulation results of Zirnstein et al. (ref. 26 ), the half-width of the 4.3 keV ENA source thickness Δl ENA is approximately 25% of the distance from the HTS to the HP over most directions of the sky used in our analysis. Because we assume that the overlap between the wave and advecting flow must be approximately Δl ENA /2 to coincide with the mean ENA response time, t ENA , our calculation of ⟨Δr ENA ⟩ and ⟨Δr HP ⟩ relies strongly on this assumption. Therefore, we assume an uncertainty of 100% for this parameter, such that the overlap region might be anywhere between 0 and twice the half-width of the ENA source region and that it cannot be greater than half of the total IHS thickness.
The uncertainties listed above are propagated through each step of the analysis. This is performed by manually varying the values of each parameter by its 1-sigma uncertainty, recalculating the desired variable with the perturbed parameter and adding the deviations of the results in quadrature to estimate the final propagated uncertainty. For example, when calculating ⟨Δr ENA ⟩, its uncertainty is calculated as where index j represents the parameter whose value was increased by its 1-sigma uncertainty before recalculating ⟨Δr ENA ⟩| j . We note that this method assumes all parameters are independent.
Article https://doi.org/10.1038/s41550-022-01798-6 Averaging results with uncertainties, such as angular smoothing over multiple pixels, is performed by weighting values by their inverse variances and calculating the uncertainty of the average. As an example, for arbitrary variable g with uncertainty σ g , its weighted average, ⟨g⟩, is calculated as The uncertainty, σ ⟨g⟩ , is calculated from one of the following equations: the propagated uncertainty or the statistical uncertainty where N eff is the effective number of measurements. In the case where the uncertainties of all points are similar, N eff approaches the actual number of points, N. Equation (21) is used to calculate the weighted average of any variables throughout the analysis. The uncertainty of the weighted average is chosen from either equation (22) or equation (23), whichever is larger.

Data availability
The results reported in the study shown in Fig. 4 and Fig. 5a-f are publicly available to download. Source data are provided with this paper.