Towards Improved Heliosphere Sky Map Estimation with Theseus

Abstract

The Interstellar Boundary Explorer (IBEX) satellite has been in orbit since 2008 and detects energy-resolved energetic neutral atoms (ENAs) originating from the heliosphere. Different regions of the heliosphere generate ENAs at different rates. It is of scientific interest to take the data collected by IBEX and estimate spatial maps of heliospheric ENA rates (referred to as sky maps) at higher resolutions than before. These sky maps will subsequently be used to discern between competing theories of heliosphere properties that are not currently possible. The data IBEX collects present challenges to sky map estimation. The two primary challenges are noisy and irregularly spaced data collection and the IBEX instrumentation’s point spread function. In essence, the data collected by IBEX are both noisy and biased for the underlying sky map of inferential interest. In this article, we present a two-stage sky map estimation procedure called Theseus. In Stage 1, Theseus estimates a blurred sky map from the noisy and irregularly spaced data using an ensemble approach that leverages projection pursuit regression and generalized additive models. In Stage 2, Theseus deblurs the sky map by deconvolving the PSF with the blurred map using regularization. Unblurred sky map uncertainties are computed via bootstrapping. We compare Theseus to a method closely related to the one operationally used today by the IBEX Science Operation Center (ISOC) on both simulated and real data. Theseus outperforms ISOC in nearly every considered metric on simulated data, indicating that Theseus is an improvement over the current state of the art.

KEYWORDS:

• Deconvolution
• Ensemble models
• Generalized additive models
• Heliosphere
• Projection pursuit regression
• Regularization

1 Introduction

Launched in 2008, the Earth-orbiting Interstellar Boundary Explorer (IBEX) mission has been continuously providing energy-resolved measurements of energetic neutral atoms (ENAs) from the boundary between our solar system and what lies beyond (i.e., the boundary between the heliosphere and interstellar space) (McComas et al. 2009b). ENAs are fast moving particles with no charge. They are produced when a particle with a net electric charge in the heliosphere undergoes charge exchange with a particle with a net neutral charge in interstellar space. The particle that previously had a net electric charge but now has a net neutral charge is called an ENA. One property of an ENA is that it no longer interacts with magnetic fields and will travel in a ballistic trajectory (i.e., a nearly straight line at the energies being observed) from the time and location it was neutralized. Some of the ENAs created in the heliosphere travel in a ballistic trajectory toward Earth and can be detected by IBEX. This is what IBEX is counting: ENAs created in the heliosphere. By continuously counting energy-resolved ENAs, IBEX is able to create spatial maps of ENA rates at different energies over time, effectively creating dynamic images of the heliosphere. These maps of ENA rates called sky maps have been successfully used to deduce the global structure and dynamics of the heliosphere, leading to groundbreaking advances in our understanding of the interstellar boundary (see Table 1 in McComas et al. 2017; McComas et al. 2020, and studies since then).

The heliosphere is a large cavity shaped like an elongated egg that encapsulates our solar system.¹ The heliosphere travels through interstellar space like a boat traveling through water; the front of the boat traveling forward with the back of the boat following along. The compressed front of the heliosphere is called the nose and the elongated back is called the tail. The ENAs observed by IBEX primarily come from two different populations. The primary population is called the globally distributed flux (GDF). These particles arrive from all directions in the sky, but arrive at higher rates toward the nose and tail of the heliosphere (e.g., McComas et al. 2009a; Schwadron et al. 2014). The second ENA population is called the ribbon. It is a relatively narrow band of enhanced ENA emissions that form a nearly complete circular band across much of the sky (Fuselier et al. 2009). The ribbon was entirely unanticipated by any model or theory prior to its discovery by IBEX: a bona fide success for discovery science (McComas et al. 2009a).

As the ribbon was unknown prior to IBEX’s launch in 2008, theories to explain its origin are actively being developed. To date, over a dozen theories have been posited to explain the ribbon’s origin (see list in McComas, Lewis, and Schwadron (2014) and more recently, Zirnstein, Heerikhuisen, and Dayeh 2018; Zirnstein et al. 2019; Zirnstein, Dayeh, and Heerikhuisen 2021), the details of which are outside the scope of this article. What is relevant is that different ribbon origin theories can lead to different ribbon profiles or shapes. For example, Zirnstein et al. (2019) shows how two different ribbon origin theories—the Weak Scattering model and the Spatial Retention model—result in distinct ribbon shapes. If those shapes could be discerned with observations, that could add support to one model over the other. The current state-of-the-art sky map rendering approach used by the IBEX Science Operations Center (ISOC), however, produces maps in which the ribbon profile is not well enough resolved to discern between competing ribbon formation theories. Furthermore, finer-scale GDF structures are also not clearly resolved (e.g., enhancements that appear at higher energies in the tailward direction (Zirnstein et al. 2016b)). Improved sky map estimates are needed.

The principal instrument on IBEX for making the observations used to construct ENA sky maps is the IBEX-Hi ENA imager (Funsten et al. 2009). Prior to the launch of IBEX, there was no expectation that ENA sources as narrow as the ribbon would be present. The IBEX-Hi ENA imager has a roughly circular instantaneous field-of-view (FOV) having a ${6.5}^{°}$ full-width at half maximum (FWHM) angular resolution (Funsten et al. 2009). What that means is when IBEX-Hi detects an ENA, even though IBEX-Hi records the direction of the boresight of the FOV in the sky at the time of detection (i.e., the location in the sky IBEX-Hi was aimed at the time of detection), that does not mean the detected ENA originated from that same boresight location. The detected ENA could have originated from any location “near” IBEX-Hi’s boresight, where “near” is probabilistically defined by the IBEX-Hi point spread function (PSF): a peaked, roughly conical-shaped distribution that describes the likelihood of ENA detection as a function of angle off the boresight. In practice, the PSF creates a blurring effect on the sky maps; a blurring effect that should be accounted for and removed in order to resolve more granular sky map features (e.g., the ribbon’s fine structure). Currently, the ISOC sky maps do not account for the PSF or attempt to remove its blurring effect.

Another limitation of ISOC sky maps is resolution. In principle, a sky map is a continuous spatial map where every location in the spatial domain has a corresponding ENA rate. In practice, a sky map is a spatial map discretized into nonoverlapping pixels where every pixel has a corresponding ENA rate. ISOC sky maps are made at a fixed $6^{°}$ pixel resolution and are unable to be reliably created at finer resolutions, resulting in a pixelated esthetic. ISOC sky maps are shown in Figure 1 for reference.

Fig. 1 Binned direct event data ecliptic longitude and latitude locations (points) for ESA steps 2 and 5 (rows) corresponding to the 2013A sky map. There are 13,716 and 13,740 binned direct event data points for ESA steps 2 and 5, respectively. From left to right, the columns show the exposure time, number of direct events, background rate, ENA rate computed as the number of direct events/exposure time—background rate, the state-of-the-art ISOC sky map, and the novel Theseus sky map estimate. The ribbon can be seen in ESA step 5’s ENA rate plot as the streak of red and yellow points starting at longitude and latitude (60, 70), extending down and right to roughly (240, –40), and extending up and right to roughly (130, 45) and more clearly seen in the ISOC and Theseus sky maps.

Our approach to sky map estimation borrows strategies from spatial statistics and image deblurring. If there were no PSF, estimating a sky map would be a straight-forward spatial statistics problem. It would require estimating a latent spatial surface from noisy, irregularly spaced count data on a sphere. Hierarchical modeling is a common approach to this type of problem, where the data are modeled as conditionally independent given the unobserved spatial surface (e.g., Wikle, Berliner, and Cressie 1998). A link function is used to connect the latent spatial surface to the expectation of the observations. The spatial surface can be estimated with a number of techniques non-exhaustively including Gaussian processes (GPs) and variants thereof (e.g., Gramacy and Lee 2008; Gramacy and Haaland 2016; Sauer, Gramacy, and Higdon 2022; Katzfuss, Guinness, and Lawrence 2022), generalized additive models (Hastie and Tibshirani 1987), projection pursuit regression (Friedman and Stuetzle 1981), or multivariate adaptive regression splines (Friedman 1991). Uncertainty quantification of the latent spatial surface can be directly estimated from some methods (e.g., GPs) or can be estimated via bootstrapping (e.g., Castillo-Páez, Fernández-Casal, and García-Soidán 2019).

If the PSF were to be ignored, the estimated spatial surface would be overly blurry relative to the unblurred sky map of interest. Given a blurred sky map, removing the blurriness and recovering the sharp sky map is an image deblurring problem. Image deblurring is essentially an ill-posed inverse or deconvolution problem. At its core, there is not enough information in a blurred image to uniquely sharpen it, therefore, additional information must be added and/or assumptions must be made. Common approaches to image deblurring include iterative methods (e.g., Biemond, Lagendijk, and Mersereau 1990), inverse filtering (e.g., MacAdam 1970; Bojarczak and Lukasik 2011), regularization (e.g., Ma, Xu, and Zeng 2017), and more recently, deep learning (e.g., Zhang et al. 2022).

In this article we introduce a two-stage sky map estimation procedure called Theseus. In Stage 1, Theseus estimates the as-observed, or blurred, sky map from the noisy and irregularly spaced binned direct event data using an ensemble approach that leverages projection pursuit regression and generalized additive models. In Stage 2, Theseus deblurs the sky map by deconvolving the PSF with the blurred map using regularization. Sky map uncertainties are computed via bootstrapping. For Theseus, the sky map resolution is user-defined, allowing for more resolved sky maps than the fixed $6^{°}$ ISOC sky maps.

The contributions made in this article are 4-fold. (a) Never before have science operations estimated sky maps been compared to simulated sky maps, allowing us, for the first time, to quantitatively evaluate different sky map estimation procedures. (b) Theseus is the first procedure used to estimate sky maps that directly accounts for and removes the blurring effect of the PSF. (c) Theseus makes maps at a user-defined resolution and executes principled interpolation; the current state of the art makes sky maps at a fixed $6^{°}$ resolution and leaves “holes” in the sky maps where data are missing. (d) We show that Theseus sky maps are better than current ISOC sky maps when compared over a variety of evaluation metrics on a collection of simulated sky maps.

