IBEX Ribbon Separation Using Spherical Harmonic Decomposition of the Globally Distributed Flux

Remote imaging of plasmas in the heliosphere and very local interstellar medium is possible with energetic neutral atoms (ENAs), created through the charge exchange of protons with interstellar neutral atoms. ENA observations collected by the Interstellar Boundary Explorer (IBEX) revealed two distinctive sources. One source is the globally distributed flux (GDF), which extends over the entire sky and varies over large spatial scales. The other source encompasses only a narrow circular band in the sky and is called the IBEX ribbon. Here, we utilize the observed difference in spatial scales of these two ENA sources to separate them. We find that linear combinations of spherical harmonics up to degree ℓmax=3 can reproduce most of the ENA fluxes observed outside the ribbon region. We use these combinations to model the GDF and the difference between the observed fluxes and the GDF yields estimation of the ribbon emission. The separated ribbon responds with a longer time delay to the solar wind changes than the GDF, suggesting a more distant source of the ribbon ENAs. Moreover, we locate the direction of the maximum plasma pressure based on the GDF. This direction is 17°.2 ± 0°.5 away from the upwind direction within the plane containing the interstellar flow and interstellar magnetic field vectors. This deflection is consistent with the expected position of the maximum external pressure at the heliopause. The maps with separated ribbon and GDF are posted concurrently with this paper and can be used to further study these two sources.


