Oblique and rippled heliosphere structures from the Interstellar Boundary Explorer – Nature Astronomy

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. ¹⁶). 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⁴⁶. 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²⁴. 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⁵⁷ 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. ²⁴), 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₀ = 1 au to distance r, plus the time for ENA propagation from r back to r₀, 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,44,45,60,61,62}

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. ⁶³), 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⁶¹

and the charge exchange cross section for H–H⁺ is⁶⁴

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 (rho = m_{mathrm{H}}n_{{mathrm{SWH}}^ + } + m_{{mathrm{He}}}n_{{mathrm{SWHe}}^{ + + }} +) (m_{mathrm{H}}n_{{mathrm{PUIH}}^ + } + m_{{mathrm{He}}}n_{{mathrm{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 (nu _{mathrm{H}} = 1.44 times 10^{ – 7},{mathrm{s}}^{ – 1}) and (nu _{{mathrm{He}}} = 1.14 times 10^{ – 7},{mathrm{s}}^{ – 1}) are the neutral H and He photoionization rates at 1 au, respectively, during Carrington rotations (CRs) 2,154–2,156 (ref. ⁶⁵) that are varied with latitude following Bzowski et al. (ref. ⁶⁶). 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₀. The SW speed and density at higher latitudes θ are extracted from IPS observations during CR 2,154–2,156 (ref. ⁴²). 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. ⁴²) 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⁻¹ 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, (left[ {rho u^2} right]_theta = left[ {rho u^2} right]_{theta = 0}), based on analyses of Ulysses observations⁶⁷. 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⁶⁸. 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⁴ 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⁶ K at 30 au from the Sun in late 2014 (ref. ⁷⁵), which is close to our model prediction of 4.0 × 10⁶ 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 + mathop {sum }limits_{i = 1}^4 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. ⁷⁹). 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⁸⁰. Interstellar neutral H distribution is assumed to be Maxwellian moving at a bulk speed of u_H = 22 km s⁻¹ 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

where (n_{{mathrm{He}},infty } = 0.015) cm⁻³ is the interstellar neutral He density far from the Sun and (lambda _{{mathrm{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_{{mathrm{He}},infty } = 25.4) km s⁻¹ is the interstellar neutral He speed with inflow direction (255°.7, 5°.1)⁸⁶, μ = 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. ⁸³ 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⁷⁵. To allow for this possibility, we assumed that (gamma _{{mathrm{SW}}} = gamma _{{mathrm{PUI}}} = gamma = 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⁸⁷. 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_{{mathrm{d}},{mathrm{tot}}} = p_{mathrm{d}}/left( {1 – f_{{mathrm{ACR}}}} right)), where f_ACR = 0.3. The fast magnetosonic wave speed is then calculated as

where B_d = B_uR.

The IHS plasma flow speed immediately downstream of the HTS is derived using the shock compression ratio from equation (13), such that (u_{{mathrm{d}},0} = u_{mathrm{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. (1)

   Define an initial search range of ({Delta}r_{{mathrm{HP}},{mathrm{min}}}^i < {Delta}r_{{mathrm{HP}}}^i < {Delta}r_{{mathrm{HP}},{mathrm{max}}}^i) (where i represents the step iteration), assuming ({Delta}r_{{mathrm{HP}},{mathrm{min}}}^i = 2) au and ({Delta}r_{{mathrm{HP}},{mathrm{max}}}^i = 70) au. The distance to the HP from the HTS is assumed to be in the middle of the range, that is, ({Delta}r_{{mathrm{HP}}}^i = left( {{Delta}r_{{mathrm{HP}},{mathrm{min}}}^i + {Delta}r_{{mathrm{HP}},{mathrm{max}}}^i} right)/2).

2. (2)

   Calculate time Δt₁ it takes for the forward-propagating wave at speed (u_{{mathrm{fw}}}left( r right) = u_{mathrm{w}} + u_{mathrm{a}}left( r right)) 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).

3. (3)

   Propagate the advecting flow at speed u_a(r) for the same amount of time Δt₁ through the IHS.

4. (4)

   Propagate the advecting flow and backward-propagating reflected wave (moving at speed (u_{{mathrm{rw}}}left( r right) = u_{mathrm{w}} – u_{mathrm{a}}left( r right))) simultaneously over time until they cross each other and their difference in position ({Delta}r_{{mathrm{af}}} – {Delta}r_{{mathrm{rw}}} = {Delta}l_{{mathrm{ENA}}}/2) (where we have assumed that ({Delta}l_{{mathrm{ENA}}}/2 cong left( {r_{{mathrm{HP}}} – r_{{mathrm{HTS}}}} right)/4) based on simulations²⁶), yielding elapsed time Δt₂. Note that Δr_af and Δr_rw are radial distances from the HTS.

5. (5)

   If the HP position r_HP is the optimal choice, the total elapsed time should equal ({Delta}t_1 + {Delta}t_2 = t_{{mathrm{ENA}}} – t_{{mathrm{HTS}}} – {Delta}t_{{mathrm{ENA}} to {mathrm{HTS}}}), where ({Delta}t_{{mathrm{ENA}} to {mathrm{HTS}}}) is the amount of time it takes ENAs at speed v_ENA to travel from the r_ENA to r_HTS. The mean ENA source distance from the HTS is then defined as ({Delta}r_{{mathrm{ENA}}} = left( {{Delta}r_{{mathrm{af}}} + {Delta}r_{{mathrm{rw}}}} right)/2), that is, the middle of the intersection of the advecting flow and reflected wave.

6. (6)

   Compare total elapsed time from Step 5 (({Delta}t_1 + {Delta}t_2)) to the time based on IBEX observations ((t_{{mathrm{ENA}}} – t_{{mathrm{HTS}}} – {Delta}t_{{mathrm{ENA}} to {mathrm{HTS}}})). If the elapsed time from Step 5 is less than the observed time, update the search range for ({Delta}r_{{mathrm{HP}}}^{i + 1}) to be (left( {{Delta}r_{{mathrm{HP}},{mathrm{min}}}^{i + 1},{Delta}r_{{mathrm{HP}},{mathrm{max}}}^{i + 1}} right) = ({Delta}r_{{mathrm{HP}}}^i,{Delta}r_{{mathrm{HP}},{mathrm{max}}}^i)). If the elapsed time from Step 5 is greater than the observed time, update the search range for ({Delta}r_{{mathrm{HP}}}^{i + 1}) to be (left( {{Delta}r_{{mathrm{HP}},{mathrm{min}}}^{i + 1},{Delta}r_{{mathrm{HP}},{mathrm{max}}}^{i + 1}} right) = ({Delta}r_{{mathrm{HP}},{mathrm{min}}}^i,{Delta}r_{{mathrm{HP}}}^i)). The new estimate for the distance of the HP from the HTS is ({Delta}r_{{mathrm{HP}}}^{i + 1} = left( {{Delta}r_{{mathrm{HP}},{mathrm{min}}}^{i + 1} + {Delta}r_{{mathrm{HP}},{mathrm{max}}}^{i + 1}} right)/2).

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,91,92,93,94,95}.

We first approximate the amount of flow deflection with help from the global heliosphere simulation⁷⁹. 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. ⁷⁹, 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⁴⁷. 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 (sqrt varphi) (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^circ = 0.7). Next, we incorporate information from Voyager observations. While Voyager 1 observations indicate slowing can be as large as roughly 50% halfway through the IHS⁹⁶, Voyager 2 observations show less slowing (roughly 25%)⁹⁷ 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⁹⁵. 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. ³³), 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. ³⁸), 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 (leftlangle {r_{{mathrm{HTS}}}} rightrangle), (leftlangle {{Delta}r_{{mathrm{ENA}}}} rightrangle) and (leftlangle {{Delta}r_{{mathrm{HP}}}} rightrangle) 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 propagates the uncertainties when calculating the distances to the HTS, mean ENA source and HP. The parameters with uncertainties are listed below:

1. (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.

2. (2)

   Initial ENA response time, t_HTS. The uncertainty of the initial ENA response time is propagated through the calculation of the heliospheric boundary distances. This uncertainty is partly composed of the propagated ENA flux uncertainties. Also, the location of the HTS based on the point of local minimum in ENA flux right before the sharp rise is not well constrained due to the potentially substantial fluctuation of ENA fluxes over time unrelated to the large SW pressure increase in late 2014. This variability, examples of which are visible in 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 (Jleft( {t_{{mathrm{HTS}}}} right) + s_J), (3) finding the point in time after t_HTS when the observed flux (J = Jleft( {t_{{mathrm{HTS}}}} right) + s_J = Jleft( {t^ ast } right)) and (4) calculating the difference (t^ ast – t_{{mathrm{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. (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.

4. (4)

   SW pressure increase start time. The start time of the SW pressure increase at 1 au is assumed to be near the time of peak gradient in the SW pressure. Specifically, we choose to average the median times of the three CRs available in IPS observations surrounding the peak gradient in SW pressure, that is, CR 2,154–2,156, with double weighting for CR 2,156 due to the abrupt rise in SW speed from a coronal hole in the southern hemisphere. The weight-averaged start time is calculated to be 2014.78. We include a 1-sigma uncertainty of one CR (27.3 days or roughly 0.075 years) in the start time and add it in quadrature to the initial and mean ENA response time uncertainties before calculating the heliosphere boundary distances.

5. (5)

   IPS SW speed, u₀, and density, n₀. The absolute accuracy of SW speeds derived from IPS observations is difficult to quantify at high latitudes, especially without in situ observations from Ulysses that were available before 2009. Sokół et al. (ref. ¹⁰¹) and Tokumaru et al. (ref. ⁴²) estimated the uncertainty in IPS-derived SW speeds by calculating the relative variance between IPS and Ulysses observations during polar orbits, yielding a relative error of roughly 10% for speed and roughly 20–30% for density. Our analysis uses SW speeds derived from IPS observations during CR 2,154–2,156, which approximately overlap the middle of the rise in SW dynamic pressure in late 2014. While this choice appeared sufficient for our purposes, and we attempted to account for all potential uncertainties in our analysis, we cannot definitively claim that averaging over these three CRs is the best choice of SW speed to use in our analysis. Second, IPS-derived SW speeds at low latitudes after 2009 are systematically higher by roughly 50–100 km s⁻¹ than speeds in the OMNI database observed in situ by spacecraft, the reason for which is not clear but appears to be caused by the evolving relationship between the SW speed and density fluctuations in different solar cycles^42,65. Because there was agreement between IPS and both OMNI and Ulysses measurements before 2009 (ref. ¹⁰¹), it is assumed that any difference between IPS and OMNI now is systematic over all latitudes and thus the IPS-derived SW speeds can be shifted by a common value at all latitudes to match OMNI data⁴². However, it is not clear how accurate this approach is without in situ measurements at high latitudes. Finally, it is important to note that IPS-derived SW speeds in CR 2,156 are higher at heliolatitudes roughly 0 to −45° compared to other latitudes, which approximately overlaps the region in the sky where the HTS and HP distances derived in our analysis are closest to the Sun. 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 year 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. On the basis of this analysis, and the uncertainty associated with needing to shift IPS-derived SW speeds to match OMNI, in our analysis we assume a relative error of 15% for IPS-derived SW speeds and 25% for IPS-derived SW density and propagate their uncertainties through the analysis.

6. (6)

   Interstellar neutral H density, n_H, and charge exchange cross section, σ. Swaczyna et al. (ref. ⁸⁰) recently updated the interstellar neutral H density within the outer heliosphere, yielding 0.127 ± 0.015 cm⁻³ 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. (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⁴ 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. (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. (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 (alpha = 1/left( {gamma – 1} right)), is roughly 2.1 halfway to the HTS and may increase to roughly 2.9 at the HTS⁷⁵. 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. (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. ¹⁰³) and predictions from particle-in-cell simulations¹⁰⁴, we include a 1-sigma relative uncertainty of 10% for the HTS compression ratio in our analysis.

11. (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 (sigma _{{{mathrm{{Gamma}}}}} = 1/sqrt 2 – {{{mathrm{{Gamma}}}}}), 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. (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¹⁸. We assume a relative uncertainty of 33% of this parameter and propagate it through the analysis.

13. (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. ²⁶), 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 (leftlangle {{Delta}r_{{mathrm{ENA}}}} rightrangle) and (leftlangle {{Delta}r_{{mathrm{HP}}}} rightrangle) 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 (leftlangle {{Delta}r_{{mathrm{ENA}}}} rightrangle), its uncertainty is calculated as

where index j represents the parameter whose value was increased by its 1-sigma uncertainty before recalculating (left. {leftlangle {{Delta}r_{{mathrm{ENA}}}} rightrangle } right|_j). We note that this method assumes all parameters are independent.

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, (leftlangle g rightrangle), is calculated as

The uncertainty, (sigma_{leftlangle g rightrangle}), 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.