The article proceeds as follows. In Section 2, we describe the binned direct event data collected by IBEX-Hi. In Section 3, we describe the sky map making methodology, Theseus. In Section 4, we compare Theseus to the current state-of-the-art sky map estimation procedure used by the ISOC on simulated maps capturing different ribbon models. In Section 5, we fit Theseus and the ISOC method to real data collected by IBEX-Hi. Finally in Section 6, we provide a discussion of future work.

2 The IBEX-Hi Binned Direct Event Data

IBEX collects temporally-, spatially-, and energy-resolved ENA measurements from the heliosphere (i.e., IBEX images the entire sky). IBEX rotates at $\sim 4$ revolutions per minute around a Sun-pointing spin axis. The IBEX-Hi imager’s boresight is aligned perpendicular to the spacecraft spin axis, thus, over the course of a spin, IBEX-Hi images ENAs arriving from a roughly ${6.5}^{°}\times {360}^{°}$ circular band of the sky. IBEX orbits the Earth with an orbital period of $\sim 9$ days, and the spin axis is repointed toward the sun roughly every 4.5 days, leading to bands of sky viewing centered $\sim {4.5}^{°}$ apart.² Two different bands of sky viewing occur during an orbit, one on the ascending arc of the orbit, and one on the descending arc; thus, each band is referred to as an arc of data collection. The intensity of ENAs of heliospheric origin does not change appreciably over the course of months, and thus the average ENA rate arriving from a specific direction averaged over the course of an arc is considered a true measure of the ENA rate from that location in the sky. It takes IBEX roughly 180 days to image the whole sky. As such, sky maps are made in 6-month cadences.³ The maps corresponding roughly to the first part of the calendar year are labeled “A” maps; the latter half “B” maps.

The IBEX-Hi instrument is continuously detecting energy-resolved ENAs in the range of 500 eV to 6 keV. The instrument covers this range by stepping through six discrete energy settings of its electrostatic analyzer (ESA). ESA step 2 is the lowest step used for science, and ESA step 6 is the highest.⁴ Each IBEX-Hi measurement event is referred to as a direct event, and occurs when particle detectors within the instrument are triggered in a specified sequence. When a direct event is measured, numerous pieces of information are recorded including the time the direct event occurred, the spatial location in the sky the IBEX-Hi boresight was pointed, and the ESA setting IBEX-Hi was set to at the time of the direct event measurement. The direct events amount to energy-resolved, spatial-temporal point process data (González et al. 2016).

While direct events are the raw measurements made by IBEX-Hi, that is not the data product we work with. Instead, we work with a curated data product called the binned direct event data. The binned direct event data are in essence a partitioned, aggregated, and processed version of the direct events. The process to convert direct events to binned direct events is standard for the IBEX mission. Essentially, the direct events measured in an arc are partitioned into 360 $1^{°}$ bins for each ESA step. The direct events assigned to each ESA step/$1^{°}$ bin amount to temporal point process data, where the total number of direct events would ideally equal the number of ENAs measured by the IBEX-Hi instrument when set to a specified ESA step and pointed at a $1^{°}$ bin of the sky. Unfortunately, not all direct events are ENAs of heliospheric origin. A direct event could be a heliospheric ENA or it could be a background particle (i.e., a high-energy penetrating radiation particle, a non-heliospheric ENA, or another particle local to the spacecraft environment that enters the instrument aperture mimicing a heliospheric ENA). A differentiating behavior of background direct events are their degree of spatial uniformity and the time scale of their dynamics. Much of the background of local origin varies much more rapidly in time, and varies with viewing direction. A standard IBEX data processing step segments the temporal point process direct event data for each ESA step/$1^{°}$ bin and discards the direct events in the time segments where it is determined that heliospheric ENA direct events are either overwhelmed by background particles, background rates change too quickly, or background particles are associated with unusual local solar wind conditions. Generally speaking, the time segments that are kept in this IBEX-standard culling step are times when the background rate is relatively constant across locations, slowly varying in time, and the background rates are not appreciably larger than their corresponding ENA rates. The total duration of the time segments kept after the culling step is referred to as the binned direct event datum’s exposure time.

A binned direct event datum corresponds to a location in the sky (i.e., $1^{°}$ bin) and an ESA setting for a given 6-month map. A binned direct event datum includes the IBEX-Hi boresight direction, the exposure time, the number of measured direct events, and a background rate (which is another standard IBEX data product). The background rate is technically an estimate, but with small uncertainties. In this article, the background rate is treated as known.

For illustration, Figure 1 plots the exposure time, the number of direct events, the background rate, and an estimated ENA rate (number of direct events/exposure time—background rate) for ESA steps 2 and 5 for the 2013A sky map. Furthermore, the ISOC and Theseus sky maps are presented, estimated from the binned direct event data. The plots in Figure 1 and throughout this article are centered on longitude 265, with longitudes in reverse order (i.e., larger longitude values are on the left and smaller ones on the right, before centering). The centering on longitude 265 corresponds to a “nose-centered” frame, where the nose of the heliosphere is placed roughly in the center of the map. The plotting of longitudes in reverse order before centering is standard in astronomy, as the sky is being viewed from inside the heliosphere. There are 13,716 and 13,740 binned direct event data points for ESA steps 2 and 5, respectively. The exposure time varies by location, with exposure times nominally between 0 and 180 sec.⁵ Locations that have missing binned direct event data correspond to ESA steps/bins where the entire time series of direct events was culled. This creates gaps in the binned direct event data. Those gaps are correlated across ESA steps, as the reason for culling at one energy step is often the same at other energy steps. Roughly two-thirds of the possible collection time is discarded by this standard IBEX data processing step. In general, lower energy ESA steps tend to have more data discarded than higher energy ESA steps. For ESA steps 2 and 5 of the 2013A sky map, the number of direct events per bin is between 0 and 51. The number of measured direct events is correlated with exposure time; the larger the exposure time, the more measured direct events. The background rate varies between 0.03 and 0.12 background particles per second and is typically common for all bins within an ESA and arc. Finally, the estimated ENA rate varies by ESA step. The ribbon can be seen in the ENA Rate panel for ESA step 5 as the streak of red and yellow points roughly connecting the longitude and latitude points (60,70), (240, –40), and (130, 45) and more clearly seen in the ISOC and Theseus estimated sky maps. There is appreciable noise in the crude ENA rate estimates, with many negative estimated ENA rates (which are physically impossible).

3 Theseus Methodology

In this section, we outline Theseus—our sky map estimation methodology. In Section 3.1, we describe the data model, define terms and lay out notation. Theseus is a two-stage estimation procedure that, for a specified ESA and 6-month map (e.g., the 2013A, ESA 2 sky map), takes a binned direct event dataset D as input and outputs an estimated unblurred sky map (i.e., a sky map where the blurring effect of the PSF has been removed) and blurred sky map (i.e., a sky map where the blurring effect of the PSF has not been removed). The unblurred sky map is the map of scientific interest, while the blurred sky map is a necessary intermediate step to compute the unblurred sky map and is a necessary map for model diagnostics (more on this in Section 4.3). More specifically, in Stage 1, Theseus estimates a blurred sky map from a binned direct event dataset (Section 3.2). In Stage 2, Theseus estimates an unblurred sky map from the estimated blurred sky map in Stage 1 (Section 3.3).

Sky map uncertainties are estimated via percentile bootstrap intervals (e.g., Efron and Tibshirani 1994, chap. 13). As such, Theseus is run N_B times, resulting in N_B estimated unblurred and blurred sky maps corresponding to N_B different bootstrap binned direct event datasets ${\mathit{D}}_b$ for bootstrap sample $b=1,2,.\dots ,N_B$. In this article, N_B = 1000. ${\mathit{D}}_b$ is constructed in two steps. First, we perform stratified sampling with replacement where strata are defined as arcs. Equal probability sampling with replacement is performed on the rows of D within each stratum, where a row of D corresponds to a binned direct event datum (i.e., a point in Figure 1). The stratified sampling encourages spatial coverage in ${\mathit{D}}_b$ to be reflective of the spatial coverage in D. Second, we simulate the number of direct events from its generating model (see (1)). This step helps reduce some of the visual artifacts Theseus can produce in the sky map estimation (this is discussed more in Section 3.3). While each estimated unblurred and blurred sky map from Theseus is properly indexed by ESA, 6-month map, and bootstrap sample, those indices will be suppressed throughout the remainder of Section 3 for notational simplicity. In what follows, we describe how Theseus estimates an unblurred and blurred sky map for a specified ESA, 6-month map, and bootstrap sample.

Finally, for clarity, we explicitly list some of the notational conventions used throughout this article.

• Spatial quantities indexed by i will correspond to locations of binned direct event data while spatial quantities indexed by j will correspond to locations of sky map pixel centers (recall from Section 1 that, in practice, a sky map is defined by a discrete collection of mutually exclusive and exhaustive pixels).

• Non-bolded quantities will be used for scalars, while bolded quantities will be used for non-scalars (e.g., vectors and matrices).

• Numerous sky map-related quantities will be estimated throughout the article. Final estimates of quantities will be denoted with “hats” (e.g., $\widehat{\theta }$), while intermediate estimates needed to arrive at the final estimated quantities will be denoted with “dots” (e.g., $\dot{\theta }$) or “double dots” (e.g., $\ddot{\theta }$). The order in which quantities are computed going from data to final estimates is as follows:

  $\text{Data}\to \stackrel{.}{\mathit{\theta }}\to \stackrel{..}{\mathit{\theta }}\to \widehat{\mathit{\psi }}\to \stackrel{.}{\mathit{\delta }}\to \widehat{\mathit{\alpha }}\to \widehat{\mathit{\delta }}\to \widehat{\mathbf{\Theta }}\to \widehat{\mathit{\theta }}$

To lessen the notational burden, the supplementary materials provides a table enumerating the symbols and corresponding definitions used throughout this manuscript.

3.1 Data Model