Introduction
Observations of energetic neutral atoms (ENAs), created through a charge exchange of suprathermal ions with interstellar neutral (ISN) atoms, enable studies of remote heliospheric plasmas (Gruntman 1997). The Interstellar Boundary Explorer (IBEX; McComas et al. 2009a) is an Earth-orbiting spacecraft dedicated to measurements of ENA fluxes. IBEX is equipped with two ENA imagers: IBEX-Lo covering the energy range from 10 eV to 2 keV (Fuselier et al. 2009b) and IBEX-Hi observing ENAs in the range 0.5-6 keV ). While IBEX-Lo focuses on observations of low-energy ISN atoms (Möbius et al. 2009) originating from the very local interstellar medium (VLISM) and inflowing into the heliosphere, ENAs measured by IBEX-Hi originate from neutralized protons in the outer heliosphere.
The first IBEX observations provided ENA sky maps showing an arc-like structure of enhanced ENA fluxes called the IBEX ribbon (Fuselier et al. 2009a;McComas et al. 2009b;Schwadron et al. 2009). The ribbon is a relatively narrow feature (∼20°wide), which extends over a nearly complete circle in the sky with a half cone angle of ∼75° (Funsten et al. , 2015. This structure had not been anticipated before IBEX was launched, and thus this discovery led to more than a dozen hypotheses explaining this phenomenon (see review by McComas et al. 2014b), which placed the source of the ribbon across various regions of the heliosphere or even at a nearby boundary layer in the VLISM. Several studies of the ribbon properties support the so-called secondary ENA mechanism, in which the neutral solar wind (primary ENAs) escapes the heliosphere, where it is ionized and subsequently re-neutralized forming secondary ENA fluxes (Heerikhuisen et al. 2010). The secondary ENAs travel back to the inner heliosphere, where they are measured, only if the magnetic field direction in the source region is almost perpendicular to the line of sight as viewed by IBEX (Schwadron et al. 2009). This process can be explained if pickup ions from the neutral solar wind are almost scatter-free and thus form a stable ring distribution where their velocities are perpendicular to the magnetic field lines in the VLISM (Chalov et al. 2010;Möbius et al. 2013;Zirnstein et al. 2018a), or if these pickup ions are spatially retained and create a higher density within the source region Isenberg 2014).
The ENA fluxes measured in any direction of the sky may contain more than one emission source. Thus, quantitative study of the ribbon as a separate source region must include identifying and separating multiple superposed emission sources. For example, we have strong evidence of a ubiquitous, slowly varying ENA emission (the globally distributed flux, or GDF) from the heliosheath, which is present in all directions in the sky ). Thus, identifying and quantifying the spatial, spectral, and temporal properties of the ribbon requires careful disentanglement of the GDF contribution to the measured ENA emission. In order to separate ribbon flux from other emission sources, Funsten et al. (2013Funsten et al. ( , 2015, Swaczyna et al. (2016a), Dayeh et al. (2019), and Reisenfeld et al. (2021) assumed various functional forms of the ribbon profile, which allowed for separations of the ribbon flux from the GDF. However, since the ribbon profile is different for various considered models , the ribbon estimated from these studies may be biased by assuming specific functional forms of this profile. Furthermore, since the ribbon location in the sky is influenced by the slow-fast latitudinal structure of the solar wind (Swaczyna et al. 2016b) and the draping of the interstellar magnetic field around the heliopause (Zirnstein et al. 2016a;Dayeh et al. 2019), a model-independent separation of the ribbon from the underlying ENA emission sources is also crucial to determine the direction of the interstellar magnetic field (Zirnstein et al. 2016b). Schwadron et al. (2011Schwadron et al. ( , 2014Schwadron et al. ( , 2018) pursued a fundamentally different approach to ribbon separation in which the GDF is interpolated across the ribbon region based on the surrounding fluxes using a system of weighting factors, which were found through an iterative process. This procedure estimates the GDF inside the ribbon based on the GDF outside the ribbon. However, as an initial guess of the ribbon, the method uses the Gaussian fit to the ribbon profile and is, therefore, susceptible to a model-dependent bias, as discussed above. Furthermore, the technique was only applied to time-combined maps and depended on a selection of weights.
In this study, we describe a new separation method based on the spherical harmonic decomposition of the GDF outside the ribbon-covered region and its interpolation within this region. This approach eliminated the model-dependent biases in the ribbon shape that influence separation results. Spherical harmonics provide an orthonormal basis of functions on a sphere (Arfken et al. 2012). Consequently, any function defined on a sphere can be approximated with any desired precision as a sum of spherical harmonics, and the function can be characterized by coefficients of this sum. Such data representation has been used, e.g., for observations of the cosmic microwave background (e.g., Kamionkowski et al. 1997) and Earth gravity models (e.g., Pavlis et al. 2012). Furthermore, spherical harmonic representations are also helpful in connection with more universal pixelization schemes such as HEALPix (Górski et al. 2005), which may be used in the future for ENA observations from the upcoming Interstellar Mapping and Acceleration Probe (IMAP; McComas et al. 2018).
The small-scale features of the signal require higher degrees of spherical harmonics, or in other words, spherical harmonics up to some maximum degree can only represent variations over larger scales. In general, it is possible to fully describe the IBEX data using spherical harmonics. However, here we use only the low-degree spherical harmonics to reconstruct the GDF, which after subtraction from the total IBEX signal, provides an estimation of the ribbon flux. The methodology is described in Sections 2 and 3. Section 4 presents the separation results for the full solar cycle of IBEX observations. Discussion of the observed GDF and ribbon structures evolving over the solar cycle is described in Sections 5 and 6. The resulting separated fluxes are posted as derivative products, concurrently, with the paper and are intended for use in follow-up studies. The format of these products is discussed in the Appendix.

Methodology
While perfect separation of the IBEX ribbon signal from the background GDF is not possible, as we cannot determine the source of the individual ENAs measured by IBEX, ENA sky maps can be decomposed into components with different spatial scales of the source in the celestial sphere. Such spatial scales do not necessarily correspond to the volume of the source region since the ENA observations provide line-of-sight (LOS) integrated ENA production. We utilize this decomposition to separate the GDF and the ribbon flux separately for each year and for each energy step. The separation method developed in this study is based on two postulates emerging from prior studies of the ribbon properties (Funsten et al. 2009b(Funsten et al. , 2015Fuselier et al. 2009a): I. The GDF is slowly spatially varying ENA emission, whereas the ribbon ENA emission is a notably sharper feature with distinctly higher flux gradients. II. The GDF is ubiquitous across the entire sky map, but the ribbon emission is a narrow feature that lies within a localized region of the sky. The width of this region is smaller than the scales of the GDF spatial variation.
These two postulates represent key differences between the ENA fluxes of the GDF and ribbon that we exploit using spherical harmonics to separate the ribbon and GDF fluxes without a priori model bias. Additionally, because each sky map is analyzed independently, it allows us to explore ribbon variability across energies and over time. This section provides the methodology of the separation and how the uncertainties are estimated. The robustness of the method and the postulates are further discussed in Section 3.

IBEX Data
Standard IBEX sky maps are obtained by summing the ENA observations in 1800 pixels with 6°× 6°resolution in the ecliptic coordinate system. The centers of the pixels are located at ecliptic longitudes λ p = 3°+ 6°× p, where p = 0, 1,...,59, and latitudes β q = − 87°+ 6°× q, where q = 0, 1,...,29. IBEX data are sometimes presented using different coordinate systems (e.g., galactic coordinates or ribbon-centered coordinates). Nevertheless, the ecliptic coordinates are the most natural due to the IBEX observational strategy in which the instrument boresight approximately follows the great circles perpendicular to the spacecraft spin axis orientation, which is kept close to the ecliptic plane and repointed once or twice every IBEX orbit around Earth to maintain the spacecraft axis close to the Sun (McComas et al. 2011).
The IBEX pixels can be enumerated by a single integer k defined as k(p, q) = p + 60q, where k k 0, 1 ,..., max = with k 1799 max = . A reverse transformation is p k k mod 60 ( ) = and r k k 60 ( ) ⌊ ⌋ = , where x ⌊ ⌋ denotes the floor function returning the greatest integer less than or equal to x. Using this convention, each IBEX map can be represented as a nonnegative vector j j k k Even though the IBEX data releases do not report data uncertainties correlations, it is helpful to define the covariance matrix Σ j = diag(σ 2 ), where diag(x) denotes a diagonal matrix with the vector values x placed at the matrix diagonal. Note that this matrix has variances (i.e., squared uncertainties) on its diagonal. The uncertainty in each pixel is dominated by statistical uncertainties, and therefore the uncertainty correlations are neglected.
This study uses the sky maps of ENA fluxes obtained from the IBEX-Hi observations covering the full solar cycle and released in Data Release 16 7 (McComas et al. 2020). The data were collected over 11 yrs from 2009 through 2019. We follow the year assignment provided by those authors. Moreover, we use ram data corrected for the Compton-Getting effect and the ENA survival probability from the outer heliosphere to their measurement at IBEX. The developed procedure may be easily applied to any IBEX data in this format, and some of these maps are included in the posted derivative products (see Appendix). We use ENA flux maps interpolated to the central energy of five electrostatic analyzer (ESA) energy steps from ESA 2 to ESA 6. Therefore, we analyze a total of 55 annual allsky ENA flux maps. We also include the time-combined maps covering the entire observation period. The procedure described in this study is independent for each analyzed map, therefore for simplicity, we drop the use of indices indicating energy steps or years further in this section.

Spherical Harmonics as an Orthonormal Basis
Spherical harmonics provide a complete orthonormal basis of smooth functions on the sphere (e.g., Arfken et al. 2012). Since the IBEX observed fluxes are real numbers, we use real spherical harmonics defined as where Y , ℓ m ( ) q j is the (complex) spherical harmonic of degree ℓ and order m. Note that the real spherical harmonics are denoted by both indices in the subscript. We follow the convention in which the Condon-Shortley phase is included in the definition of spherical harmonics. This is a standard convention used in popular mathematical libraries (e.g., SciPy, cmath) and in other popular programming languages (e.g., IDL, Mathematica). We utilize the ecliptic coordinate system directly, relating the azimuth with the ecliptic longitude (j = λ) and the polar angle with the ecliptic colatitude (θ = 90°− β, where β is the ecliptic latitude).
The spherical harmonics allow for representation of functions with the spatial scale of changes (defined as the distance from the function maximum to minimum) of Δ ℓ = 180°/ℓ, where ℓ is the degree of spherical harmonic. For each degree ℓ, there are 2ℓ + 1 orders of spherical harmonics numbered by m = − ℓ, − ℓ + 1,...,ℓ. Consequently, for a maximum degree of spherical harmonics equal ℓ max , there are ℓ 1 max 2 ( ) + spherical harmonics describing the function. Following Postulate I, we assume that the GDF component of an IBEX map can be represented by function j , ( ) q j defined in the sky as a sum of spherical harmonics up to degree ℓ max : where the coefficients can be calculated as: For ℓ max  ¥, any smooth function can be completely restored by this procedure, i.e., j j , , . However, for a limited number of harmonics, function j , ( ) q j is just an approximation of function j , ( ) q j , which smooths variations with spatial scales smaller than ℓ max D . Since the IBEX flux is only known for each pixel, the function value j , ( ) q j for the integral given in Equation (3) is assumed constant j k in each pixel k, and Equation (3) can be rewritten as a sum over all pixels: where the surface area of pixel k is denoted as Ω k . Defining the elements. Therefore, they are elements of two different linear spaces, and the matrices U and Y are used to transform between them.

Iterative GDF Decomposition
Postulate II provides that the IBEX observed pixels can be split into two subsets, the first one, which we denote as set A, includes pixels that we assume do not contain the ENA emission from the ribbon, while set B (the ribbon mask) consists of pixels that may have the ribbon component. Of course, set B may also include pixels that do not include any ribbon flux, but we try to define these two sets so that set A has as many pixels as possible. For this study, we define set B from a procedure described in Section 2.5, which identifies part of the sky for which the GDF modeled using the spherical harmonic decomposition cannot explain the IBEX observed fluxes because the spatial scale of this variation is much smaller.
The procedure presented in Section 2.2 requires that we know the GDF in all the IBEX pixels. However, as long as the spatial scales of the ENA flux gradients of the GDF exceed the width of the ribbon mask, it is possible to reconstruct the GDF inside this mask through interpolation. Therefore, we construct the following iterative procedure to obtain the spherical harmonics decomposition of the GDF. We iterate starting with the lowest possible degree of the spherical harmonics ℓ 0 ¢ = and continue by including higher degree harmonics when the procedure converges. In the first iteration, we start with the following zeroth-order estimation of the GDF in all pixels: where j k k A á ñ Î is an average GDF over the pixels outside the ribbon mask. While this selection of the initial conditions significantly speeds up further calculations, we test other possibilities. For example, we alternatively assume 0 flux in the set B or used actual signal values also inside this set. All variants give the same final results, and we conclude that our method does not depend on selecting this initial condition. In each iteration, we calculate the coefficients of the spherical harmonics: that are further used to determine the next order spherical harmonics decomposition: the spherical harmonic degree of ℓ¢. This decomposition is used to fill the ribbon pixels in the next iteration: We iterate Equations (9)-(11) until the following sum converges: typical value of the sum given in Equation (12) is of the order of the number of points in the set A. Therefore, the required convergency criterion is much more stringent than the range considered for the uncertainty derived from χ 2 statistics. After convergence for a given degree ℓ¢ of the spherical harmonics is achieved, we increase this degree unless it is equal to the highest degree adopted for the study ℓ ℓ max ( ) ¢ = , which is discussed later. For ℓ 0 ¢ ¹ , the starting approximation is defined as follows: where i ℓ 1 ¢-denotes the iteration for which the convergency was obtained for the maximum degree of ℓ 1 ¢ -. After the convergence is achieved for the maximum degree ℓ ℓ max ¢ = (see Section 3), we stop the procedure. The final coefficients of the GDF decomposition are given by Equation (9): , and the reconstructed GDF for the IBEX pixelization is obtained from Equation (10) This reconstruction provides an estimation of the GDF inside the ribbon mask and over the rest of the sky. If Postulate I was strictly fulfilled, we would expect that the difference between the measurements and the reconstruction outside the ribbon should only represent statistical fluctuations of the measurements. However, it is possible that small-scale ENA emission structures outside the ribbon and superposed with the GDF may exist. As we cannot uniquely identify small-scale ENA emission structures within the ribbon region, we estimate the flux magnitude from these structures outside the ribbon mask and include them in the uncertainty system. Since the timecombined maps show the smallest statistical fluctuations, we use the time-averaged contribution and apply it to all yearly maps with the same central energy. The variance of this contribution is calculated for each energy step as: where |A| denotes the cardinality of set A, i.e., number of points outside the ribbon mask. The sum includes only the points outside of the mask, as it represents structures other than the ribbon. Thus, the sum in the bracket would be approximately equal to 1 if the spherical harmonic decomposition describes the observed GDF fully. However, we find that this sum is always above 1, and the excess over 1 is used to estimate the variance h 2 . Equation (14) applies directly to the uncertainty of the GDF reconstruction; while it does not inform on the possible presence or absence of similar emission features of the ribbon, it is used as a component of uncertainty in the separated ribbon flux.

The Ribbon Flux and Uncertainty Estimation
Estimation of the ribbon flux using the obtained reconstructed GDF is straightforward, as we simply assume that the ribbon flux (r) is given as a difference of the IBEX observed flux (j) and the reconstructed GDF g  ( ): This equation may give a negative estimate of the ribbon flux both outside and inside the ribbon. However, in general, the distribution of this difference outside the ribbon should be statistically consistent with zero, which we verify in Section 3.
Uncertainties of the spherical harmonic coefficients obtained from Equation (9) can be calculated using uncertainty propagation. The covariance matrix of these coefficients is given by the following formula: where Σ g is the covariance matrix of the GDF. Based on Equation (11), within set A, g is equal to the observed ENA flux, for which the uncertainty is given simply as σ i -square root of the flux variance provided in the IBEX data release.
Since the values inside set B are reconstructed from the decomposition, their uncertainties are assessed as a standard deviation of the observed ENA fluxes in the considered IBEX map. In other words, the covariance matrix Σ g is equal to the original diagonal data covariance matrix Σ j with diagonal elements corresponding to pixels in set B replaced by the variance of ENA fluxes in the considered map. Even though the covariance matrix Σ g is diagonal, the covariance matrix of the spherical harmonic coefficients Σ c is not diagonal, and the correlations between them must be considered in further calculations.
The covariance matrix of the reconstructed GDF can be propagated from the spherical harmonic coefficients with the following formula: The smallest structures that can be reproduced by spherical harmonics up to the degree ℓ max is ℓ 180 .
ℓ max max D =  Therefore, smaller structures in GDF that are not reconstructed should be considered an additional uncertainty component. The magnitude of the variance is calculated as given in Equation (14). Still, due to possible angular correlations of this component, the covariance matrix should account for this because possible averaging over neighboring pixels should not assume that the related uncertainties in these pixels are independent. Hence, we estimate the related covariance matrix elements as follows: where h 2 is the expected variance of small-scale structures as given by Equation (14), k k , q ¢ is the angular distance between pixels k and k¢, and ℓ 180 1 is the largest spatial scale of structures that cannot be reconstructed for selected ℓ max .
Finally, the covariance matrix of the estimated ribbon flux is given as follows: This formula combines the contribution of original observational uncertainties Σ j , the uncertainties of the spherical harmonic decomposition g S , and the contributions from possible small-scale structures in the GDF Σ h . Note that the covariance matrix includes data variance (i.e., squared uncertainty) on the diagonal. Therefore, Equation (19) corresponds to a sum of squares of uncertainties from respective components.

Finding the Ribbon Mask
The part of the sky covered by the ribbon is estimated as the region where the ribbon flux is significantly larger than possible statistical fluctuation resulting from a counting statistic of the IBEX-Hi observations. Since observation times included in the IBEX maps must meet several quality criteria, the ENA fluxes in some pixels have poor statistics and, therefore, large uncertainties. Consequently, we also exclude pixels with fluxes below the estimated uncertainty, i.e., where the signal-to-noise ratio is below 1. This criterion also removes pixels without reported fluxes, for which the flux and variations are given as zero in the IBEX data releases. These pixels are excluded similarly to the ribbon mask, but unlike the ribbon mask, they are determined on a map-by-map basis. While this method may preferentially remove pixels with lower fluxes over those with higher fluxes with similar exposure times, it additionally provides that the number of counts is high enough to use Gaussian approximation of the flux uncertainty in these pixels. Nevertheless, the number of pixels removed by this criterion is much smaller than the number of pixels removed by the main ribbon mask, and thus this selection does not bias our results.
The main ribbon mask is obtained from an iterative procedure that starts with the mask that excludes only the pixels with low signal-to-noise ratios as described above. Then, the spherical harmonic decomposition is performed as described in Section 2.3 for the time-combined maps for each energy range. After this decomposition is obtained, the sum of the normalized residual signal is calculated for each pixel: where e enumerates central energies from ESA 2 to ESA 6. For the next iteration mask, we select pixels with the highest value of this sum. We start with 200 pixels (out of 1800 pixels in the IBEX maps) containing the highest normalized residuals and keep only contiguous regions covering at least 30 pixels. This criterion means that we exclude from the ribbon mask pixels that are scattered over the map. The decomposition is repeated with the new mask, and the procedure is iterated until the mask converges. The mask is determined based on the timecombined maps to minimize the impact of statistical variations and thus assumes that while the intensity of the ribbon may vary over time (McComas et al. 2014a(McComas et al. , 2017, the location does not ). Since the ribbon may not be limited to 200 pixels in the sky, we increase the number of pixels considered in the mask by 20 in subsequent iterations.
The iteration stops when the mask is one compact region and less than 15 out of 20 is added to the largest compact mask, i.e., a significant portion of the newly added pool of pixels is outside the main ribbon region. We find a separate mask for each maximum degree of spherical harmonics ℓ max applied in our procedure, but the same mask is used for all energy steps. Note that the variation of the ribbon position within the IBEX energy range does not significantly widen this region. Indeed, in all considered cases, the ribbon mask is narrow enough to meet the requirement given in Postulate II. The number of pixels included in the ribbon mask decreases with the maximum degree ℓ max , and is equal to 572, 423, 404, 397, 396, and 299 for ℓ max from 1 to 6, respectively. The mask, together with the value of normalized residuals, is shown in Figure 1. The significant decrease in the ribbon mask size between ℓ 1 max = and ℓ 2 max = suggests that at least the maximum order of 2 is needed to reconstruct the GDF. Additionally, the change between ℓ 5 max = and ℓ 6 max = shows that spherical harmonics of degree 6 reconstruct some of the ribbon flux near the upwind direction. The spherical harmonics are supposed to model the GDF only; therefore, the optimal maximum degree is likely between 2 and 5 (see further discussion in Section 3). After the ribbon mask is determined based on the time-combined maps, the procedure is performed for all considered years and energies, but the mask is not changed in these iterations.

Maximum Degree of Spherical Harmonics
To select the maximum degree ℓ max that is optimal for the separation of the IBEX ribbon, we calculate the normalized spectrum of spherical harmonic coefficients obtained from the time-combined maps without applying any mask, i.e., accounting for the ENA fluxes from the entire sky. The resulting spectrum is presented in Figure 2. The spectra show for all energy steps, except the highest, local minima for ℓ = 4. Therefore, the structures with spatial scales of Δ 4 = 45°have a smaller amplitude than those reproduced by spherical harmonics of degree ℓ = 5 (Δ 5 = 36°). This finding quantifies the observation that the ribbon structures are narrower than those in the GDF (Postulate I). The GDF should be well represented by spherical harmonics with ℓ 3, while the ribbon requires ℓ 5. Consequently, to reproduce most of the GDF, we need spherical harmonics up to ℓ 3 max = .
To confirm the above finding, the methodology presented in Section 2 is applied with ℓ max ranging from 1 to 6. Figures 3 and 4 present the GDF decomposition and the resulting ribbon flux for the time-combined maps for central energies 1.1 and 2.7 keV (corresponding to ESA steps 3 and 5), respectively, with ℓ max ranging from 1 to 6. The figures also show the residual flux outside the ribbon mask, defined as the difference between the observed flux and the reconstructed GDF. Therefore, the ribbon and residual fluxes are both obtained from Equation (15); the ribbon denotes the signal inside the ribbon mask (set B), while outside (set A), it is called here the residual flux, which includes the statistical variation and smallscale structures of the GDF. The statistical significance, denoted as signal-to-noise ratio (S/N), is the ratio of the ribbon flux or residual flux to the uncertainty. For the ribbon, the uncertainty is obtained from the ribbon flux covariance matrix given in Equation (19), and for the residual signal, only the data flux variance is used.
The GDF reconstructed using spherical harmonics with ℓ 1 max = does not show the main features of the IBEX signal, such as the increased flux in the nose and tail regions or the two lobes of lower emission in the port and starboard sides. This conclusion agrees with the discussion from Section 2.5, which suggested that this degree is too low to represent the GDF. Oppositely, for ℓ 6 max = , the spherical harmonic decomposition suggests the maximum GDF emission is within the ribbon region, and the structures seem to be clearly aligned with the  ribbon. The maps for ℓ max from 2 to 5 show a gradual increase of the smaller-scale features that can be included in the GDF, but at the same time, the reconstructed signal becomes aligned with the ribbon, suggesting that the decomposition reconstructs part of the ribbon signal.
Selection of the optimal ℓ max is a purposeful balance that distinguishes the ribbon from the GDF. While a higher value of ℓ max more accurately captures the sharp, smaller-scale structure of the ribbon, a lower value of ℓ max within this range ensures that less ribbon flux is incorrectly identified as GDF and also results in a larger ribbon mask, which can better accommodate possible temporal changes of the ribbon location over the time span of the annual maps. Inspection of Figures 3 and 4 indicates that ℓ 3 max = provides a good balance as these maps do not show significant structures aligned with the ribbon mask but can capture most of the GDF structure. Moreover, this selection allows the margin to include possible larger spatial features within the ribbon flux that are larger than those characterized by spherical harmonics with ℓ 5 max = , i.e., where the ribbon-related maximum the power spectrum is observed. This finding is consistent with the observation based on the spectrum of spherical harmonics presented in Figure 2. Table 1 shows the maximum GDF and ribbon flux obtained from the decomposition of the time-combined maps and the expected magnitude of small-scale variations of the GDF, as estimated from Equation (14) for all considered ℓ max . The table shows that the magnitude of small-scale variations decreases for higher ℓ max , except for ℓ 6 max = and central energies of 0.7 and 1.1 keV, for which it increases because the ribbon mask does not cover the entire ribbon region.
To verify that the developed uncertainty system is complete, we check whether the residual signal outside the ribbon mask, where we do not expect any ribbon flux, is consistent with no ribbon flux. For this purpose, we calculate the following χ 2 statistics:  (19). The sums are over pixels outside the ribbon mask (set A) and over all the energy steps and years of IBEX observations (marked as a sum over the maps). If the system is complete, the χ 2 value should be approximately equal to the number of degrees of freedom. Figure 5 shows the histograms of normalized residuals for these two χ 2 statistics and the values of the reduced statistics 2 c n , i.e., divided by the number of degrees of freedom (ν = N − K ), which are expected to be close to 1 for a complete uncertainty system. The number of data points is a sum of the number of points outside the mask for all maps (N = ∑ maps |A|), and the number of parameters is a product of the numbers of energy steps, years of observations, and spherical harmonic coefficients: K ℓ energy steps years 1 Statistic j 2 ( ) c S n significantly exceeds one but decreases for higher degrees of spherical harmonics from 4.08 for ℓ 1 max = to 1.22 for ℓ 6 max = . Therefore, higher degrees of spherical harmonics are needed to fully represent the GDF signal. However, the statistic r 2 ( ) c S n is close to one (0.86-1.06) for all ℓ max . Moreover, the normalized residuals with this uncertainty system r ( ) S approximately follow the normal distribution with standard deviation equal to 1. Thus, the uncertainty system, including the statistical distribution of residual signal from the time-combined maps, is complete regardless of the maximum spherical harmonics degree. The selection of ℓ max should be thus driven by other factors.
Based on the above arguments, we conclude that the ribbon separation from the IBEX data using the spherical harmonic decomposition of the GDF is possible. We use ℓ 3 max = as the optimal maximum degree of spherical harmonics describing the GDF. The GDF includes some smaller-scale components, whose magnitudes are provided in Table 1 and included in the overall uncertainty system. While Postulate I from Section 2 is not strictly fulfilled, the magnitude of the small-scale structures is small enough that even including them as an uncertainty allows for studies of the ribbon flux. For the optimum value ℓ 3, max = the width of the ribbon mask is small enough to reconstruct the GDF within the mask, as requested in Postulate II.

Results: GDF and Ribbon Flux
The results of the procedure described in Section 2 with the maximum degree ℓ 3 max = are presented in Figures 6-10 McComas et al. 2011). Nevertheless, the structures reconstructed in this region should be interpreted with caution. As indicated in Section 2.5, we exclude two types of pixels from finding the spherical harmonic decomposition of the GDF. First, we do not use pixels in the ribbon mask, which is the same for all maps (the white outline in the figures). Moreover, we exclude pixels with high relative uncertainty (the yellow outline), different for each map. The distribution of poor statistic pixels is not random. Since the IBEX-Hi boresight follows the great circle perpendicular to the IBEX spin axis, which is located close to the ecliptic plane and repointed every few days, the vertical groups of pixels in yearly maps correspond to the same observation periods. Measurements with periods of high backgrounds are systematically observed seasonally during poor heliospheric viewing (e.g., when the instrument field of view is obstructed by the magnetosphere or a majority of the IBEX orbit lies inside the magnetosphere) or episodically (e.g., during the passage of solar wind transients).
The GDF presented in the second column is reconstructed from the spherical harmonic decomposition, and therefore, it shows a signal with spatial scales of variation larger than Δ 3 = 60°. However, the amplitude of smaller-scale structures (represented by h in Table 1) is ∼10-20 times smaller than the magnitude of the GDF. Or in other words, the represented reconstruction should account for 90%-95% of the distributed ENA emission. Therefore, even though we do not expect a significant ribbon flux outside the ribbon mask region, we suggest using the entire decomposition map to represent the GDF rather than just replacing the observed pixel flux values within the ribbon mask with the reconstructed flux. Otherwise, the constructed data set would show different statistical properties inside and outside of the mask. Moreover, the reconstructed GDF does not include most of the statistical variation existing in the IBEX observations. Differently from the GDF, the ribbon flux, shown in the third column of these figures and given as a difference between the actual IBEX observations and the reconstructed GDF, includes the statistical variation of the original data. Therefore, these maps are much less smooth and show variation between pixels. Identical color bars for the first three columns in these figures enable direct comparison of these components. The ribbon flux is the highest inside of the ribbon mask, as expected. To estimate the statistical significance of the ribbon signal, we divide it by the square root of the diagonal elements of the complete uncertainty matrix; these are shown in the fourth column. This column indicates the statistical significance of the observed fluxes and confirms that most pixels with S/N > 3σ are also inside the ribbon region. As discussed in Section 3, the distribution of the normalized residuals outside of the ribbon does not show significant departures from the expected distribution.
Finally, the two last columns show the residual signal outside of the ribbon region and the residual signal normalized by the standard deviation of the original data. Note that we do not include uncertainties of the GDF decomposition or contribution from small-scale structures for the last column.
The first rows in Figures 6-10 show the results from the time-combined IBEX maps. These maps have the lowest statistical uncertainties and thus allow for a special evaluation of the separation process with high S/N information. Even within the ribbon mask, the separated ribbon is a comparatively narrow feature that appears much broader and brighter in the measured sky maps due to the superposition of the underlying emission structures in the GDF. This fact is noticeable especially close to the upwind direction where a local maximum in the GDF partially overlaps with the ribbon region. Due to the minimization of the statistical uncertainties, the time-combined maps have the highest S/N for the ribbon map. At the same time, the residual flux is much dimmer in the combined maps compared to the individual years, which indicates that most of the observed variation is statistical.
While the significance of the residual signal is slightly enhanced in the combined maps, especially in the two highest energy steps, it is not as much pronounced as the ribbon flux (note different scales in the fourth and sixth columns).
Nevertheless, some small-scale structures that cannot be reproduced by the spherical harmonic decomposition are observed close to the tail direction, especially in the highest energy steps. The two maxima south and north from the downwind direction appear to have steeper gradients than can be reproduced by spherical harmonics with ℓ 3. These features are likely related to the latitudinal structure of the solar wind confined in the heliotail (Zirnstein et al. , 2020Baliukin et al. 2020;Kornbleuth et al. 2020).

The Ribbon and GDF Evolving with the Solar Cycle
Observations in the two lowest energy steps (Figures 6 and  7) show similar emission structures within the GDF, which only slightly change over the solar cycle. These maps show systematic enhancements of the ENA emission close to the upwind direction as well as close to the downwind direction. The region of enhanced emission in the tail direction moves during the solar cycle from a direction slightly to the starboard side of the heliosphere during the earlier years to the downwind direction in the later years. However, this feature may result from poor statistics due to the magnetospheric obstruction, as discussed in Section 4. The time changes of the GDF for these two energy steps are relatively small. However, the ribbon flux is slightly more prominent in the first three years; later, it decreases, and recently it partially recovers but does not reach the magnitude observed in 2009.
Observations in energy steps 1.7 and 2.7 keV (Figures 8 and  9) show much more significant time changes in the GDF, but the changes are not uniform in the sky. Especially for energy 1.7 keV, the enhanced ENA emission close to the upwind direction shows a decrease followed by an increase of the flux, while the enhancement in the downwind direction (i.e., on the left and right edge of these maps) monotonically decreases over the covered period. The ribbon flux also seems to reduce over time, except for a slight increase south from the upwind direction in the most recent years. However, the ribbon flux in a region closer to the North ecliptic pole is much smaller in 2019 compared to 2009. At the highest energy step, shown in Figure 10, the ribbon flux is dim compared to the GDF. In previous studies (Funsten et al. , 2015Dayeh et al. 2019), the geometrical width of the ribbon at this energy is interpreted to be much wider than at lower energies; however, the approaches used in these studies included model bias and also were not optimized to identify or quantify a dim ribbon superposed on brighter GDF and the enhanced emission in the nose region. It is also possible that the decomposition with spherical harmonics up to ℓ 3 max = can reproduce some of the ribbon flux, yielding a dimmer ribbon after subtraction of the reconstructed GDF. However, the general structure of the GDF does not appear to be aligned with the ribbon but resembles the GDF reconstructed in the other The methodology presented in this paper allows us to track the time evolution of the GDF and ribbon flux separately inside the ribbon region. Figure 11 shows the time series of the flux averaged over the six regions covered by the IBEX ribbon. First, we select pixels with centers falling within ±9°from the circle centered at ecliptic (λ, β) = (218°.33, 40°.38) and radius 74°.81, which corresponds to the weighted mean position as determined by Dayeh et al. (2019). Then, these pixels are split into six regions Ri, i = 1, 2,...6, based on the ecliptic longitude of pixel centers falling inside the respective ranges λ ä (24°+ i × 30°, 54°+ i × 30°). The regions are marked in the top panel of Figure 11. Since the IBEX signal significantly varies among pixels within these regions, the flux is only included in the plot if at least 85% of the selected pixels in the source maps are covered in a given year. This criterion excludes data points in 2009-2010 (region R6), 2015 (regions R4 and R5), and 2016 (region R5). The combination of the fluxes inside the considered regions results in small statistical uncertainties. However, the presented uncertainty bars for the total flux and the ribbon does not represent variability inside of these regions. For the GDF, we include results obtained for reconstructions with ℓ 2 max = and 4, but only uncertainties for the optimal reconstruction with ℓ 3 max = are shown. The combination of the GDF and ribbon flux takes into account correlations between uncertainties of individual pixels.
The plot shows that contributions of the GDF and ribbon fluxes evolve differently. The two lowest energy steps show smaller changes over time, and the GDF and ribbon flux amplitudes are comparable. However, as the energy increases, the evolutions of the GDF and ribbon become much different. The amplitudes are much stronger for the GDF and show a significant increase in the two latest years, mostly due to the heliosphere response to the increase in the solar wind dynamic pressure in 2014 as discussed by McComas et al. (2019McComas et al. ( , 2020. During this time, we do not see a response to this pulse from the ribbon in most of these observations, except for a slight increase in energy step 4.3 keV in region R5. This is a part of the ribbon close to the upwind direction, which should show the first response since the heliopause is relatively close to this position. Nevertheless, the response is delayed compared to the GDF response. This indicates that the GDF and the ribbon ENAs are produced by two different processes, likely in different source regions.

Maximum Signal and Pressure Positions
The ENA emission structures of the GDF result from the global structure and thermodynamic properties of the heliospheric plasma. The IBEX ribbon is a unique discovery reflecting the large-scale draping of the interstellar magnetic field around the heliosphere but acts as an obstruction when attempting to interpret the structure of the GDF behind the ribbon. Nevertheless, with the careful separation of the ribbon, the underlying GDF can inform us of the plasma pressure in the heliosheath. Using their ribbon separation technique, Schwadron et al. (2014) calculated the stationary LOS-integrated pressure of plasma proton from GDF separated flux (see their Equation (1)). These calculations revealed that the maximum of the LOS-integrated pressure is not aligned with the upwind direction but is located about 20°south from this direction. McComas & Schwadron (2014) advocated that this explains the plasma flow velocities observed by Voyager 2 as being directed away from this maximum pressure point and not from the upwind direction. However, Grygorczuk et al. (2015) showed that this plasma flow vector is located inside of the hydrogen deflection plane. Since Voyager 2 was (coincidentally) located very close to this plane for most of the time in the inner heliosheath, this finding suggests flow away from the upwind direction and not from the maximum pressure as indicated by the pressure maximum derived by Schwadron et al. (2014).
The two main external pressures shaping the heliosphere come from the flow of the VLISM relative to the Sun and the local interstellar magnetic field (e.g., Parker 1961). While the solar wind is characterized by its slow-fast latitudinal structure around solar minimum, the dynamic pressure of the solar wind is uniform over all heliographic latitudes, as shown by the Ulysses observations (McComas et al. 2008). Consequently, the variation of the solar wind speed in latitude should not significantly impact the product of the LOS and pressure, which is estimated here. Figure 12 presents the LOS-integrated pressure obtained from the time-combined reconstructed GDF flux. The maximum pressure position is marked with a dark blue ellipse. The ellipse represents the 3σ uncertainty region obtained based on the developed uncertainty system. The maximum pressure direction is centered at ecliptic , 266 . 7 1 . 6, 8 . 2 0 . 9 ( ) ( ) l b =    -   and is aligned with the so-called B-V plane, defined by the directions of the ISN helium flow (Swaczyna et al. 2018) and the pristine interstellar magnetic field far from the heliosphere (Zirnstein et al. 2016b). This plane also includes the directions of the ISN hydrogen flow (Lallement et al. 2005), the Warm Breeze flow (Kubiak et al. 2016), and the IBEX ribbon center . Therefore, the plasma flow observed by Voyager 2 in the inner heliosheath is away from the maximum pressure (McComas & Schwadron 2014) and from the inflow direction (Grygorczuk et al. 2015). Furthermore, this observation strengthens the finding that the B-V plane is an approximate symmetry plane of the heliosphere.
The pressure shown in Figure 12 is obtained from the timecombined observations. Additionally, we calculated this pressure for the reconstructed GDF from individual years. Moreover, we also checked the position of maximum flux in the individual energy steps. Figure 13 shows positions of these directions, characterized by the angle from the upwind direction and the angular distance from the B-V plane. Negative and positive angular distances from the B-V planes correspond to the hemisphere containing the southern and northern ecliptic poles, respectively.
The figure shows that changes from year to year are statistically significant since they exceed the uncertainties. Zirnstein et al. (2018b) showed that changes in the solar wind Figure 12. Stationary LOS-integrated pressure of the proton plasma forming the ENA observed by IBEX, obtained from the time-combined GDF reconstructed maps. The following directions are marked with 3σ ellipses: the maximum pressure (this study, dark blue), the ISN He flow (Swaczyna et al. 2018; purple), the ISN H flow (Lallement et al. 2005; gray), the Warm Breeze (Kubiak et al. 2016;cyan), the local interstellar magnetic field (Zirnstein et al. 2016b;orange), and the IBEX ribbon center yellow). In addition, the approximate position of the ribbon is marked with gray dashed lines, and the B-V plane is marked with the white line. White arrows indicate positive and negative inclinations to the B-V plane, and the angular distance of the maximum pressure from the nose is marked with a green arrow. dynamic pressure can be observed in ENA fluxes with a delay of few years, which depends on the position in the sky. Therefore, the single-year ENA maps should not be interpreted as a snapshot of the heliosphere at a specified time. However, the combination of the entire solar cycle observation gives the best possible representation of average heliosphere conditions. The solar wind dynamic pressure pulse seen in this solar cycle may, however, affect this result (McComas et al. 2019). Nevertheless, note that the right panel of Figure 13 shows that over the first 5 yrs, the maximum pressure was shifted to the south, consistent with the result obtained by Schwadron et al. (2014) and interpreted by McComas & Schwadron (2014).
The maximum pressure is shifted from the upwind direction due to the interstellar magnetic field, which contributes to the interstellar pressure outside the heliopause (Fahr et al. 1986;Ratkiewicz 1992;Fuselier & Cairns 2013). To demonstrate this interpretation of the observations, we present a simple calculation of the maximum pressure offset first provided by Fahr et al. (1986Fahr et al. ( , 1988. They estimated the deflection of the maximum pressure from the upwind direction using a Newtonian approximation model, given as follows: ) r = ć m −3 (m p is the mass of a proton), as found by Zirnstein et al. (2016b), is 21.3 ± 1.3 km s −1 in the VLISM. With the interstellar flow speed estimation of 25.85 ± 0.40 km s −1 (Swaczyna et al. 2021), the Alfvénic Mach number is thus 1.21 ± 0.08. Using the inclination of the magnetic field θ BV = 39°.5 ± 0°.7 from Zirnstein et al. (2016b), the deflection given by Equation (23) is θ s = 18°.7 ± 2°.0, and agrees with the deflection of the maximum pressure found empirically in this study: 17°.2 ± 0°.5. Therefore, we conclude that the maximum pressure position obtained in this study is consistent with the expected ones based on the influence of the interstellar dynamic and magnetic pressure at the heliopause.

Summary
The separation of the ENA fluxes originating from the ribbon and the GDF in the IBEX-Hi observations is crucial to understanding the global heliosphere plasma structure and its dynamic interaction with the VLISM. The new approach used in this study exploits two key differences in the properties of the emission structures of the ribbon and of the GDF to clearly separate their superposed fluxes: the smaller characteristic spatial dimension of the ribbon compared to the GDF and the spatially localized emission of the ribbon relative to the global emission of the GDF. Outside of the ribbon mask (within which the ribbon flux resides), we find that the spherical harmonics decompositions of the GDF with spherical harmonics up to maximum degree of ℓ 3 max = allow for accurate reproduction of the global GDF. We also estimate the magnitude of small-scale structures within the GDF that would likely be better represented using higher degree spherical harmonics ( Table 1) for inclusion in the uncertainty calculated for the ribbon flux.
The ribbon flux is then estimated as a difference between the total observed ENA flux and the spherical harmonic decomposition of the GDF interpolated in the ribbon mask, which is determined solely based on the residual flux. Importantly, this method does not rely on a spatial model of the ribbon emission and thus is not susceptible to model bias. The ribbon and GDF emission structures derived in this study show distinctively different temporal variations, which confirms that these two components are likely created in different source regions. Specifically, the evolution of the ENA flux in the highest energy step in response to the increased solar wind dynamic pressure in 2014 is visible primarily in the GDF but not in the ribbon flux. Furthermore, only a slight increase in the ribbon flux can be observed in a part of the ribbon located close to the heliospheric nose in the most recent year and the highest energy steps, consistently, with results reported by McComas et al. (2020). This suggests that the ribbon source region is the closest near the nose but still is located beyond where the GDF originates, i.e., beyond the heliopause.
The GDF reconstructed from the spherical harmonic decomposition allowed us to find the stationary (with respect Figure 13. Position of the maximum pressure (black symbols) and maximum flux (with color-coded energy steps) relative to the upwind direction (left panel) and the B-V plane (right panel). The time series shows how the position changes with time, and the subpanels on the right part of each panel show the time-combined results. Missing data points correspond to times when the maximum in the proximity of the upwind direction was not present or was not statistically significant (see .