Let ${\mathit{s}}_i=({\text{lon}}_i,{\text{lat}}_i)$ be the ecliptic longitude/latitude coordinates of location ${\mathit{s}}_i$ for binned direct event datum $i=1,2,\dots ,N_o$ where N_o is the number of binned direct event data observations, ${\text{lon}}_i\in [0,360)$, and ${\text{lat}}_i\in [-90,90]$. For the basic data model, we have (1) $$y({\mathit{s}}_i) | t({\mathit{s}}_i),\theta ({\mathit{s}}_i),b({\mathit{s}}_i)\sim \text{Poisson}(t({\mathit{s}}_i)[\theta ({\mathit{s}}_i)+b({\mathit{s}}_i)])$$(1) where $y({\mathit{s}}_i)\in 0,1,2,\dots$ is the number of direct events (ENAs + background particles) counted by IBEX-Hi, $t({\mathit{s}}_i)>0$ is the known exposure time in seconds, $b({\mathit{s}}_i)\ge 0$ is the known background rate (background particles/second), and $\theta ({\mathit{s}}_i)\ge 0$ is the unknown ENA rate (ENAs/second) corresponding to the blurred sky map. As expressed in (1), the number of direct events counted by IBEX-Hi when being pointed at location ${\mathit{s}}_i$ for $t({\mathit{s}}_i)$ seconds is modeled as an ENA rate $\theta ({\mathit{s}}_i)$ plus a background rate $b({\mathit{s}}_i)$ times the exposure time $t({\mathit{s}}_i)$. On average, $t({\mathit{s}}_i)b({\mathit{s}}_i)$ direct events are background particles and $t({\mathit{s}}_i)\theta ({\mathit{s}}_i)$ are ENAs.

While $\theta ({\mathit{s}}_i)$ is the ENA rate corresponding to the IBEX-Hi instrument at location ${\mathit{s}}_i$, $\theta ({\mathit{s}}_i)$ is not the ENA rate of the heliosphere corresponding to location ${\mathit{s}}_i$ (i.e., is not the value of the unblurred sky map at location ${\mathit{s}}_i$). This is because not all ENAs counted when the IBEX-Hi instrument is pointed at location ${\mathit{s}}_i$ originated from location ${\mathit{s}}_i$ due to the IBEX-Hi PSF. The blurred sky map at location ${\mathit{s}}_i$ is related to the unblurred sky map $\mathbf{\Theta }$ as follows: (2) $$\theta ({\mathit{s}}_i)=\mathit{K}({\mathit{s}}_i)\mathbf{\Theta }$$(2) where $\mathbf{\Theta }$ is an unknown, nonnegative $N_p\times 1$ column vector of ENA rates, N_p is the number of mutually exclusive and exhaustive pixels the sky map has been partitioned into⁶, and $\mathit{K}({\mathit{s}}_i)$ is a known $1\times N_p$ row vector. $\mathit{K}({\mathit{s}}_i)$ is the blurring operator encoding the properties of the IBEX-Hi PSF corresponding to the boresight location ${\mathit{s}}_i$. $\mathit{K}({\mathit{s}}_i)$ has two properties: every entry of $\mathit{K}({\mathit{s}}_i)$ is nonnegative (i.e., $\mathit{K}{({\mathit{s}}_i)}_j\ge 0$ for all $j\in 1,2,\dots ,N_p$) and $\mathit{K}({\mathit{s}}_i)$ sums to 1 (i.e., ${\sum }_{j=1}^{N_p}\mathit{K}{({\mathit{s}}_i)}_j=1$). $\mathit{K}({\mathit{s}}_i)$ probabilistically describes the mapping between the true ENA origin location to the observed ENA location ${\mathit{s}}_i$. Specifically, $\mathit{K}{({\mathit{s}}_i)}_j$ is the probability an ENA originated from within pixel j of the sky map given it was counted by the IBEX-Hi instrument when pointed at location ${\mathit{s}}_i$. The PSF is a property of the IBEX-Hi instrument and must be accounted for in the statistical model. How $\mathit{K}({\mathit{s}}_i)$ is constructed is detailed in the supplementary materials.

To reiterate, the unblurred sky map $\mathbf{\Theta }$ is the quantity of inferential interest. Our goal is to estimate $\mathbf{\Theta }$, with uncertainty, given $\mathit{y}=\{y({\mathit{s}}_1),y({\mathit{s}}_2),\dots ,y({\mathit{s}}_{N_o})\},\mathit{t}=\{t({\mathit{s}}_1),t({\mathit{s}}_2),\dots ,t({\mathit{s}}_{N_o})\}$, and $\mathit{b}=\{b({\mathit{s}}_1),b({\mathit{s}}_2),\dots ,b({\mathit{s}}_{N_o})\}$. We estimate $\mathbf{\Theta }$ in a two stage procedure. In Stage 1 (Section 3.2), we estimate the blurred sky map from the binned direct event data. In Stage 2 (Section 3.3), we estimate the unblurred sky map given the Stage 1 estimate of the blurred sky map.

3.2 Theseus Stage 1: Estimating the Blurred Sky Map $\mathit{\theta }$

The first stage of Theseus is to estimate the blurred sky map. Let ${\mathit{s}}_j$ be the center of pixel j for pixels $j=1,2,\dots ,N_p$. The value of the blurred sky map at pixel j is computed as $\theta ({\mathit{s}}_j)=\mathit{K}({\mathit{s}}_j)\mathbf{\Theta }$. Let K be an $N_p\times N_p$ matrix where the jth row of K is $\mathit{K}({\mathit{s}}_j)$. We define $\mathit{\theta }$ to be the blurred sky map, where $\mathit{\theta }=\mathit{K}\mathbf{\Theta }$.

The way we estimate $\mathit{\theta }$ is by fitting a surface to the N_o locations ${\mathit{s}}_i$ corresponding to the binned direct event data. With a fitted surface, we can predict the surface at any location s, including the N_p pixel center locations that define $\mathit{\theta }$.

Recall from (1) that $$\mathrm{E}(y({\mathit{s}}_i) | t({\mathit{s}}_i),\theta ({\mathit{s}}_i),b({\mathit{s}}_i))=t({\mathit{s}}_i)[\theta ({\mathit{s}}_i)+b({\mathit{s}}_i)].$$

Let $$z({\mathit{s}}_i)=y({\mathit{s}}_i)/t({\mathit{s}}_i)-b({\mathit{s}}_i)$$be the method of moments estimate of ENAs per second at location ${\mathit{s}}_i$, treated as an initial estimate of $\theta ({\mathit{s}}_i)$.

For Theseus Stage 1, let ${\mathit{s}}_i$ be the ecliptic longitude/latitude location in spherical coordinates. Theseus Stage 1 uses spherical coordinates rather than ecliptic coordinates for estimation to naturally account for spherical wrapping and prevent edge effects. We model $z({\mathit{s}}_i)$ as a noisy realization from $\theta ({\mathit{s}}_i)$: $$z({\mathit{s}}_i)=\theta ({\mathit{s}}_i)+ϵ({\mathit{s}}_i),$$where $ϵ({\mathit{s}}_i)$ is a mean 0 error term. Our goal is to estimate the spatial surface, $\mathit{\theta }$. We estimate $\mathit{\theta }$ in three steps:

Step 1: Estimate candidate blurred sky map ${\stackrel{.}{\mathit{\theta }}}_l$ for $l=1,2,\dots ,N_L$, where N_L is the number of candidate blurred sky maps.

Step 2: Estimate a single, initial blurred sky map $\stackrel{.}{\mathit{\theta }}$ as an ensemble of the candidate blurred sky maps ${\stackrel{.}{\mathit{\theta }}}_l$.

Step 3: Estimate the residual-adjusted blurred sky map $\stackrel{..}{\mathit{\theta }}$ as a residual-adjusted version of $\stackrel{.}{\mathit{\theta }}$.

These three steps provide the high-level roadmap for Theseus Stage 1. There are, however, many subjective decisions that must be made in order to operationalize these steps. It is difficult to make “optimal” decisions in a discovery science setting where minimal assumptions about the underlying structure of the sky map are available. Therefore, in this article we focused on developing procedures that protect against sub-optimal decisions rather than seeking optimal ones. We will show that the collective set of decisions we made ultimately led to improved sky map estimates.

Step 1: Estimate ${\stackrel{.}{\mathit{\theta }}}_l$. We anticipate the blurred sky maps to vary continuously over space. The goal of this step is to estimate a diverse and computationally inexpensive collection of candidate blurred sky maps using readily available software. We consider two different regression approaches to estimate $\mathit{\theta }$: projection pursuit regression (PPR) (Friedman and Stuetzle 1981) and generalized additive models (GAMs) (Hastie and Tibshirani 1987). PPR and GAMs are both general purpose, nonparametric regression models, useful for estimating continuous surfaces of varying degrees of smoothness. They do so by estimating $\mathit{\theta }$ as the sum of smooth functions of inputs (e.g., spherical locations ${\mathit{s}}_i$). Numerous choices are required before fitting can commence, including the choice of smoothing functions and observation weighting, where the observation weighting addresses the non-constant variance of $z({\mathit{s}}_i)$. The details of the PPR and GAM fits can be found in the supplementary materials.

To assist in describing Theseus, we will present a working example in Sections 3.2 and 3.3 corresponding to a single bootstrapped dataset of a simulated binned direct event dataset generated from a simulated Spatial Retention sky map (details in Section 4). Figure 2 shows N_L = 8 candidate blurred sky map estimates: 4 PPR-generated and 4 GAM-generated. At a high-level, each candidate estimates a similar blurred sky map: they all estimate the highest ENA rate around ecliptic longitude/latitude (265, –15), and two other lobes of elevated ENA rates near (100, 50) and (90, –45). The specifics of the estimated blurred sky maps, however, noticeably vary from candidate to candidate. For instance, Candidates 5 through 8 (the GAM-generated candidates) tend to be smoother but do not obviously estimate a ribbon, while Candidates 2 and 4 (the PPR-generated candidates using exposure time weighting) do visually estimate a ribbon. Banding artifacts are present in some of the PPR-generated candidate sky maps (e.g., Candidates 1 and 2) as a result of the internal PPR rotations and ridge function estimates. As Figure 2 makes clear, the particular choices of smoothing functions, weightings, and regression approaches are consequential and challenging to select a priori. Rather than try to determine what the “optimal” smoothing function/weighting/regression approach is, we combine the candidate blurred sky maps to form a single, initial blurred sky map estimate.

Fig. 2 The N_L = 8 candidate blurred sky map estimates and the simulated binned direct event data, $z({\mathit{s}}_i)$. The gray points in the simulated data panel correspond to $z({\mathit{s}}_i)$ outside the color bar range.

Step 2: Estimate $\stackrel{.}{\mathit{\theta }}$. The goal of Step 2 is to combine the N_L candidate blurred sky maps ${\stackrel{.}{\mathit{\theta }}}_l$ into a single, initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}$. A straightforward way to do this would be to learn a weight $w_l\in [0,1]$ for each candidate map such that ${\sum }_{l=1}^{N_L}{w}_l=1$ and the resulting initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}={\sum }_{l=1}^{N_L}{w}_l{\stackrel{.}{\mathit{\theta }}}_l$ minimizes some distance between the binned direct event data and $\stackrel{.}{\mathit{\theta }}$. This approach, however, is restricted to assign a single weight to each candidate blurred sky map rather than assign different weights to different parts of candidate blurred sky maps. From Figure 2, it appears different candidate blurred sky maps may capture some parts of the blurred sky map better than others. For instance, the PPR-generated candidates might be better suited to capture the ribbon, while the GAM-generated candidates might be better suited to the GDF.

Instead of weighting each candidate sky map, we choose to combine the candidate blurred sky maps into a single, initial blurred sky map estimate by treating the candidate blurred sky maps as covariates in a second GAM regression model. Specifically, we fit the model (3) $$\begin{array}{c}y({\mathit{s}}_i) | t({\mathit{s}}_i),{\dot{\theta }}_1({\mathit{s}}_i),\dots ,{\dot{\theta }}_{N_L}({\mathit{s}}_i),b({\mathit{s}}_i)\\ \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\sim \text{Poisson}(t({\mathit{s}}_i)[\theta ({\mathit{s}}_i)+b({\mathit{s}}_i)])\end{array}$$(3) (4) $$\theta ({\mathit{s}}_i)={\alpha }_0+\sum _{l=1}^{N_L}{f}_l({\dot{\theta }}_l({\mathit{s}}_i))$$(4) (5) $$f_l({\dot{\theta }}_l({\mathit{s}}_i))=\sum _{r=1}^{N_{lr}}{\alpha }_{lr}{g}_{lr}({\dot{\theta }}_l({\mathit{s}}_i)).$$(5)

Cubic regression spline basis functions were used for $g_{lr}()$. We use functions within the mgcv package (Wood, Pya, and Säfken 2016) in R (R Core Team 2022) to carry out the estimation. Figure 3 shows the contribution each candidate blurred sky map makes to the initial blurred sky map estimate, $\stackrel{.}{\mathit{\theta }}$. We see that the initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}$ largely resembles the ribbon of Candidates 2 and 4 and the smooth GDF of Candidates 5 through 8. By combining the candidate blurred sky maps via the GAM-regression, the banding artifacts over the GDF from PPR-generated Candidates 1 and 2 (Figure 2) are greatly diminished.

Fig. 3 The candidate blurred sky map components $f_l({\dot{\theta }}_l({\mathit{s}}_i))$ (left) and the initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}$ (right). Following (4), the initial blurred sky map estimate is calculated as the sum of the candidate blurred sky map components $f_l({\dot{\theta }}_l({\mathit{s}}_i))$ plus the estimated constant ${\alpha }_0\approx 0.15$ ENAs/sec.

Step 3: Estimate $\stackrel{..}{\mathit{\theta }}$. The final step is to correct any systematic lack of fit that may remain between the initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}$ and the observed binned direct event data y. This step is included as a layer of protection. If $\stackrel{.}{\mathit{\theta }}$ is consistent with y (as is desired), then the residual-adjusted blurred sky map $\stackrel{..}{\mathit{\theta }}$—which is an estimate of $\mathit{\theta }$, but not Theseus’s final estimate of $\mathit{\theta }$—will be similar to $\stackrel{.}{\mathit{\theta }}$. If, however, $\stackrel{.}{\mathit{\theta }}$ is inconsistent with y, that inconsistency can be remedied.

Computing $\stackrel{..}{\mathit{\theta }}$ proceeds as follows. First, we compute the residuals for each observation: (6) $$\dot{r}({\mathit{s}}_i)=y({\mathit{s}}_i)-\dot{\lambda }({\mathit{s}}_i)$$(6) (7) $$\dot{\lambda }({\mathit{s}}_i)=t({\mathit{s}}_i)[\dot{\theta }({\mathit{s}}_i)+b({\mathit{s}}_i)].$$(7)

We then fit a GAM regression model to $\dot{r}({\mathit{s}}_i)$ as a function of $\dot{\theta }({\mathit{s}}_i)$ and ${\mathit{s}}_i$, again using cubic regression spline basis functions, resulting in GAM-fitted residuals $\widehat{r}({\mathit{s}}_i)$. The final step is to fit another GAM model like in Step 2 ((3)–(5)), replacing ${\dot{\theta }}_l({\mathit{s}}_i)$ in (4) and (5) with $\dot{\theta }({\mathit{s}}_i)$ and $\widehat{r}({\mathit{s}}_i)$.

The right panel of Figure 4 shows the residual-adjusted blurred sky map estimate, $\stackrel{..}{\mathit{\theta }}$, as the sum of the initial blurred sky map estimate component (left) and the fitted residual surface component (middle). The right panel of Figure 4 shows the residual-adjusted blurred sky map estimate $\stackrel{..}{\mathit{\theta }}$. The residual adjustment was quite minor, as the initial blurred sky map estimate $\stackrel{.}{\mathit{\theta }}$ was already consistent with the binned direct event data y. In general, it is overwhelmingly common but not guaranteed that $\stackrel{.}{\mathit{\theta }}\approx \stackrel{..}{\mathit{\theta }}$.

Fig. 4 The initial blurred sky map component (left), the residual-adjustment component (middle) and the residual-adjusted blurred sky map estimate $\stackrel{..}{\mathit{\theta }}$ (right). $\stackrel{..}{\mathit{\theta }}$ is computed by adding the left and middle panels. Relatively minor residual adjustments were made to the initial blurred sky map component, as noted by the scale of the residual-adjustment component compared to the scale of the initial blurred sky map component.

3.3 Theseus Stage 2: Estimate an Unblurred Sky Map $\mathbf{\Theta }$

Recall that the $N_p\times 1$ blurred sky map $\mathit{\theta }$ is related to the $N_p\times 1$ unblurred sky map $\mathbf{\Theta }$ through the $N_p\times N_p$ PSF matrix K, where $\mathit{\theta }=\mathit{K}\mathbf{\Theta }$ and $\mathit{K}({\mathit{s}}_j)$—the jth row of K—is the blurring operator applied to the jth pixel center ${\mathit{s}}_j$. For this application, ${\mathit{K}}^{-1}$ exists, suggesting that the unblurred sky map $\mathbf{\Theta }$ can be easily computed given the blurred sky map $\mathit{\theta }$ as follows: $\mathbf{\Theta }={\mathit{K}}^{-1}\mathit{\theta }$. While this relationship holds for the true unblurred and true blurred sky maps, there is a problem when $\mathit{\theta }$ is replaced with an estimate $\widehat{\mathit{\theta }}$.⁷ This problem is related to the large condition number of the PSF matrix K (over 10 million), thus, the inverse problem is ill-conditioned (Rice 1966). A large condition number for the matrix K indicates that small differences between an estimated blurred sky map $\widehat{\mathit{\theta }}$ and the true blurred sky map $\mathit{\theta }$ can result in large differences between their inverses (i.e., between $\widehat{\mathbf{\Theta }}={\mathit{K}}^{-1}\widehat{\mathit{\theta }}$ and $\mathbf{\Theta }={\mathit{K}}^{-1}\mathit{\theta }$).

This phenomenon is illustrated in Figure 5. Four similar but not identical blurred sky maps are shown on the right and their unblurred counterparts are shown on the left. For each pair of unblurred/blurred sky maps, it is true that the K times the unblurred sky map equals the blurred map, and ${\mathit{K}}^{-1}$ times the blurred sky map equals the unblurred sky map. Figure 5 makes clear, however, that even if an estimated blurred sky map $\widehat{\mathit{\theta }}$ is nearly identical to the true blurred sky map $\mathit{\theta }$, computing the unblurred sky map as $\widehat{\mathbf{\Theta }}={\mathit{K}}^{-1}\widehat{\mathit{\theta }}$ is unreliable. This immediately raises challenges for unblurring $\stackrel{..}{\mathit{\theta }}$, the Theseus Stage 1 estimate of $\mathit{\theta }$.

Fig. 5 A diagram illustrating the blurring/unblurring (convolution/deconvolution) process. Unblurred sky map space is on the left; blurred sky map space is on the right. The true unblurred/blurred sky map is outlined in red, while the estimated pairs of unblurred/blurred sky maps are in black with hats. For all unblurred/blurred sky map pairs (true or estimates), it is true that $\mathit{\theta }=\mathit{K}\mathbf{\Theta }$ and ${\mathit{K}}^{-1}\mathit{\theta }=\mathbf{\Theta }$. The point spread function matrix K has a large condition number of over ten million. A large condition number poses a challenge for deconvolution, as small differences between to $\widehat{\mathit{\theta }}$ and $\mathit{\theta }$ can result in large differences between $\widehat{\mathbf{\Theta }}$ and $\mathbf{\Theta }$.

To address the challenges with unblurring (also referred to as “unfolding” (Lyons 2011)), we turn to regularization (e.g., Kuusela and Panaretos 2015; Kuusela and Stark 2017). Regularization requires making some additional assumption about the unblurred sky map. This can be challenging in general, but particularly so in our discovery science setting where we are willing to make few assumptions about the sky maps. What we do assume about the unblurred sky map is that it is “close” to the blurred sky map. That is, the blurred sky map is just a less crisp version of the unblurred sky map as is the case with ${\widehat{\mathbf{\Theta }}}_1$ and $\mathbf{\Theta }$ in Figure 5. By making this assumption, we could estimate $\mathbf{\Theta }$ by minimizing the following equation: (8) $$\widehat{\mathbf{\Theta }}=\underset{\mathbf{\Theta }}{\arg \min }\left|\right|\stackrel{..}{\mathit{\theta }}-\mathit{K}\mathbf{\Theta }|{|}_2^2+\varphi ||\stackrel{..}{\mathit{\theta }}-\mathbf{\Theta }|{|}_2^2,$$(8) where $\left|\right|\stackrel{..}{\mathit{\theta }}-\mathit{K}\mathbf{\Theta }|{|}_2^2$ encourages the convolution of the unblurred sky map $\mathbf{\Theta }$ with the PSF matrix K to be close to a blurred sky map estimate $\stackrel{..}{\mathit{\theta }}$, while $\varphi \left|\right|\stackrel{..}{\mathit{\theta }}-\mathbf{\Theta }|{|}_2^2$ encourages the difference between the unblurred sky map and a blurred sky map estimate to be small. Equation (8)(8) $$\widehat{\mathbf{\Theta }}=\underset{\mathbf{\Theta }}{\arg \min }\left|\right|\stackrel{..}{\mathit{\theta }}-\mathit{K}\mathbf{\Theta }|{|}_2^2+\varphi ||\stackrel{..}{\mathit{\theta }}-\mathbf{\Theta }|{|}_2^2,$$(8) corresponds to a generalized ridge regression problem (van Wieringen 2015, sec. 3.4).

The relationship between the blurred and unblurred sky map can be rewritten. Let $\mathbf{\Theta }=\mathit{\theta }+\mathit{\delta }$, where $\mathit{\delta }$ is an $N_p\times 1$ vector called the sharpening map that equals $\mathbf{\Theta }-\mathit{\theta }$. It is called the sharpening map because it plays the role of transforming the blurred sky map $\mathit{\theta }$ into the unblurred sky map $\mathbf{\Theta }$ by sharpening the blurred sky map. It follows that (9) $$\mathit{\theta }=\mathit{K}\mathbf{\Theta }=\mathit{K}(\mathit{\theta }+\mathit{\delta })\Rightarrow \mathit{\psi }\equiv \mathit{\theta }-\mathit{K\theta }=\mathit{K\delta }.$$(9)

The vector $\mathit{\psi }$ is the once blurred sky map $\mathit{\theta }$ minus the twice blurred sky map $\mathit{K\theta }$ and is estimated as $\widehat{\mathit{\psi }}=\stackrel{..}{\mathit{\theta }}-\mathit{K}\stackrel{..}{\mathit{\theta }}$. The PSF matrix K plays the role of dampening and smearing out the sharp, peaked features in the blurred sky map (e.g., the ribbon) but leaves the more gradually varying parts of the blurred sky map relatively unchanged (e.g., the GDF). Consequently, the outline of the ribbon can be seen by plotting $\widehat{\mathit{\psi }}$, as is done in the left panel of Figure 6.

Fig. 6 (left) The difference between the once and twice blurred residual-adjusted blurred sky map estimate: $\widehat{\mathit{\psi }}=\stackrel{..}{\mathit{\theta }}-\mathit{K}\stackrel{..}{\mathit{\theta }}$. The ribbon is present as the largest deviations from 0. (middle) $\widehat{\mathit{\psi }}$ versus $\mathit{K}\stackrel{.}{\mathit{\delta }}$ and (right) $\widehat{\mathit{\psi }}$ versus $\mathit{K}\widehat{\mathit{\delta }}$. The red line is the y = x line. There is a misfit between $\mathit{K}\stackrel{.}{\mathit{\delta }}$ and $\widehat{\mathit{\psi }}$, where $\mathit{K}\stackrel{.}{\mathit{\delta }}$ underestimates large, positive values of $\widehat{\mathit{\psi }}$ and overestimates large, negative values of $\widehat{\mathit{\psi }}$. This misfit is a result of the regularization imposed by ridge regression. $\widehat{\mathit{\delta }}$ corrects this misfit, resulting in a $\mathit{K}\widehat{\mathit{\delta }}$ that is better matched with $\widehat{\mathit{\psi }}$.

Under this new parameterization where $\mathbf{\Theta }=\mathit{\theta }+\mathit{\delta }$, assuming $\left|\right|\mathit{\theta }-\mathbf{\Theta }|{|}_2^2$ is small is equivalent to assuming that $\left|\right|\mathit{\delta }|{|}_2^2$ is small. Thus, (8) can be rewritten as a standard ridge regression problem (Hoerl and Kennard 1970). Specifically, $\mathbf{\Theta }$ could be estimated as (10) $$\stackrel{.}{\mathbf{\Theta }}=\stackrel{..}{\mathit{\theta }}+\stackrel{.}{\mathit{\delta }}$$(10) where (11) $$\stackrel{.}{\mathit{\delta }}=\underset{\mathit{\delta }}{\arg \min }\left|\right|\widehat{\mathit{\psi }}-\mathit{K\delta }|{|}_2^2+\varphi |\left|\mathit{\delta }\right|{|}_2^2.$$(11)

Equation (10)(10) $$\stackrel{.}{\mathbf{\Theta }}=\stackrel{..}{\mathit{\theta }}+\stackrel{.}{\mathit{\delta }}$$(10) says that the unblurred sky map can be estimated as the sum of the residual-adjusted blurred sky map estimate $\stackrel{..}{\mathit{\theta }}$ and an initial sharpening map estimate $\stackrel{.}{\mathit{\delta }}$. The initial sharpening map estimate $\stackrel{.}{\mathit{\delta }}$ is the solution to the ridge regression problem. The regularization part of (11), namely $\varphi \left|\right|\mathit{\delta }|{|}_2^2$, encourages $\mathit{\delta }$ to be small (i.e., encourages the unblurred and blurred sky maps to be “close”). The tuning parameter $\varphi \ge 0$ controlling the degree of regularization in (11) is learned via cross-validation. We use the glmnet package (Friedman, Hastie, and Tibshirani 2010) in R to perform the ridge regression, imposing constraints such that $\stackrel{..}{\mathit{\theta }}+\stackrel{.}{\mathit{\delta }}\ge \mathit{0}$.⁸ To the best of our knowledge, the rearranging of terms resulting in the sharpening map in (9) and estimating the sharpening map via ridge regression constitutes a novel approach to image deblurring.

Estimating $\mathit{\delta }$ via ridge regression, however, creates a misfit between $\widehat{\mathit{\psi }}$ and $\mathit{K}\stackrel{.}{\mathit{\delta }}$. From (9), we have that $\mathit{\psi }=\mathit{K\delta }$. The middle panel of Figure 6 plots $\widehat{\mathit{\psi }}$ versus $\mathit{K}\stackrel{.}{\mathit{\delta }}$ where a clear misfit can be seen. $\mathit{K}\stackrel{.}{\mathit{\delta }}$ underestimates $\widehat{\mathit{\psi }}$ for large, positive values and overestimates $\widehat{\mathit{\psi }}$ for large, negative values. The implication of using $\stackrel{.}{\mathit{\delta }}$ in (10) to estimate the unblurred sky map $\mathbf{\Theta }$ is that the base and peak of the ribbon may be over- and underestimated, respectively.

To correct this misfit, we estimate the sharpening map as a misfit-adjusted version of $\stackrel{.}{\mathit{\delta }}$. We assume the sharpening map $\mathit{\delta }$ can be written as (12) $$\mathit{\delta }=\mathit{A\alpha }$$(12) where ${\mathit{\alpha }}_{3\times 1}=({\alpha }_0,{\alpha }_1,{\alpha }_2)\prime ,{\mathit{A}}_{N_p\times 3}=[\mathit{1},{\stackrel{.}{\mathit{\delta }}}_{-},{\stackrel{.}{\mathit{\delta }}}_{+}]$, 1 is a vector of 1s, and the vectors ${\stackrel{.}{\mathit{\delta }}}_{-}$ and ${\stackrel{.}{\mathit{\delta }}}_{+}$ are hinge function representations of the vector $\stackrel{.}{\mathit{\delta }}$, where all entries of $\stackrel{.}{\mathit{\delta }}$ greater than 0 (${\stackrel{.}{\mathit{\delta }}}_{-}$) or less than 0 (${\stackrel{.}{\mathit{\delta }}}_{+}$) are set to 0. Equation (12)(12) $$\mathit{\delta }=\mathit{A\alpha }$$(12) just says that the relationship between $\mathit{\delta }$ and $\stackrel{.}{\mathit{\delta }}$ is linear, but a different linear relationship is possible between $\mathit{\delta }$ and the negative part of $\stackrel{.}{\mathit{\delta }}$ (${\stackrel{.}{\mathit{\delta }}}_{-}$) than with the positive part of $\stackrel{.}{\mathit{\delta }}$ (${\stackrel{.}{\mathit{\delta }}}_{+}$). The different linear relationships, however, are required to match when an entry of $\stackrel{.}{\mathit{\delta }}$ equals 0, preventing a jump discontinuity at 0 in (12). We estimate $\mathit{\alpha }$ and $\mathit{\delta }$ as follows: $$\begin{array}{c}\widehat{\mathit{\alpha }}=\underset{\mathit{\alpha }}{\arg \min }\left|\right|\widehat{\mathit{\psi }}-\mathit{KA\alpha }|{|}_2^2\\ \widehat{\mathit{\delta }}=\mathit{A}\widehat{\mathit{\alpha }}.\end{array}$$

The right panel of Figure 6 plots $\widehat{\mathit{\psi }}$ versus $\mathit{K}\widehat{\mathit{\delta }}$. As can be seen, the misfit that was present when using $\stackrel{.}{\mathit{\delta }}$ has been corrected when using $\widehat{\mathit{\delta }}$.

Figure 7 shows the final unblurred and blurred sky map estimates, where the final unblurred sky map estimate is $$\widehat{\mathbf{\Theta }}=\stackrel{..}{\mathit{\theta }}+\widehat{\mathit{\delta }}$$and the final blurred sky map estimate is $$\widehat{\mathit{\theta }}=\mathit{K}\widehat{\mathbf{\Theta }}.$$

Fig. 7 From left to right, the residual-adjusted blurred sky map estimate $\stackrel{..}{\mathit{\theta }}$, the sharpening map estimate $\widehat{\mathit{\delta }}$, the final unblurred sky map estimate $\widehat{\mathbf{\Theta }}=\stackrel{..}{\mathit{\theta }}+\widehat{\mathit{\delta }}$, and the final blurred sky map estimate $\widehat{\mathit{\theta }}=\mathit{K}\widehat{\mathbf{\Theta }}$.

The effect of adding the estimated sharpening map $\widehat{\mathit{\delta }}$ to the estimated blurred sky map $\stackrel{..}{\mathit{\theta }}$ is that the lower latitude side of the ribbon in $\stackrel{..}{\mathit{\theta }}$ (i.e., the blue streak in $\widehat{\mathit{\delta }}$) increases in ENA rate while the upper latitude side of the ribbon in $\stackrel{..}{\mathit{\theta }}$ (i.e., the red streak in $\widehat{\mathit{\delta }}$) decreases in ENA rate. As these ribbon increases and decreases to $\stackrel{..}{\mathit{\theta }}$ are in close spatial proximity, $\widehat{\mathit{\delta }}$ creates a more dramatic gradient from higher ENA rates to lower ENA rates on the upper latitude side of the ribbon.

Figure 7 (and all the illustrative figures in Section 3) correspond to one bootstrapped binned direct event dataset ${\mathit{D}}_b$. That is, for a given ESA step (e.g., ESA step 2), a given 6-month map (e.g., 2013A), and a given bootstrapped binned direct event dataset ${\mathit{D}}_b$, Theseus estimates an unblurred sky map ${\widehat{\mathbf{\Theta }}}_b$ and a blurred sky map ${\widehat{\mathit{\theta }}}_b$. The final point estimate unblurred sky map $\widehat{\mathbf{\Theta }}$ and point estimate blurred sky map $\widehat{\mathit{\theta }}$ for a given ESA and 6-month map are computed as averages over all N_B = 1000 bootstrap estimates: $$\widehat{\mathbf{\Theta }}=N_B^{-1}\sum _{b=1}^{N_B}{\widehat{\mathbf{\Theta }}}_b\text{and}\widehat{\mathit{\theta }}=N_B^{-1}\sum _{b=1}^{N_B}{\widehat{\mathit{\theta }}}_b.$$

Estimating $\widehat{\mathbf{\Theta }}$ as the average over all bootstrap samples (i.e., bagging (Breiman 1996)) rather than a single run of Theseus applied to the original binned direct event dataset D provides stability to the final sky map estimate and can greatly reduce the banding artifacts initially caused by the PPR regression in step 1 of Theseus Stage 1 that can remain in ${\widehat{\mathbf{\Theta }}}_b$ (Figure 7). This is because, though PPR very often has banding artifacts, those artifacts are not always found in the same location across bootstrap samples. By averaging over all N_B unblurred sky map estimates, the banding artifacts are often (but not always) blended out of the final unblurred sky map estimate $\widehat{\mathbf{\Theta }}$. Further strategies for addressing the banding artifacts are discussed in Section 6.

In Sections 4 and 5, each sky map is estimated with N_B = 1000 bootstrap samples. Bootstraps were run in parallel using R version 4.1.2 on a 2.6 GHz AMD EPYC 7513 processor with 2 32-core sockets having 2 threads per core (128 threads total). Each sky map took 10–15 min to compute. Code to reproduce aspects of this manuscript for one bootstrap sample is available at https://github.com/lanl/Theseus.

4 Simulation Study

To investigate the inferential properties of Theseus and compare Theseus to a near facsimile of the current state-of-the-art sky map estimation procedure used by the ISOC, we perform a simulation study. The ISOC sky map estimation method is a local estimation method, where each pixel depends only on “nearby” binned direct event data. The pixel resolution is fixed at 6 degrees. As such, ISOC maps have a coarser and noisier esthetic than do Theseus maps. A complete description of the implemented ISOC sky map estimation methodology is provided in the supplementary materials.

4.1 Simulation Study Setup

We consider four simulated sky maps with different ribbon models motivated by some of the most advanced theoretical ribbon models available in the literature:

• GDF-only: a sky map with no ribbon, just globally distributed flux (Zirnstein et al. 2017). Elevated ENA rates are present at the nose and tail.

• Weak Scattering: the GDF-only sky map plus a ribbon with heterogeneous brightening and a relatively smooth, symmetric ribbon profile (Zirnstein, Heerikhuisen, and Dayeh 2018).

• Spatial Retention: the GDF-only sky map plus a ribbon with heterogeneous brightening and an asymmetric ribbon profile (Schwadron and McComas 2013; Zirnstein et al. 2019).

• Varying Ribbon Profile: the GDF-only sky map plus a ribbon that transitions from the Weak Scattering ribbon to the Spatial Retention ribbon with highly-concentrated pockets of higher ENA rates added at locations on the ribbon as well as a disk over the nose that appears suddenly. This model is intentionally unrealistic, but is useful for testing the limits of Theseus and ISOC.

The different simulated unblurred sky maps are shown in the top row of Figure 8. Each model is simulated for ENAs at ∼1.7 keV, corresponding to IBEX-Hi ESA step 4. The Weak Scattering and Spatial Retention ribbon model results shown in Figure 8 are extracted from a single cross-sectional cut across the model ribbon and then copied across the sky to mimic a ribbon with consistent cross-section shape. The change in ribbon intensity across the sky is artificially varied to mimic reality. Binned direct event counts y were simulated from (1) and (2). Specifically, for a given unblurred sky map $\mathbf{\Theta }$ and the locations, exposure times, and background rates of ESA steps 2 through 6 of the 2018B map, we draw $y({\mathit{s}}_i)\sim \text{Poisson}(t({\mathit{s}}_i)[\theta ({\mathit{s}}_i)+b({\mathit{s}}_i)])$ where $\theta ({\mathit{s}}_i)=\mathit{K}({\mathit{s}}_i)\mathbf{\Theta }$ for $i=1,2,\dots ,N_o$. Five binned direct event datasets D were simulated for each ribbon model. This amounts to five different simulated binned direct event datasets for each of the four ribbon models, resulting in 20 total ESA/ribbon model combinations.

Fig. 8 The true and estimated sky maps from Theseus and ISOC (rows) for the different simulated ribbon models (columns) for ESA 6. Visually, there appears to be better agreement between Theseus and the true sky maps than ISOC and the true sky maps.

4.2 Visual Comparison

Figure 8 plots unblurred sky map estimates for the different ribbon models based on ESA 6 binned direct event data. Visually, there appears to be better agreement between Theseus and the true sky maps than ISOC and the true sky maps. ISOC sky maps appear grainier than Theseus because ISOC makes sky maps at a fixed $6^{°}$ resolution, while the resolution is a user-defined parameter for Theseus ($2^{°}$ maps are shown in Figure 8). ISOC sky map estimates are done on a pixel-by-pixel basis, and thus can suffer from erratic pixel estimates in regions with little to no binned direct event data. Theseus, on the other hand, borrows information over space via the PPR and GAM smoothing, resulting in less erratic pixel-to-pixel sky map estimates. The knots present in the Varying Ribbon Profile near ecliptic longitude and latitude (120, 40) and (300, 20) are not captured by either Theseus or ISOC. A similar comparison between Theseus and the simulated blurred sky maps can be found in the supplementary materials.

4.3 Agreement with Simulated Binned Direct Event Data

To assess how well Theseus-estimated blurred sky maps and the ISOC-estimated unblurred sky maps agree with the simulated binned direct event data, we examine the probability integral transform (PIT) histograms (Diebold, Gunther, and Tay 1998; Gneiting, Balabdaoui, and Raftery 2007). More concretely, let PIT($y({\mathit{s}}_i)$, $\widehat{\lambda }({\mathit{s}}_i)$) be a random draw from Uniform($l({\mathit{s}}_i),u({\mathit{s}}_i)$), where $l({\mathit{s}}_i)$ is the cumulative distribution function (cdf) of a Poisson evaluated at $y({\mathit{s}}_i)-1$ with mean parameter $\widehat{\lambda }({\mathit{s}}_i)=t({\mathit{s}}_i)[\widehat{\theta }({\mathit{s}}_i)+b({\mathit{s}}_i)]$ and $u({\mathit{s}}_i)$ is the cdf of a Poisson evaluated at $y({\mathit{s}}_i)$. Note that, if $y({\mathit{s}}_i)=0$, then $l({\mathit{s}}_i)=0$. The uniform draw is required because of the discreteness of the counts $\mathit{y}({\mathit{s}}_i)$ but is not required in the continuous setting. If the binned direct event observations were conditionally independent draws from their generating data model with the blurred sky map estimate $\widehat{\mathit{\theta }}$ plugged in, we would expect the collection of PIT($y({\mathit{s}}_i),\widehat{\lambda }({\mathit{s}}_i)$) to follow a Uniform(0,1) distribution.

In addition to PIT histograms, we assess the spatial distribution of the PIT values. If PIT($y({\mathit{s}}_i),\widehat{\lambda }({\mathit{s}}_i)$) $\stackrel{\text{iid}}{\sim }$ Uniform(0,1) and were randomly assigned to binned direct event locations ${\mathit{s}}_i$, we would expect the mean and variance of the difference between nearest neighbor PIT values to be 0 and 1/6 (the mean and variance of the difference of two independent, Uniform(0,1) random variables).

For each simulated map and ESA, PIT values and nearest neighbor differences of PIT values were bootstrapped 1000 times, recording their means and variances plotted in Figure 9. The results are largely consistent with expectations under theory, indicating no serious lack-of-fit for either ISOC or Theseus. That said, Theseus PIT histograms have better agreement with nominal values than do the ISOC PIT histograms, especially for the PIT variances. The ISOC PIT variances are systematically lower than nominal, suggesting ISOC maps are over-fitting the binned direct event data.

Fig. 9 Bootstrapped PIT histogram means (PIT Mean) and variances (PIT Var) and nearest neighbor PIT histogram difference means (PIT Diff Mean) and variances (PIT Diff Var) for the different simulated sky maps (columns) between the binned direct event data and the ISOC (dark green) and Theseus (mustard) estimated blurred sky maps. The nominal mean and variance for the PIT histograms are 1/2 and 1/12, respectively. The nominal mean and variance for differences of nearest neighbor PIT histograms are 0 and 1/6, respectively. Nominal values are denoted by red, horizontal lines. Results show no serious lack-of-fit for ISOC or Theseus, though Theseus has PIT histograms closer to the nominal means/variances than does ISOC.

4.4 Unblurred Estimated Sky Map Comparisons to Simulated Sky Maps

To assess the quality of Theseus and ISOC unblurred sky map estimates, we compute the mean absolute percent error (MAPE) to assess the accuracy of the point estimates, the continuous rank probability score (CRPS) to assess the skill (Matheson and Winkler 1976; Gneiting and Raftery 2007), the 95% confidence interval widths to asses the sharpness, and the 95% empirical coverage to assess the reliability. For details on how these metrics are calculated, see the supplementary materials.

Figure 10 presents the results comparing the unblurred sky map estimates to the true, unblurred sky maps. As shown in Figure 10, on average, ISOC MAPE is approximately 85% larger than Theseus MAPE, ISOC CRPS is approximately 80% larger than Theseus CRPS, and ISOC 95% confidence interval widths are approximately 65% wider than Theseus CI widths. Furthermore, Theseus has empirical coverage closer to nominal coverage than ISOC for nearly all ribbon models and ESAs. Theseus’s empirical coverages were within 2.5% of nominal for all ribbon models and ESAs except for ESAs 3, 4, and 6 of the Varying Ribbon Profile model. ISOC’s empirical coverages were within 2.5% of nominal for less than half of the 20 ESA/ribbon model combinations.

Fig. 10 The mean absolute percent error (MAPE), 95% confidence interval width (CI Width), 95% empirical coverage (Coverage), and continuous rank probability score (CRPS) for ISOC and Theseus for each simulated sky map (columns) and ESAs (x-axis). For all metrics other than Coverage, the best value is 0. For Coverage, the best value is 0.95 denoted by horizontal, red line. Theseus outperforms ISOC across all metrics and simulated sky map scenarios.

The unrealistic Varying Ribbon Profile ribbon model was the most challenging for ISOC and Theseus to estimate. For Theseus, this is likely related to the assumption made in Stage 2 of Theseus that the unblurred and blurred sky maps are “close” to each other: an assumption that is most significantly violated in the Varying Profile Ribbon model with the sharpest and most asymmetric ribbon shape.

The results presented in Figure 10 provide strong evidence that Theseus sky map estimates are better than ISOC sky map estimates. The same analysis was performed, partitioning the sky map into the pixels that contain the ribbon and those that do not with similar findings (see the supplementary materials for more details).

4.5 Comparing the Ribbon Shapes

The shape of the ribbon is of scientific interest. In this section, we isolate and compare the ribbon profiles from the Theseus and ISOC sky maps to the truth. Recall that the ribbon forms a nearly complete circular band across the sky. In order to isolate and compare the ribbon profiles, we first project the sky maps into a ribbon-centered frame (details available in the supplementary materials). In ribbon-centered frame, the ribbon appears as a horizontal streak (see the supplementary materials). We are able to do this ribbon-centering for simulated maps without error because we know the ribbon center location.⁹

Figure 11 plots the ribbon profiles for the ESA 6 sky maps (i.e., the ENA rates for the vertical cross-sections between the horizontal dashed lines in Figure 5 in the supplementary materials). It appears that Theseus better visually captures the shape of the ribbon than ISOC. Theseus captures the smoothly varying ribbon shapes from adjacent cross-sections, while ISOC ribbon shape estimates vary more erratically as a result of ISOC’s pixel-by-pixel estimation procedure.

Fig. 11 Ribbon cross-sections for each ribbon model (columns) and sky map estimation method (rows). Each line within a panel corresponds to a ribbon-centered longitude (i.e., a vertical slice in Figure 5 of the supplementary materials). We can see that Theseus better approximates the shape of the ribbon than does ISOC. Neither ISOC nor Theseus, however, completely captures the asymmetry of the ribbon for the Spatial Retention model or the Varying Ribbon Profile model.

Not only does Theseus appear to capture the ribbon shape better than ISOC, it does a better job quantitatively. Figure 12 shows the mean absolute error between the skew of the simulated ribbon cross-sections and the ISOC and Theseus-estimated skews (see the supplementary materials for details on how skew is computed). Figure 12 provides further evidence that Theseus better captures relevant features of the ribbon than does ISOC.

Fig. 12 Mean absolute error (MAE) between the true ribbon skew and the Theseus and ISOC estimated ribbon skew, averaged over all cross-sections. 1000 bootstrap samples are shown. Theseus has smaller average MAE than ISOC for all ribbon model/ESA comparisons indicating Theseus better captures the skew of the ribbon than ISOC.

While Theseus does a better job capturing the ribbon shape than does ISOC, both Theseus and ISOC fail to capture some of the finer features of the ribbon, such as the sharp peak/pronounced asymmetry of the Spatial Retention model and the strong asymmetric features of the Varying Ribbon Profile. ISOC’s failure to capture these sharp, asymmetric ribbon shapes is in part because it does not account for the IBEX-Hi PSF and is estimating a blurred sky map but treating it as an unblurred sky map. Theseus’s failure to capture these sharp, asymmetric ribbon shapes is in part related to the assumption made in Theseus Stage 2 that the unblurred sky map is “close” to the blurred sky map. While this assumption prevents unblurred sky map estimates from being wildly unreasonable (recall ${\widehat{\mathbf{\Theta }}}_2$ and ${\widehat{\mathbf{\Theta }}}_3$ in Figure 5), it also prevents Theseus from capturing the parts of the unblurred sky map that deviate most dramatically from the blurred sky map: namely the ribbon. To better capture the ribbon in the Spatial Retention model and Varying Ribbon Profile, additional assumptions would likely need to be made about the shape of the ribbon; assumptions intentionally avoided for the discovery science setting we are operating within.

5 Real Binned Direct Event Data Example

We fit Theseus and ISOC to the real binned direct event data collected by the IBEX-Hi instrument corresponding to the ESA 4 “A” maps between the years 2010 to 2021 to illustrate how the current state-of-the-art ISOC sky maps compare to Theseus sky maps. Unblurred sky maps are shown in Figure 13. Theseus sky maps appears noticeably smoother than ISOC sky maps. Theseus sky maps are at a $2^{°}$ pixel resolution and interpolate over regions of the sky while ISOC sky maps are made at a $6^{°}$ resolution on a pixel-by-pixel basis and leave “holes” in the sky maps where no binned direct event data are present.

Fig. 13 ISOC (left) and Theseus (right) unblurred sky map estimates fit to real binned direct event data corresponding to ESA 4 collected over roughly the first half each calendar year (the “A” maps). Gray pixels in the ISOC sky maps correspond to pixels ISOC does not estimate because no binned direct event data was collected sufficiently close to the pixel’s center.

While we do not know if ISOC or Theseus sky maps better reflects the true, latent unblurred sky maps, we have reason to believe it is Theseus due to the similar relative performance between Theseus and ISOC we see on real data as we did on simulated data (recall Section 4). For simulated data where Theseus sky maps better captured simulated sky maps than did ISOC, we saw that Theseus had narrower 95% confidence interval widths and PIT histograms better reflecting samples from a Uniform(0,1). We have the same findings with real data. For instance, ISOC’s 95% confidence interval widths are, on average, ∼50% wider than Theseus’s (Figure 14). Furthermore, Theseus PIT histograms are consistent with a Uniform(0,1) for all sky maps, while ISOC PIT histograms show a lower frequency for PIT values near 1 for all sky maps than would be expected, and a similar mounded shape for 2010A and 2011A that we saw with the simulated sky maps (Figure 15).

Fig. 14 ISOC (left) and Theseus (right) 95% confidence interval widths (ENAs per second) corresponding to ESA 4. Theseus has narrower 95% confidence intervals widths than ISOC overall. Gray pixels in the ISOC sky maps either represent pixels with no estimate or pixels with extremely large 95% confidence interval widths that were truncated for ease of viewing.

Fig. 15 Probability integral transform (PIT) histograms by years from 2010 to 2021 for ESA 4. The top set of histograms are for Theseus sky maps; the bottom for ISOC sky maps. No indication of non-uniformity is visible for Theseus PIT histograms. Slightly mounded shapes exists for 2010A and 2011A for ISOC maps. For later years, fewer PIT values near 1 are observed relative to what would be expected.

Artifacts of the Theseus and ISOC sky map estimation approaches are present. In the ISOC sky maps, we see coarse pixelation which can obscure finer scale features. We can also see some vertical banding in locations where binned direct event data exposure times are lower (e.g., in 2014A near ecliptic longitude 190 and ecliptic latitude range 0 to 90). If no binned direct event data are sufficiently close to a pixel’s center, ISOC will not provide an estimate for that pixel; those pixels appear as gray in Figure 13. For Theseus, some of the artificial banding caused by the PPR regression in step 1 of Theseus Stage 1 can be seen in the light blue bands running along the ribbon in the 2010A sky map and less so in 2011A sky map. While the artificial banding in Figure 13 is relatively minimal, it can be more pronounced, especially in ESA 2 where less binned direct event data is available to estimate the sky maps (not shown).

There are numerous scientific studies that will directly benefit from the demonstrated improvements found in Theseus sky maps over the standard ISOC sky maps. In terms of ribbon science, as already discussed there is the improved ability to discern the shape of the ribbon profile, and knowing this has the potential to differentiate between competing ribbon formation theories. Beyond this, the higher resolution of Theseus maps more accurately locates points along the ribbon in the sky. This will allow a more accurate application than previously done of the astronomical method of parallax to determine the distance to the ribbon (Swaczyna et al. 2016). More accurately demarcating the circular path of the ribbon around the sky will allow a better determination of the location of the ribbon center than previously done (Funsten et al. 2013; Dayeh et al. 2019). This is important because the ribbon center is associated with the direction of the interstellar magnetic field (e.g., Schwadron et al. 2009; Zirnstein et al. 2016c), which helps govern the shape of the heliosphere (Zirnstein et al. 2016a) and is also a topic of great interest in galactic astronomy (e.g., Frisch et al. 2012, 2022). And finally, the improved rendering of the GDF in Theseus sky maps will allow for more accurate determinations of the size and time variability of GDF features (e.g., McComas et al. 2013; Dayeh et al. 2022), which in turn helps constrain simulations of the heliosphere’s dynamics and evolution, and allows for a more accurate determination of the three-dimensional shape of the heliosphere than done previously (Reisenfeld et al. 2021).

6 Discussion

In this article, we presented a novel sky map estimation procedure called Theseus. Theseus produces sky map estimates in two stages. In Theseus Stage 1, a blurred sky map is estimated from binned direct event data. In Theseus Stage 2, an unblurred sky map is estimated from the previously estimated blurred sky map. Bootstrapping is used to quantify pixel uncertainties. Over a variety of realistic simulation scenarios, we demonstrated that Theseus is better than the current state-of-the-art sky map estimation procedure used by ISOC with respect to half a dozen metrics. Specifically, Theseus had more uniform PIT histograms, smaller MAPE, narrower 95% confidence interval widths, better empirical coverage, smaller CRPS, and more accurate ribbon skewness than ISOC sky maps. Finally, we fit Theseus and ISOC to real binned direct event data and showed that, like with the simulated data, Theseus had more uniform PIT histograms and narrower 95% confidence interval widths, providing some assurance that Theseus sky maps are likely better estimates of the true, latent unblurred sky maps than are ISOC’s.

While Theseus is a clear improvement over the current state-of-the-art used by the ISOC, it is by no means perfect. For instance, while there is no indication the blurred sky maps are over-smoothed based on the PIT histogram analysis of Section 4.3, visually there is a suggestion that Theseus blurred sky maps may not be capturing some fine scale spatial structure. Furthermore, though efforts to reduce banding artifacts in the sky map estimates were helpful, they did not completely remove them in all sky maps. The banding artifacts are caused by the PPR-generated candidates in step 1 of Theseus Stage 1 and can be amplified by the sharpening map of Theseus Stage 2. Even though they are relatively minor, the banding artifacts can be problematic for science as they might be confused as genuine heliospheric features rather than unfortunate methodological byproducts.

In the future, there are numerous options to try to reduce the impacts of these artifacts and address possible over-smoothing, both through methodological modifications and data enhancements. No single option is perfect, but all are worth exploring. Methodologically, we could replace the PPR-candidates in step 1 of Theseus Stage 1 with smoothed versions of themselves (e.g., we could replace ${\stackrel{.}{\mathit{\theta }}}_l$ with $g(\mathit{K}{\stackrel{.}{\mathit{\theta }}}_l)$; a debiasing function $g()$ applied to the smoothed/blurred version of ${\stackrel{.}{\mathit{\theta }}}_l$). This should reduce the banding artifacts in the candidate blurred sky maps but possibly at the expense of capturing the sharpness of the ribbon peaks. We could also consider alternatives (or additions) to PPR to include in step 1 of Theseus Stage 1 such as multivariate adaptive regression splines (Friedman 1991), neural networks (Rosenblatt 1958; Rumelhart, Hinton, and Williams 1986), local regression and likelihood methods (Loader 2006), or fast Vecchia approximations of Gaussian processes (Katzfuss and Guinness 2021; Katzfuss, Guinness, and Lawrence 2022). Providing a more diverse collection of potentially higher variance candidate blurred sky maps might address the perception of over-smoothing.

Alternatively, there is evidence to suggest that as the number of available ENA counts for binned direct event data increases, the banding artifacts fade (e.g., sky maps 2010A and 2011A have the lowest total number of counts of all sky maps in Figure 13). From a data enhancement perspective, we could incorporate a new data product. There is another data product related to the IBEX mission not used in this article that would effectively triple the event rate for the binned direct event data. Additionally, we could improve the data pre-processing step. Recall that the standard IBEX data pre-processing step discards roughly two-thirds of the direct event data when constructing binned direct event data. Improving the pre-processing step in a way that converts more direct event time segments into binned direct event data will increase total exposure time and may help reduce the banding artifacts.

Kuusela and Panaretos (2015) and Kuusela and Stark (2017) discuss how uncertainty estimation for regularized problems can go poorly if the regularization bias is not accounted for, especially in contexts where the unblurred and blurred functions/surfaces can be quite different from each other. It is worth emphasizing how the heliosphere sky map estimation problem avoids some of these more general concerns. First, the debiasing step in Theseus Stage 2 attempts to remove some of the misfit (loosely bias) introduced in the regularization step. Second, and most importantly, the degree of bias introduced by the point spread function is relatively small, which is why we are able to assume the unblurred and blurred sky maps are “close” to each other. Third, the empirical coverage estimates in the simulation study with realistic unblurred sky maps gives us confidence that our uncertainty estimates are reasonable.

One of the most challenging and scientifically important features for Theseus (and ISOC) to capture is the ribbon shape. Recall that the ribbon was completely unexpected and only discovered after IBEX’s launch. The IBEX-Hi’s ${6.5}^{°}$ FWHM PSF was anticipated to be, and is narrow enough to resolve GDF features. Resolving ribbon features, however, has proven more challenging. Assessing if Theseus sky maps are “good enough” to discern between competing ribbon origin theories will depend on the theories and how similar their resulting ribbon profiles are. If, for instance, competing theories can be differentiated by their skew (as is the case with the Weak Scattering and Spatial Retention models: recall Figure 11), then Theseus sky maps should be good enough to contribute to the scientific process now. If, however, two competing theories are quite similar and only deviate in, say, how peaked their ribbons are (imagine two slightly different versions of the Spatial Retention model), Theseus sky maps are likely not presently good enough for that level of discernment.

The challenge of improving ribbon shape estimation is further amplified in our discovery science setting, where few assumptions can be made a priori about the shape of the ribbon. Were we willing and able to make additional assumptions about the shape of the ribbon, however, we would likely be able to unblur the blurred sky maps more accurately. For instance, Kuusela and Stark (2017) provides an illustration of how unblurring can be improved when shape constraints can be assumed. Alternatively, as heliosphere simulation capabilities continue to mature, efforts to link computer model output directly to binned direct event data via computer model calibration should be explored.

In the absence of additional ribbon shape assumptions, improved instrumentation with narrower PSFs may be needed to accurately resolve the ribbon. Lucky for us, the upcoming Interstellar Mapping and Acceleration Probe (IMAP) (McComas et al. 2018) mission will have instrumentation with a FWHM field-of-view roughly half that of IBEX-Hi. Theseus is well-positioned to analyze IMAP data upon arrival after its scheduled 2025 launch.

Supplementary Materials

In the online supplementary materials of this article, we provide details on how the point spread function matrix K is constructed, additional GAM and PPR details, details on how ISOC maps are constructed, additional analyses supplementing Section 4, and a detailed table of the notation presented in the article. Additionally, the zip file contains R code and datasets to recreate Figures 2, 3, 4, 6, and 7. The R code and datasets can also be found at https://github.com/lanl/Theseus.

Acknowledgments

We thank the editor, associate editor, and two anonymous reviewers for their constructive and insightful feedback. We also thank C.C. Essix for her assistance and encouragement. Approved for public release: LA-UR-22-30879.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

Research presented in this manuscript was supported by the Laboratory Directed Research and Development (LDRD) program of Los Alamos National Laboratory (LANL) under project number 20220107DR and by the NASA IBEX Mission as part of the NASA Explorer Program (80NSSC20K0719). Additionally, L.J.B. was supported by LANL’s LDRD Richard Feynman Postdoctoral Fellowship (20210761PRD1), while E.J.Z. acknowledges support from NASA grant 80NSSC17K0597 in the development of the ribbon models.

Notes

1 For context, the outer boundary of the heliosphere, called the heliopause, is over 100 astronomical units from the Sun, where 1 astronomical unit is the mean distance between Earth and the Sun.

2 For the first few years of the IBEX mission, IBEX reorientation occurred approximately every 7 days, or once per orbit, not 4.5 days. In June 2011, the IBEX orbit was changed to a 9-day period orbit and since then, repoints are done twice per orbit.

3 Heliospheric ENA rates vary on the order of months and are treated as constant for each 6-month sky map.

4 Due to data quality issues, ESA 1—the lowest energy step—is not used.

5 180 sec × 6 ESA steps × 360 $1^{°}$ bins $=$ 388,800 sec $=$ 4.5 days, which is the duration of an arc.

6 For a $2^{°}$ sky map, $N_p=(180/2)*(360/2)=16,200$ pixels.

7 Recall the notational conventions of Section 3 and note here that the “hat” notation of $\widehat{\mathit{\theta }}$ is being used as a generic estimate of $\mathit{\theta }$, not as Theseus’s final estimate of $\mathit{\theta }$.

8 Further regularization was imposed on pixels proportional to their size. While all pixels are of the same size in degrees, they have different surface area sizes on a sphere as a function of latitude. Pixels near latitude 0 are the largest while pixels near latitudes –90 and 90 are the smallest. Smaller pixels were more severely regularized by using the penalty.factor argument in the glmnet() function. Without this extra regularization, the pixel estimate in $\dot{\mathit{\delta }}$ near latitudes –90 and 90 could wander appreciably from 0.

9 It is also of scientific interest, but outside the scope of this work, to estimate the ribbon center for non-simulated settings.