SkyFACT: High-dimensional modeling of gamma-ray emission with adaptive templates and penalized likelihoods

TL;DR

SkyFACT introduces a novel high-dimensional modeling approach for gamma-ray emission analysis, combining adaptive templates with penalized likelihoods to account for model imperfections and improve decomposition accuracy.

Contribution

The paper presents a new hybrid method that integrates image reconstruction and adaptive template regression with regularization, enabling efficient handling of large nuisance parameter spaces in gamma-ray modeling.

Findings

Successfully applied to synthetic and real Fermi-LAT data

Effectively models complex diffuse gamma-ray emission components

Provides a flexible framework for future diffuse emission studies

Abstract

We present SkyFACT (Sky Factorization with Adaptive Constrained Templates), a new approach for studying, modeling and decomposing diffuse gamma-ray emission. Like most previous analyses, the approach relies on predictions from cosmic-ray propagation codes like GALPROP and DRAGON. However, in contrast to previous approaches, we account for the fact that models are not perfect and allow for a very large number ($\gtrsim10^5$) of nuisance parameters to parameterize these imperfections. We combine methods of image reconstruction and adaptive spatio-spectral template regression in one coherent hybrid approach. To this end, we use penalized Poisson likelihood regression, with regularization functions that are motivated by the maximum entropy method. We introduce methods to efficiently handle the high dimensionality of the convex optimization problem as well as the associated semi-sparse covariance matrix, using the L-BFGS-B algorithm and Cholesky factorization. We test the method both on synthetic data as well as on gamma-ray emission from the inner Galaxy, $|\ell|<90^\circ$ and $|b|<20^\circ$, as observed by the Fermi Large Area Telescope. We finally define a simple reference model that removes most of the residual emission from the inner Galaxy, based on conventional diffuse emission components as well as components for the Fermi bubbles, the Fermi Galactic center excess, and extended sources along the Galactic disk. Variants of this reference model can serve as basis for future studies of diffuse emission in and outside the Galactic disk.

Tables

Table 1. Table 1 : Summary of components and associated regularization hyper-parameters used for the models Run1 to Run5 in this work. For each component, we show in the compact matrix notation the adopted regularization hyper-parameters for the spatial ( λ 𝜆 \lambda ), spectral ( λ ′ superscript 𝜆 ′ \lambda^{\prime} ) and overall ( λ ′′ superscript 𝜆 ′′ \lambda^{\prime\prime} ) modulation parameters, and for the spatial ( η 𝜂 \eta ) and spectral ( η ′ superscript 𝜂 ′ \eta^{\prime} ) smoothing parameters. The definition of hyper-parameters is shown in Eq. ( 2.7 ). For the spatial and spectral hyper-parameters, an ‘ ∞ \infty ’ indicates that the corresponding modulation parameters are kept fixed to one, while for the smoothing parameters, a ‘0’ indicates no smoothing on that component. We also show the total number of fit parameters N param subscript 𝑁 param N_{\text{param}} , the naive number of DOF ( N ebin × N pix − N par subscript 𝑁 ebin subscript 𝑁 pix subscript 𝑁 par N_{\text{ebin}}\times N_{\text{pix}}-N_{\text{par}} ), the effective number of DOF as estimated from mock data along with the effective number of model parameters and data bins, the values of the Poisson and regularization parts of the likelihood function, and the model fidelity, ℱ ℱ \mathcal{F} as an indicator for the goodness-of-fit (see section  6.1 ). Effective fit parameters are rounded to the nearest hundred.

Components		Run1	Run2	Run3	Run4	Run5
		Regularization hyper-parameters:
IGRB
3FGL PSC
Gas (0–19 kpc)				—	—	—
Gas ring I (0–3.5 kpc)		—	—
Gas ring II (3.5–6.5 kpc)		—	—
Gas ring III (6.5–19 kpc)		—	—
Extended sources		—	—	—
Inverse Compton
Fermi bubbles		—	—	—
511 keV template		—	—	—	—
Naive model parameters, $N_{\text{param}}$		20253	78573	97838	104596	107639
Naive DOF		708747	650427	631162	624404	621361
Eff. model parameters, $N_{\text{param}}^{\text{eff}}$		1900	11800	10200	12700	12800
Eff. data bins, $N_{\text{data}}^{\text{eff}}$		624700	622700	620200	618800	619000
Eff. DOF, $k$		622800	610900	610000	606100	606200
$-2\ln {ℒ}_P$		1016041	637742	633334	627206	626998
$-2\ln {ℒ}_R$		14152	24652	23709	21640	20988
Model fidelity, $ℱ$		627	164	153	145	144

Table 2. Table 2 : Hyper-parameters for mock 2-component runs.

Components	free	constrained	smoothed
	Regularization hyper-parameters:
Gas
Inverse Compton

Equations

$${\varphi }_{pb}=\sum _k{T}_p^{(k)}{\tau }_p^{(k)}\cdot S_b^{(k)}{\sigma }_b^{(k)}\cdot {\nu }^{(k)},$$

$${\sigma }_b^{(k)},{\tau }_p^{(k)},{\nu }^{(k)}\ge 0.$$

$${\mu }_{pb}^{\mathrm{D}}=\sum _{p^{\prime }}{\mathcal{P}}_{bp{p}^{\prime }}{\mathcal{E}}_{b{p}^{\prime }}{\varphi }_{p^{\prime }b},$$

$${\mu }_{pb}^{\mathrm{P}}=\sum _s{\mathcal{P}}_{bp}({\mathbf{\Omega }}_{\mathbf{s}}){\mathcal{E}}_b({\mathbf{\Omega }}_{\mathbf{s}})\cdot S_b^{(s)}{\sigma }_b^{(s)}\cdot {\nu }^{(s)},$$

$$\ln \mathcal{L}=\ln {\mathcal{L}}_{\mathrm{P}}+\ln {\mathcal{L}}_{\mathrm{R}},$$

$$\ln {\mathcal{L}}_{\mathrm{P}}=\sum _{pb}{c}_{pb}-{\mu }_{pb}+c_{pb}\ln \frac{{\mu }_{pb}}{c_{pb}},$$

$$\begin{gathered} -2\ln {\mathcal{L}}_R=\sum _k{\lambda }_k{\mathcal{R}}_X({\mathbf{\tau }}^{(k)})+{\lambda }_k^{\prime }{\mathcal{R}}_X({\mathbf{\sigma }}^{(k)})+{\lambda }_k^{\prime \prime }{\mathcal{R}}_X({\nu }^{(k)})+{\eta }_k{\mathcal{S}}_1({\mathbf{\tau }}^{(k)})+{\eta }_k^{\prime }{\mathcal{S}}_2({\mathbf{\sigma }}^{(k)}) \\ +\sum _s{\lambda }_s^{\prime }{\mathcal{R}}_X({\mathbf{\sigma }}^{(s)})+{\lambda }_s^{\prime \prime }{\mathcal{R}}_X({\nu }^{(s)})+{\eta }_s^{\prime }{\mathcal{S}}_2({\mathbf{\sigma }}^{(s)}), \end{gathered}$$

$$\lambda {\mathcal{R}}_{{\ell }_2}(\mathbf{x})=\lambda \sum _i(x_i-1{)}^2,$$

$$\lambda {\mathcal{R}}_{{\ln }^2}(\mathbf{x})=\lambda \sum _i{\ln }^2{x}_i,$$

$$\lambda {\mathcal{R}}_{MEM}(\mathbf{x})=2\lambda \sum _{i}1-x_i+x_i\ln x_i.$$

$$\eta {\mathcal{S}}_1(\mathbf{x})=\eta \sum _{(p,p^{\prime })\in \mathcal{N}}(\ln x_p-\ln x_{p^{\prime }}{)}^2,$$

$$\eta {\mathcal{S}}_2(\mathbf{x})=\eta \sum _b(\ln x_{b-1}-2\ln x_b+\ln x_{b+1}{)}^2.$$

$$({\mathbf{\varphi }}^D{)}_i=(A^{(1)}\mathbf{\theta }{)}_i(A^{(2)}\mathbf{\theta }{)}_i(A^{(3)}\mathbf{\theta }{)}_i.$$

$${\mathbf{\mu }}^D=\sum _j{P}_{ij}({\varphi }^D{)}_j(\mathbf{E}{)}_j,$$

$$\text{EDM}\equiv \frac{||\nabla f|{|}_2^2}{2}\frac{||\delta \mathbf{x}|{|}_2}{||\delta \nabla f|{|}_2},$$

$${\mathcal{I}}_{ij}=-{⟨\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\ln \mathcal{L}⟩}_{\mathcal{D}(\mathbf{\theta })},$$

$$-\sum _{ij}\delta {\theta }_i\delta {\theta }_j\frac{{\partial }^2}{\partial {\theta }_i\partial {\theta }_j}\ln {\mathcal{L}}_P=\int d\Omega dE\frac{\mathcal{C}}{{\Phi }^2}\delta {\Phi }^2+\frac{\Phi -\mathcal{C}}{\Phi }\sum _k\frac{\delta {\tau }_k}{{\tau }_k}\frac{\delta {\sigma }_k}{{\sigma }_k}{\Phi }_k.$$

$$N_{\mathrm{H}}\simeq N_{\mathrm{HI}}+2N_{{\mathrm{H}}_2}=N_{\mathrm{HI}}+2X_{\mathrm{CO}}{W}_{\mathrm{CO}},$$

$$N_{\text{data}}=N_{\text{pix}}\times N_{\text{ebin}}.$$

$$N_{\text{data}}^{\text{eff}}\equiv {⟨-2\ln {\mathcal{L}}_P(\mathbf{\theta })⟩}_{\mathcal{D}(\mathbf{\theta })}=\sum _{i=1}^{N_{\text{data}}}2{⟨c\ln \frac{{\mu }_i(\mathbf{\theta })}{c}⟩}_{c\sim P({\mu }_i(\mathbf{\theta }))}=\sum _{i=1}^{N_{\text{data}}}d({\mu }_i(\mathbf{\theta })),$$

$$N_{\text{param}}^{\text{eff}}\equiv N_{\text{data}}^{\text{eff}}-{⟨-2\ln {\mathcal{L}}_P(\hat{\mathbf{\theta }})⟩}_{\mathcal{D}(\mathbf{\theta })},$$

$$\mathcal{F}=\sqrt{-2\ln {\mathcal{L}}_P(\hat{\mathbf{\theta }})-k}.$$

Full text

1. 1 Introduction

2. 2 SkyFACT – Concepts and definitions

3. 2.1 Multi-linear modeling with adaptive templates

4. 2.2 Poisson likelihood and regularization terms

5. 2.3 Parameter optimization with L-BFGS-B and convergence conditions

6. 2.4 Error estimation with Cholesky decomposition

7. 2.5 Convex optimization and non-degeneracy conditions

8. 3 Modeling Galactic diffuse emission

9. 3.1 Data selection

10. 3.2 Foreground modeling

11. 4 Towards a simple model for the Galactic disk emission

12. 4.1 Conventional diffuse emission components with nuisance parameters

13. 4.2 Additional spectral and spatial components

14. 5 Characteristics of the reference model

15. 5.1 Spatial modulation parameters

16. 5.2 Spectral modulation parameters

17. 5.3 The Galactic Center Excess

18. 6 Discussion

19. 6.1 The goodness-of-fit

20. 6.2 Future directions and potential applications

21. 7 Conclusions

22. A Tests with synthetic data

23. A.1 Image reconstruction, MEM and smoothing

24. A.2 Component separation and Poisson-induced bias with two components

25. A.3 Tests with Run5 mock data

1 Introduction

Two important activities in the analysis of astronomical images are parametric modeling and image reconstruction. Identifying adequate parametric models plays a central role in the astrophysical interpretation of observations, in particular for model discrimination and parameter regression. In the case of gamma rays, modeling observational data is very often done using template regression techniques, where fits with linear combinations of physically or observationally motivated spatial templates are used to perform component separation, often as function of energy [1, 2, 3, 4, 5, 6]. An orthogonal approach uses spectral templates to infer spatial characteristics of the various emission components [5, 7, 8]. Template regression does in general allow a straightforward incorporation of observations at other frequencies or sky regions.

A more ambitious approach is the full spatio-spectral modeling of gamma-ray emission, using numerical codes that simulate the generation, propagation and interaction of a multitude of cosmic-ray species in the Galaxy. This approach was pioneered with GALPROP [see, e.g., 9, 10, 11]; modern independent codes like DRAGON [12, 13, 14] and PICARD [15, 16] build on that success (for a description of similarities and differences between the codes see ref. [14]). Both template regression and the full modeling of gamma-ray emission are used for parameter regression, usually of up to a few dozen parameters [10]. The overall performance of the models can be considered as reasonably good, with residuals typically at the level of $\sim 30\%$ [10]. However, given the uncertainties in various emission components, which are, e.g., related to errors in estimated gas densities [17, 18, 19, 20], poorly constrained interstellar radiation field [21, 22, 23], simplifying assumptions in the physical cosmic-ray propagation models [e.g., 12, 24, 25], it is not a surprise that the models are far away from providing formally acceptable fits to the data [6].

One possible consequence of imperfect modeling is that estimators can be severely biased. At high Galactic latitudes, where the spatial characteristics of the background are reasonably well constrained and the relevant emission components (with the exception of, e.g., Loop I [26, 27]) are reasonably well known [see, e.g. 28], one could argue that biases are under control and comparable to statistical uncertainties. However, close to and in the Galactic disk, where modeling the Galactic emission becomes increasingly important and statistical errors are small, bias can be rather significant. Strong residuals along the Galactic disk [see, e.g. 6, 25] can potentially lead to wrong conclusions in particular about extended emission features. One possible approach to incorporate the characteristic magnitude of residuals as systematic error in statements about extended emission components was presented in [6], and based on a principal component analysis.

Model-agnostic image reconstruction is a central tool to aid the visual interpretation of complex data. It is hence no surprise that the use of image reconstruction techniques spans disciplines. Their relevance extends across many research fields, including optical interferometry [e.g., 29], the analysis of gamma rays [e.g., 30, 5], and image reconstruction in medical positron emission tomography (PET) [e.g., 31]. Interestingly, the underlying formal problems are often very similar, which allows the transfer of algorithms and techniques. Image reconstruction is usually based on penalized likelihood maximization (or, in a Bayesian interpretation, maximum-a-posteriori reconstruction), often using Poisson or Gaussian likelihoods or one of their bias-reducing variants [32].

The number of parameters is usually very high ($\geq 10^{4}$, often $\gg 10^{5}$), as parameters correspond to individual pixels in the image reconstruction domain. This requires additional regularization conditions to (a) reduce shot noise, (b) prevent over-fitting of the data, and (c) define unique solutions in potentially under-constrained problems [e.g., 33, 34, 35]. The regularization has the form of penalty terms in the likelihood, and its impact on the fit is controlled by various hyper-parameters. In a Bayesian interpretation, they correspond to priors on the regression parameters. Depending on the application, on the desired image properties, and on the available optimization algorithms, the regularization terms can take very different forms  [34, 36]. The arguably simplest regularizer is the $\ell_{2}$ norm, sometimes called ‘least square error method’ or ‘energy penalty’, which directly penalizes the variance of the model parameters [34, 33]. A traditional regularization technique popular in astronomy is the maximum entropy method (MEM) [30, 37, 38, 39, 40, 41, 1, 33]. It aims to minimize the configurational entropy, providing a reconstructed image that is as ‘uninformative’ and featureless as possible, while still being compatible with the data. However, non-smooth and/or sparsity-enhancing regularizers have become increasingly important, since they aid in feature selection. In particular the simple $\ell_{1}$ regularization enhances the sparsity of features in the image, and is connected to compressed sensing [42]. The similar ‘total variation regularization’ enhances the sparsity of the image in the gradient domain [34, 36]. More complex variants of this approach enforce sparsity in the wavelet space [43]. However, explicit smoothing regularizers, based on finite differences and the $\ell_{2}$ norm, are also used. The induced effective smoothing has an exponential profile and is in general non-uniform, since it depends on the varying signal-to-noise ratio of the image [44, 45, 46, 29].

In gamma-ray data analysis, a combination of various aspects of image reconstruction and parametric modeling is desirable for a number of reasons. In the case of existing spectral templates or models, fits in the spatial domain are practically an image reconstruction problem; a typical example is the Fermi bubbles [47, 48]. But even in situations where spatial information about an emission component is available (from, e.g., observations at other frequencies), the associated predictions have in general correlated uncertainties that can be fully modeled only with a very large number of nuisance parameters. An analysis of these nuisance parameters resembles again an image reconstruction problem. Indeed, various approaches that effectively combine aspects of image reconstruction and modeling exist in the literature. One example is the D3PO algorithm [5, 7], which is based on a Bayesian analysis using information field theory. The authors of this technique perform component separation, based on spectral information as well as priors on the smoothness of diffuse emission, to separate point source and diffuse components. Another example is ref. [1], where the authors perform an analysis of multi-frequency WMAP data, using spectral models and the MEM for spatial regularization.

In most gamma-ray analyses to date, nuisance parameters for spatial and/or spectral components of the models are not accounted for. Instead, a common approach is to ‘bracket uncertainties’ by focusing on a discrete set of plausible astrophysical configurations and analysis decisions. Each of these scenarios leads then to a different result for the observable of interest, and the envelope of all results defines the ‘systematic uncertainty band’ [e.g., 10, 28, 6, 49]. This approach is problematic in a number of ways. Usually, none of the models provides a statistically acceptable fit to the data. This raises the question why the consideration of various deficient models should lead to correct conclusions. Furthermore, it greatly hinders any rigorous model comparison, because the typically adopted statistical approaches (comparing the goodness-of-fit, considering likelihood ratios, or even Bayesian model comparison) exhibit undefined behaviour. In many cases, heuristic arguments can be made to support conclusions, but this does not replace a robust statistical approach. Lastly, the nuisance parameters are not only accounting for uncertainties in the model predications, but usually provide valuable information about the underlying astrophysical cause for the model deviations. This information is otherwise lost, or has to be guessed from residual maps. With the present work, we aim at providing a first step to overcome these limitations.

In this article, we propose a new hybrid approach, which attempts to generalize and combine the benefits of many of the above techniques. Our approach incorporates intrinsic uncertainties of various model components as penalty terms on spatial and spectral template modulation parameters. It effectively allows the application of image reconstruction techniques to excesses above partially constrained backgrounds. In particular, we aim to (a) facilitate component separation in scenarios where only partial knowledge about the spatio-spectral characteristics of the components is available, and (b) introduce a sufficient number of nuisance parameters in the analysis such that we can obtain formally good fits. A good fit to the data is a necessary criterion for reliable model comparison and bias-free parameter estimation. Our approach is based on a penalized Poisson likelihood [e.g., 46, 50] optimization with adaptive spatio-spectral (2-dim + 1-dim) templates. The spatial and spectral modulation parameters for the adaptive templates can be in general regularized with MEM, $\ell_{2}$ or explicit first- or second-order smoothing terms. Baring avoidable component degeneracies, we show that the underlying optimization problem is convex and hence has a unique (local and global) solution, despite the potentially millions of parameters. We perform analytic calculations of the gradient, which allows the use of the efficient optimization algorithm L-BFGS-B [51, 52, 53]. Uncertainties and covariances of the fit parameters are sampled from the inverse of the Fisher information matrix, using sparse Cholesky decomposition [54]. We study the properties of our approach with synthetic data.

As the first application of our method, we analyze gamma-ray data from the Fermi Large Area Telescope (LAT), focusing on a Region of Interest (RoI) along the Galactic plane including the Galactic center. Our analysis is based on predictions from the cosmic-ray propagation code DRAGON [14], as well as a number of ad hoc templates that absorb residual emission in the inner Galaxy. We focus here on the technical aspects of obtaining a good fit to gamma-ray data, and will present a detailed physical interpretation of the results elsewhere. One of the purposes of this work is hence to define a simple reference Galactic diffuse emission model that we can refer to (and further improve) in future work.

This paper is organized as follows: In section 2, we describe the adaptive template model, regularization terms in the likelihood function, the fitting strategy as well as our treatment of the covariance matrix of the problem. In section 3, we describe some details of the modeling of Galactic diffuse emission that we use in this work. In section 4 we provide a step-by-step example for how to construct a model for gamma-ray emission in our RoI. In section 5 we will discuss selected aspects of this model. In section 6, we discuss the fit quality of the reference model. In section 7, we present our conclusions. In appendix A we present several studies using synthetic data to estimate possible biases of the method.

2 SkyFACT – Concepts and definitions

SkyFACT (Sky Factorization with Adaptive Constrained Templates) is based on penalized maximum likelihood regression, using spatial and spectral templates with additional nuisance degrees of freedom that account for uncertainties of these templates. Penalization terms act as priors on the nuisance parameters. Errors are estimated by sampling from the inverse Fisher information matrix of the best fit point. We discuss model definition, penalization terms, the optimization procedure and error calculations in detail in this section.

2.1 Multi-linear modeling with adaptive templates

The diffuse component of the photon flux (per time and area) in energy bin $b$ and spatial pixel $p$ is given by a tri-linear model,

[TABLE]

where $T_{p}^{(k)}$ and $S_{b}^{(k)}$ refer, respectively, to the spatial and spectral template for emission component $k$. Furthermore, $\tau_{p}^{(k)}$ and $\sigma_{b}^{(k)}$ are spatial and spectral modulation parameters, respectively, and $\nu^{(k)}$ is an overall component normalization factor. Note that, in general, the overall normalization factor cannot be absorbed into the spatial and/or spectral modulation parameters: In some cases the spatial and spectral uncertainties might be very small, while only the overall normalization remains as unconstrained fitting parameter.

The only strict physical constraint on the modulation parameters and the normalization is non-negativity,

[TABLE]

Further constraints are introduced by regularization terms in the likelihood, and will be discussed below. For diffuse components, the expected photon count in pixel $p$ and energy bin $b$ is given by

[TABLE]

where $\mathcal{E}_{bp}$ denotes the exposure in the energy bin $b$ and $\mathcal{P}_{bpp^{\prime}}$ the effect of the instrument Point Spread Function, PSF (the probability that a photon from pixel $p^{\prime}$ is measured in pixel $p$). We neglect here the effects of energy dispersion, although there is no conceptual difficulty in adding this when required.

The modeling of point sources does not require spatial modulation parameters, so the model is only bi-linear. The expected number of signal photons in this case is given by

[TABLE]

where the sum is over point sources $s$ at angular locations $\boldsymbol{\Omega_{s}}$, and the PSF and exposure are evaluated at the source position. Note that the PSF matrix depends on the exact position of the source within the analysis pixel (whether it is at the center or at the edges). The total expected photon count is given by $\mu_{pb}=\mu_{pb}^{\rm D}+\mu_{pb}^{\rm P}$. To keep the notation compact, we use the symbol of the upper index ($s$ vs $k$) to discriminate between parameters related to point sources and related to diffuse components.

2.2 Poisson likelihood and regularization terms

The likelihood function used in this work has two major parts,

[TABLE]

where $\mathcal{L}_{\rm P}$ corresponds to the Poisson likelihood that confronts model predictions with observations, and $\mathcal{L}_{\rm R}$ contains various regularization terms for the modulation parameters.

The Poisson likelihood reads

[TABLE]

where $\mu_{pb}$ and $c_{pb}$ refer, respectively, to the expectation value and the observed number of events. We emphasize that we use the full Poisson likelihood in our fits, and do not assume Gaussian errors. The definition is slightly different from the commonly adopted Cash-statistic [55], since we are actually considering the likelihood ratio conditioned on the observed data (in the sense that $\ln\mathcal{L}_{P}=0$ for $\mu_{pb}=c_{pb}$). This version of the Poisson likelihood is commonly used in X-ray astronomy and is called the “c-statistic” in the widely-used spectral fitting tool XSPEC111https://heasarc.nasa.gov/docs/xanadu/xspec/ [56]. It facilitates goodness-of-fit calculations, but has no effect on parameter regression.

A central ingredient of our method are the large number of modulation parameters. They can either act as nuisance parameters that effectively account for systematic uncertainties in the adopted templates, or they permit functionality that is akin to image reconstruction. However, even in the case of image reconstruction, it remains important to regularize the modulation parameters in order to reduce the effects of Poisson noise and avoid over-fitting of the data. In the present work, the regularization terms are defined as

[TABLE]

where $X$ refers to the type of regularization, for each term separately. The first three regularization hyper-parameters, $\lambda_{k}$, $\lambda^{\prime}_{k}$ and $\lambda^{\prime\prime}_{k}$, control respectively the constraints on the spatial, spectral and overall modulation parameters, for the diffuse components. The regularization hyper-parameters $\eta_{k}$ and $\eta^{\prime}_{k}$ control the strength of the spatial and spectral smoothing, respectively. Explicit smoothing is for instance useful to reduce the impact of Poisson noise on the spatial modulation parameters, or to enforce a physically motivated smoothness of a given emission component. The point source regularization hyper-parameters, $\lambda^{\prime}_{s}$, $\lambda^{\prime\prime}_{s}$ and $\eta^{\prime}_{s}$, have an analogous functionality for point sources. We use the vector notation, $\boldsymbol{\tau}^{(k)}$ etc, to indicate that the regularization and smoothing terms depend in general on all parameters simultaneously. We briefly discuss three choices for the form of the regularization function $\mathcal{R}_{X}$, although only one of them (MEM) is used in the rest of this paper. The three choices are equivalent for large values of the regularization parameters, but differ in the tails, which is relevant if the regularization hyper-parameters are small and therefore the modulation parameters not tightly constrained.

The $\ell_{2}$-norm regularization is quite commonly used in the literature [34, 33]. For our parameter choices, it has the simple form

[TABLE]

where $i$ runs over the elements of $\boldsymbol{x}$. This regularization is a convex function, with the minimum at $x_{i}=1$; the associated variance is $1/\sqrt{\lambda}$. It can be used to constrain the element-wise variance of $\boldsymbol{x}$, and corresponds to Gaussian priors in a Bayesian interpretation. However, it behaves poorly when parameters are supposed to lie within a logarithmic range, e.g., one or two dex. For that reason we do not use this regularization in the present work.

A naive modification of the above $\mathcal{R}_{\ell_{2}}(\cdot)$ regularization that constrains parameters on logarithmic scales is obtained by applying the $\ell_{2}$ in log-space of the modulation parameters. This yields,

[TABLE]

and has a global minimum at $x_{i}=1$. The second derivative is given by $2\lambda(1-\ln x_{i})/x_{i}^{2}$. This means that again for large values of $\lambda$, this regularization is equivalent to the above linear $\ell_{2}$ one. For small values of $\lambda$, it corresponds to a quadratic penalty term in the log-space of $x_{i}$, which is useful to constrain parameter ranges within a few dex, if necessary. However, a disadvantage of the regularization term is that it is not convex (the second derivative becomes non-positive for $\ln x_{i}\geq 1$). This spoils the convexity of the entire optimization problem, potentially introducing multiple local minima. For this reason, we do not consider this regularizer further in the present work.

When the magnitude of an emission component is largely unknown and should be inferred from the data, a well motivated regularization is based on the maximum entropy method [30, 38, 39]. It reads

[TABLE]

The MEM regularizer is convex, and has a minimum at $x_{i}=1$. The variance for large values of $\lambda$, as derived from the second derivative at $x_{i}=1$, is $1/\sqrt{\lambda}$. It hence becomes equivalent to the $\ell_{2}$ and $\ln^{2}$ regularization in the large-$\lambda$ limit. For smaller values of $\lambda$, the main effect of this regularization term is to smooth out regions with low flux levels, while still permitting pronounced bright regions if they are preferred by the data. A good intuition can be obtained by noticing that its first and second derivatives are given by $2\lambda\ln x_{i}$ and $2\lambda/x_{i}$, respectively. For very small values of $x_{i}$, $\mathcal{R}(\boldsymbol{x})_{MEM}$ hence steepens infinitely and prevents $x_{i}$ from completely approaching zero. For very large values of $x_{i}$, on the other hand, the Poisson part of the likelihood can take control, such that compact bright regions remain relatively unconstrained. Due to its numerous practical properties, we use the MEM regularizer in the current work, not only for the spatial modulation parameters, but for all parameters. In figure 1, we show how the MEM regularizer behaves relative to the other regularizers discussed above.

The smoothness of components can be imposed by constraining the gradient of the modulation parameters instead of their range. For the spatial modulation parameters, this is enforced by terms of the form

[TABLE]

where $\mathcal{N}$ refers to the set of nearest-neighbour pixel pairs. Here, only pairs are included where both pixels are actually part the template (which can only cover a fraction of the sky). Note that we constrain the gradient of the logarithm of the modulation parameters instead of the modulation parameters directly. This has the advantage that the overall normalization of the modulation parameters drops out. The disadvantage is a loss of convexity, like for the above $\ln^{2}$ regularization. This will be discussed below.

The spectral smoothing terms that we adopt in this work, on the other hand, constrain the second derivative of the logarithm of the modulation parameter, and not the gradient. The corresponding regularization term has the form

[TABLE]

Assuming that the energy bins are logarithmically separated, this regularization favours power-law spectra. Since many astrophysical spectra actually have the form of a (potentially rolling) power-law, this is a reasonable behaviour that favours in many cases ‘more physical’ spectra. However, it has to be used with care since it will cause an automatic ‘power-law extrapolation’ when the shot-noise level increases at high energies.

The strength of all regularization terms are controlled with the regularization hyper- parameters $\lambda_{k}$, $\lambda^{\prime}_{k}$, $\lambda^{\prime\prime}_{k}$, $\eta_{k}$ and $\eta^{\prime}_{k}$ for diffuse components, and $\lambda^{\prime}_{s}$, $\lambda^{\prime\prime}_{s}$, $\eta^{\prime}_{s}$ for point sources. We will discuss below, and in appendix A, a few examples of typical values and their effects on the resulting reconstruction.

2.3 Parameter optimization with L-BFGS-B and convergence conditions

Internally, the tri-linear model for the calculation of expected photon counts from diffuse emission components is represented by three matrices $A^{(1)}$, $A^{(2)}$ and $A^{(3)}$. These map model parameters $\boldsymbol{\theta}\equiv(\boldsymbol{\tau}^{(k)},\boldsymbol{\sigma}^{(k)},\nu^{(k)},\boldsymbol{\sigma}^{(s)},\nu^{(s)})^{T}$ onto the expected photon fluxes $\boldsymbol{\phi}^{D}\equiv(\phi_{bp})$ via

[TABLE]

The transformed vectors are multiplied element-wise, as indicated by the index $i$. Although the matrices $A^{(1,2,3)}$ are very large (in our example of the order $10^{5}\times 10^{5}$), they are rather sparse (in our example with a sparsity well above $99\%$) and can be efficiently stored in the compressed sparse column format. The use of optimized sparse matrix routines leads then to very fast evaluation times. Another advantage of the above representation is its simple analytical form, which facilitates the calculation of gradients and the Fisher information matrix. The final expected number of photons, $\boldsymbol{\mu}^{D}=(\mu_{bp}^{D})$ are obtained from the matrix product

[TABLE]

where $P$ and $\boldsymbol{E}$ represent the PSF and the exposure map convolution in data space respectively. For the PSF convolution, we assume periodic boundary conditions for simplicity, which leads to minor artefacts at the borders of the RoI, but is not relevant for the present discussion (although it will some impact on the formal goodness-of-fit). Point sources are calculated in a similar way, but here two matrices $B^{(1)}$ and $B^{(2)}$ are sufficient, since the corresponding model is only bi-linear. For the models below, we find that the evaluation of the likelihood function and the gradient takes less than a second on one core of a typical modern computer. Sparse matrix multiplication is not straightforward to parallelize on shared or distributed memory systems, which makes it somewhat challenging (but certainly not impossible [e.g., 57]) to reduce the computation time further, even on multi-core systems.

Parameter optimization in SkyFACT is performed with the L-BFGS-B method [51, 52, 53]. The L-BFGS-B algorithm is based on the quasi-Newtonian method of Broyden, Fletcher, Goldfarb and Shanno (BFGS) [58]. Its limited-memory version (L-BFGS) does not directly store numerical approximations to the Hessian matrix, but rather a number of previously calculated gradients to construct approximate Hessian matrices on the fly. This makes it possible to use it for very large number of parameters ($10^{6}$ and more). Furthermore, the version used here (L-BFGS-B) supports parameter boundaries, using the active-set method. This is essential in order to be able to impose the non-negativity constraints. We use here the implementation from Refs. [52, 53, 59]. The algorithm takes both the objective function and its derivative as input, which are both straightforward to calculate based on the above analytical expressions.

There are two main convergence criteria that we consider in this work. In the simplest case, one can just require that the vertical step-size of the last iteration is below a certain threshold value (this is one of the criteria of the adopted L-BFGS-B algorithm). A more informative quantity is however the estimated (vertical) distance to the minimum (EDM), which is known from the popular Minuit algorithm [60]. To estimate the EDM properly requires information about the full Hessian, which is however not easily available outside the L-BFGS-B implementation that we use. We therefore adopt an extremely simple heuristic definition of an ‘effective EDM’ that is motivated by the one-dimensional case. It is defined as

[TABLE]

where $\delta\boldsymbol{x}$ refers to the parameter changes in the last step, $\delta\boldsymbol{\nabla}f$ to the change of the gradient, and $||\cdot||_{2}$ is the $\ell_{2}$ norm. We find that although it is only a crude (noisy, not monotonically decreasing) estimate of the vertical distance to the minimum, it is a useful dimensionless estimator to determine how close the fit is to the global minimum.

2.4 Error estimation with Cholesky decomposition

Estimating parameter uncertainties requires knowledge of the full covariance matrix of the model. Here, we approximate the covariance matrix by the inverse of the Fisher information matrix. The Fisher matrix encodes the expected error, and is independent of the data (it is numerically faster to calculate than the full Hessian matrix of the likelihood function, and numerically more stable even when the minimum of the system has not been exactly reached). Formally, it is given by

[TABLE]

where the average is taken over mock data from the best-fit model, $\mathcal{D}(\boldsymbol{\theta})$. The average can be taken analytically, and the resulting expressions are relatively simple to handle (see, e.g., ref. [61] for a detailed discussion in context of Poisson likelihoods). Since the Fisher matrix is composed of sparse matrices, it is sparse as well, at the 99% level in the examples discussed below.

We estimate the errors of individual model parameters as well as predicted fluxes and other quantities using a sampling method. This circumvents the need to explicitly calculate the covariance matrix (this would be computationally very costly, since the inverse of a sparse matrix is not sparse itself). The procedure is based on sparse Cholesky decomposition [54]: The Fisher information matrix is decomposed in the form $P\mathcal{I}P^{T}=LDL^{T}$, where $L$ is a lower unit triangular matrix, $D$ is a positive-definite diagonal matrix, and $P$ is a fill-reducing permutation matrix that maximizes the sparsity of $L$. Based on the matrices $L$, $D$ and $P$, one can efficiently sample model parameter vectors $\delta\boldsymbol{\theta}$ such that their covariance is given by the inverse, $\mathcal{I}^{-1}$. These model parameters describe then deviations from the best-fit model, similar to samples from Bayesian posterior distributions.

The sampling procedure works as follows. We select a vector of standard normal distributed parameters, $\boldsymbol{x}\sim\mathcal{N}(0,1)$, with the same length as the number of model parameters. Next, we obtain $\boldsymbol{y}$ as solution of $L^{T}\boldsymbol{y}=D^{-1/2}\boldsymbol{x}$, which can be done by back-insertion. Finally, we calculate $\delta\boldsymbol{\theta}=P^{T}\boldsymbol{y}$. We repeat the process as many times as necessary to obtain a reasonably-sized sample. As mentioned above, the full covariance matrix, $\Sigma_{ij}$, of the model parameters $\boldsymbol{\theta}$ can be estimated from $\langle\delta\theta_{i}\delta\theta_{j}\rangle\approx\Sigma_{ij}$, where the average is taken over many samples. In practice, we are usually interested in mean values and variances of model fluxes, which are functions of $\boldsymbol{\theta}$. Mean values for fluxes and other model predictions are derived from the best-fit value of $\boldsymbol{\theta}$. Variances are derived by calculating model predictions for $\boldsymbol{\theta}+\delta\boldsymbol{\theta}$ and averaging over many samples.

2.5 Convex optimization and non-degeneracy conditions

Convex optimization is the minimization of convex functions, and has a number of useful properties. In particular, any local minimum is also a global minimum, meaning that there is only one best-fit solution to the problem. We will show here that – as long as the model components remain sufficiently non-degenerate – the objective function is indeed convex. We will point out qualitatively under what conditions problems can occur (one can also derive quantitative conditions, but the resulting expressions are somewhat lengthy and we do not reproduce them here).

We first consider the Poisson likelihood $\mathcal{L}_{P}$. To simplify the notation, it is useful to write the $\ln\mathcal{L}_{P}$ in the unbinned form (we refer to ref. [61] for a discussion of unbinned Poisson likelihoods). For simplicity, we neglect point source degrees of freedom (which are very sparse in the sky and not expected to be a problem in the present scenario), as well as the overall normalization factor. We are in particular interested in the differential version of the Hessian. It is given by

[TABLE]

Here, we defined the ‘unbinned counts map’ $\mathcal{C}(\Omega,E)\equiv\sum_{i}\delta^{(2)}(\Omega-\Omega_{i})\delta(E-E_{i})$, where $E_{i}$ and $\Omega_{i}$ are the energy and position of photon $i$, and $\delta(\cdot)$ refers to Dirac delta distributions. Furthermore, $\delta\Phi\equiv\delta\boldsymbol{\theta}\cdot\boldsymbol{\nabla}_{\theta}\Phi$, and $\delta\tau_{k}$ and $\delta\sigma_{k}$ are defined analogously. Finally, $\Phi_{k}$ is the contribution to $\Phi$ from component $k$.

It is clear that the first term in Eq. (2.17) is always non-negative, and it is strictly positive unless two or more model components are exactly degenerate (in that case some change in component $k$ can be exactly undone by the opposite change in component $k^{\prime}$, leading to $\delta\Phi=0$ in that particular direction of the parameter vector $\delta\boldsymbol{\theta}$). The second term can, depending on $\mathcal{C}$, however become negative and spoil the convexity of the objective function. However, it depends on the simultaneous change of spectral and spatial modulation parameters of the same component. It is exactly zero for any component that has either spatial or spectral degrees of freedom, but not both. As long as the model consists only of components that each have either a fixed distinct morphology or spectrum, the optimization problem is strictly convex.

However, if both spectral and spatial variations are allowed for at least two components in the fit, one can find in general degenerate transformations with $\delta\boldsymbol{\Phi}=0$ while $\delta\tau_{k}\neq 0$ and $\delta\sigma_{k}\neq 0$. In that case, the second term of Eq. (2.17) can potentially become negative, and the objective function could have multiple local minima. This corresponds to the case where components have so much spectral and spatial freedom that they can effectively exchange their roles in the fit. Situations like that are usually not desired, and should be prevented by the regularization of the modulation parameters.

As discussed above, both the MEM and $\ell_{2}$ regularizations are completely convex, and hence cannot spoil the convexity of the problem. Since we use the MEM regularization, increasing the value of the corresponding hyper-parameters $\lambda$ can only increase the convexity of the problem, while breaking the degeneracy between model components.

Note, however, that the smoothing terms $\mathcal{S}$ are not generally convex. This is caused by the fact that the logarithm of the modulation parameters (e.g., $\ln\tau_{p}^{(k)}$) enters their definition, which is similar to the above $\ln^{2}$ regularization. For this reason, it is advisable to check explicitly that the global minimum is actually reached. This can be done by confirming that reducing the smoothing hyper-parameters $\eta$ and $\eta^{\prime}$ does not qualitatively change the result (although it will reduce the smoothness of the components in the fit; see appendix A.3.)

3 Modeling Galactic diffuse emission

3.1 Data selection

We use $7.6$ years of Fermi-LAT Pass 8 data from 2008 August 4 and 2016 March 3. We apply standard data cuts to events recommended by the Fermi Science Support Center (“(DATA_QUAL>0) && (LAT_CONFIG==1)” and zenith angles $>90^{\circ}$). We select only ULTRACLEAN events (evclass=512) and use both FRONT and BACK converted events (evtype=3). We use the Fermi Science Tools222http://fermi.gsfc.nasa.gov/ssc/data/analysis/ (v10r0p5) to build count and exposure maps binned into 25 log-spaced energy slices in the energy range $0.34-228.65$ GeV. Our RoI spans $|\ell|\leq 90^{\circ}$ and $|b|\leq 20.25^{\circ}$. The data are binned in Cartesian coordinates with a pixel size of $0.5^{\circ}$. In figure 2, we show the data used in this analysis. Note that the adopted pixelization ensures that the bright Galactic disk falls onto the central row of pixels. Since we use a relatively small number of pixels in the analysis (the number of spatial nuisance parameters would otherwise explode in our current treatment), using a pixelization scheme that takes into account the characteristic emission gradient towards the disk is particularly important.

3.2 Foreground modeling

We model the foreground components following standard methods and codes. In particular, we rely on the prescriptions of the GALPROP [62] and DRAGON [63] codes, and refer to the Fermi-LAT Collaboration works, e.g. [10]. Details of the modeling are explained below.

The gamma-ray emission from interactions of cosmic rays (CR) with the interstellar medium (ISM) and low-energy ambient photons is the so-called Galactic diffuse emission (GDE). At GeV energies, the GDE originates mainly from three processes: (1) the decay of neutral pions produced by collisions of CR protons with ISM; (2) the bremsstrahlung radiation of electrons in the electric field of ISM charged particles; and (3) the up-scattering to GeV energies of the low-energy ambient photons of the interstellar radiation field (ISRF) due to Inverse Compton scattering (ICS) off CR electrons.

The GDE depends on the production of CR is the Galaxy, their distribution and spectra, as well as on the propagation mechanisms in the ISM and the interactions with the environment. In the present work, we model the GDE in a relatively minimalistic way, and explore what ingredients we need to obtain a reasonably good fit to the data (after nuisance parameters are taken in account).

Hadronic $\pi^{0}$ emission.

The spatial distribution of the hadronic emission from $\pi^{0}$ decay traces closely the target material distribution, i.e. the distribution of hydrogen in the Galaxy. The hydrogen in the Galaxy is in the form of atomic (HI), molecular (H2), and ionized (HII) gas, traced by different observables, namely the 21 cm line, the 2.6 mm CO emission line, and pulsar dispersion measurements respectively [64]. The most abundant component is represented by the molecular and atomic gas (almost 70% by mass), while ionized H is subdominant in mass (only few percent on average).

The total H column density (in cm*-2*) can be approximated by the sum of atomic (HI) and molecular hydrogen (H2):

[TABLE]

where the molecular H2 column density is assumed to be proportional to the integrated line intensity of CO, $W_{\rm CO}$ (in K km/s ), through the $X_{\rm CO}$ coefficient (in cm*-2*/(K km/s)). The $X_{\rm CO}$ is typically considered to be constant thorough the Galaxy. However, there are no a priori reasons for it to be such and, indeed, gamma-ray analyses have shown that, most likely, it is spatial dependent [65, 20, 66, 25].

In the present work, we build the $\pi^{0}$ template as the sum of the gas column densities for atomic and molecular hydrogen, available within the GALPROP public release.333http://galprop.stanford.edu The HI column density and $W_{\rm CO}$ integrated line intensity maps are based on the corresponding 2D analytical models developed in ref. [62], renormalised to the Leiden-Argentine-Bonn 21 cm line survey [67] and to the 2.6 mm CO line survey [68], respectively. For the HI, the spin temperature used to correct for opacity of the 21 cm line is $T_{S}=125$ K. No dust-reddening correction is applied to the column density. Furthermore, we do not apply a dark gas correction, since we are instead interested in seeing how dark gas affects the spatial modulation parameters associated with hadronic emission [69, 70, 71]. The GALPROP gas maps are binned in Galacto-centric radial annuli, using the distance information extracted from the gas velocity and rotation curve (for details about the methodology see appendix B of ref. [10]). In practice, we obtain the gas template by adding the hydrogen column densities as: $N_{\rm HI}+2N_{\rm H_{2}}$, assuming a constant $X_{\rm CO}=1.9\times 10^{20}$ cm*-2*/(K km/s). We build three gas templates using the following radial binning: the first gas ring corresponds to the emission within 0–3.5 kpc from the GC, the second ring to 3.5–6.5 kpc, and the third one to 6.5–19 kpc. The advantage of using radial rings is that we can combine them arbitrarily to get maps for different spherical shells about the GC and study how this affects our conclusions.

Leptonic ICS emission.

To model the ICS emission we use the publicly available DRAGON code [63] and GammaSky [24, 72], a custom numerical package interfaced to DRAGON to compute the gamma-ray emission from the CR distribution output of DRAGON. We use the “Ferrière” model for the source distribution from the primary components [64], and propagation parameters corresponding approximately to the “KRA4” model in ref. [72].444 We use the run_2D_KRA.xml model file available at https://github.com/cosmicrays/DRAGON/examples. The ISRF model adopted by DRAGON/GammaSky is documented in [73] and available for download from the GALPROP public website (v54).

Leptonic bremsstrahlung emission.

Starting from the same DRAGON model file and using the GALPROP gas maps as target for interactions, we model the bremsstrahlung emission of electrons and positrons off interstellar ions and nuclei. The gamma-ray emission from bremsstrahlung is typically suppressed with respect to that from $\pi^{0}$ and ICS. However, it might play a significant role in the very inner region of the Galaxy; see for example ref. [74]. We tested the effect of including this additional template in our baseline model, without improvement of the fit that would be relevant for the current work (for our purposes, it is enough that bremsstrahlung can be absorbed by the $\pi_{0}$ templates). Therefore, we do not include the bremsstrahlung template in our baseline model.

The spatial templates used as starting points in the fits are shown in figure 3. The energy range shown is 0.96 – 12.98 GeV, chosen to highlight both the $\pi^{0}$ bump as well as a significant fraction of the Galactic center excess emission (discussed below).

Although DRAGON also produces the energy spectra of the various emission components, we do not use them in the current work. The energy spectra depend critically on the energy spectra of hadronic and leptonic CRs, and hence on propagation and injection. Instead, we use the $\pi^{0}$ and ICS component spectra from [10] as references. These spectra are obtained from the gamma-ray analysis of CR propagation models that are also in agreement with local CR measurements, and hence represent for us a good starting point. However, CR spectra in the Galaxy might be in general quite different from what we measure locally. For this reason, we use the reference spectra only as starting points in the fit, and allow considerable freedom such that they can adapt to the gamma-ray data if needed.

4 Towards a simple model for the Galactic disk emission

In this section, we derive a simple reference model for Galactic diffuse emission in our main RoI. Since even this reference model will have $10^{5}$ parameters, the meaning of the term ‘simple’ is not automatically evident. What we actually aim to find is a minimal set of physical components, with distinct spatial or spectral characteristics (these contribute a very low number of parameters), augmented with realistically sized nuisance parameters (this number is necessarily very large). The goal is to remove most of the residuals in our RoI in a physically motivated way. Our reference model is not meant to be the final answer to modeling all of the gamma-ray emission in our RoI, but rather as a base for a number of follow-up studies.

The starting point for our reference model will be the numerical predictions for hadronic and leptonic emission in the Galaxy from DRAGON, based on locally measured CR fluxes and HI and CO gas maps, as discussed in section 3. We will study how much of the observed gamma-ray emission is already explained by these a priori predictions, both directly and when augmented with additional nuisance parameters that account for expected spatial and spectral uncertainties in the predictions. For a number of irreducible residuals, particularly in the bulge region, we will introduce additional templates.

4.1 Conventional diffuse emission components with nuisance parameters

In table 1, we list five increasingly complex models, their components and the associated regularization parameters (Run1 to Run5). They provide a step-by-step illustration for how additional components and nuisance parameters account for residuals when modeling the gamma-ray emission observed by Fermi-LAT in our RoI. We will discuss each of the five models separately.

Run1. This is our initial model, mostly based on the a priori information that enters the modeling with DRAGON. It includes just one gas and one ICS component, along with a component for the isotropic gamma-ray background (IGRB) as measured in ref. [28] and the 3FGL point sources [75]. In the case of the IGRB, we allow spectral uncertainties with a regularization of $\lambda=16$ (this corresponds to $1/\sqrt{\lambda}=25\%$ uncertainty) that are somewhat smaller than the systematic errors quoted in ref. [28]. This prevents large deviations from the measured spectrum. The spectra of the gas and ICS components are constrained to be close to the spectra corresponding to the locally measured CRs (again within $25\%$), and are taken from ref. [10]. For all of the components of Run1, except the IGRB, we keep the overall normalization, $\nu^{(k)}$, free in the fit. The same is the case for the point sources. Furthermore, we use a weak regularization to keep the source spectrum close to its original spectrum in the 3FGL (to within $20\%$). The total number of model parameters is 20253 (77 for the diffuse components, and the remaining 20176 for the 776 3FGL point sources). We show the resulting spectra and the spectral residuals in figure 4. It can be clearly seen that very significant residuals remain in particular in the central region. The overall spectrum integrated over the full RoI is better recovered, however only at the expense of a strongly distorted IGRB spectrum.

Furthermore, in figures 5 and 6, we show in the top panels the significance of the residuals when subtracting the best-fit Run1 model from Fermi-LAT data. As expected, we find that this simple scenario reproduces observations, with fractional residuals around $\sim 30\%$ in many regions of the sky. In particular, large and very significant residuals remain along the bright Galactic plane, and at higher energies in the Fermi bubble regions and towards Cygnus X (at $\ell\sim 80^{\circ}$). A histogram of the residuals is shown in figure 7, and will be discussed further below in context of Run5.

Run2. As a next step, we introduce also spatial modulation for all components, except the IGRB, with values for the regularization hyper-parameters listed in table 1. We introduce $33\%$ ($\lambda=10$) uncertainties for the morphology of the gas component, and $100\%$ ($\lambda=1$) uncertainties for the morphology of the ICS component. We also now include smoothing regularization terms, with stronger smoothing for the ICS component than for the gas. The ICS component is smoothed with a hyper-parameter of $\eta=100$, which corresponds to roughly $10\%$ variations between neighboring pixels, while the gas component is smoothed with $\eta=25$ ($20\%$ variations). The total number of parameters is now 78548. As in Run1, 20176 of these parameters correspond to point sources and their spectra. Most of the remaining parameters correspond to spatial modulation parameters for the ICS and gas components.

Looking at figures 5 and 6, it is apparent that accounting for spatial and spectral nuisance parameters in the fit reduces the amount of residuals in the sky drastically. Typical fractional residuals are now of the order $\lesssim 10\%$ at energies above 1.24 GeV. Fractional residuals above the minimum energy of 0.34 GeV are practically zero (this can also be seen in figure 5), because there are enough degrees of freedom in the fit to remove the most significant residuals. They happen to be at the lowest energies due to the larger number of low-energy photons.

One effect of the nuisance parameters is that reducible residuals, in the sense that they can be absorbed by small changes of the local intensity of model components without significantly changing their spectra, are largely removed, whereas irreducible residuals remain. Structures that become now evident are the bubble-shaped positive residuals that appear from the second energy band on, and, at high energies, numerous localized residuals along the Galactic plane.

Run3. This model is similar to the previous one, except that we split the gas component into three rings that cover the radial ranges 0–3.5 kpc, 3.5–6.5 kpc and 6.5–19.0 kpc. Details about the construction of gas rings can be found in section 3. The main effect is that residuals along the Galactic disk are further reduced, as can be seen in the third row of figures 5 and 6. However, localized residuals remain along the Galactic plane, both close to and at the Galactic center, as well as further away. Most of these residuals grow stronger at higher energies, indicating rather hard fluxes. Indeed, many of the residuals correspond to extended emission associated with the positions of various extended sources observed by Fermi-LAT. Below, we will be able to absorb most of these sources with extended source templates from the Fermi 3FGL.

In figure 8, we show the modulation parameters that correspond to the four diffuse components in Run3. For the gas ring I and gas ring II, the modulation parameters stay close to one except for a few bright spots along the Galactic disk. Most of these correspond to extended sources in the Fermi 3FGL. However, in gas ring III, which reaches up to high latitudes, we can see the effects of ‘dark gas’, which is not included in our gas templates and typically forms extended halos around CO and HI emission regions [71]. This gas component is not captured by the usual CO line or 21 cm observations, but can be inferred from gas extinction maps to good accuracy. We did not include such a correction in our initial maps shown in figure 3. Observing that the nuisance parameters are in fact recovering the dark gas from the gamma-ray data alone indicates that our method works as expected. In the ICS template, we also clearly see an additional component in the modulation parameters: the Fermi bubbles. The spectrum for the Fermi bubbles is somewhat similar to the ICS spectrum, but it is, of course, a distinct component that we will include as such in subsequent fits.

4.2 Additional spectral and spatial components

Run4. Up to now, we considered only the most conventional ingredients in the modeling of the GDE, although augmented with nuisance parameters. We found in Run3 that various significant residuals still remain. We will now add additional components that account for these residuals. We note that this does not mean that this is the only valid way to decompose the observed gamma-ray emission from the Galaxy in individual components.

We first add a spectral template for the Fermi bubbles. The Fermi bubbles component has an unconstrained morphology, while the spectrum is fixed, to within $5\%$ ($\lambda=400$), to the Fermi bubble spectrum measured in ref. [48]. We furthermore use a weak spatial smoothing regularization with $\eta=4$. In this case, we use non-zero smoothing regularization terms to reduce the effect of Poisson noise on the reconstructed morphology, although this is also expected to somewhat wash out potentially sharp features in the bubble edges. We constrain the bubble template to be non-zero only in an hour-glass shaped region defined by $|\ell|+4\cos(b/6.4^{\circ})<16^{\circ}$.

Furthermore, we add templates for the many spots of localized extended high-energy emission that we saw in Run3 along the Galactic disk. Many of these regions are colocated with the extended sources in the 3FGL. 18 such sources, out of 25 total in the 3FGL catalog, fall within our RoI. We therefore introduce new diffuse components that are only non-zero in regions that correspond to the extended sources in the 3FGL catalog that fall within our RoI. We initialize the spectra for these components to the extended source spectra from ref. [75], but leave considerable freedom for the spectra in the fit. Furthermore, the morphologies of the extended sources are left unconstrained, except for some smoothing. The templates of the extended sources are restricted to circular regions centered on their 3FGL locations with radii varying from 1.5–5*∘*, depending on the size of the residuals left by the source in Run3. We use a total of 16 templates for the extended sources. We exclude two sources, $\gamma$-Cygni and HESS J1616-508, because their templates overlap with the templates of other extended sources, and they are therefore completely degenerate with them, given the freedom allowed in both the spectral and spatial modulation parameters. Details about the adopted hyper-parameters can be found in table 1.

We can see in figures 5 and 6 that the additional components remove much of the residual emission that remained in Run3. Besides removing the Fermi bubbles at high latitudes, most of the bright high-energy emission along the Galactic disk is removed. However, not all of the residuals disappear. The most notable one that remains at high energies is localized at $\ell\approx-48^{\circ}$ is extended by a few degrees. This residual is most likely associated with an extended source or sources not present the 3FGL catalog. In fact, there are two extended sources within $|b|<0.5^{\circ}$ near $\ell=-48^{\circ}$ in the recently released 3FHL catalog [76], 3FHL J1420.3-6046e and 3FHL J1409.1-6121e (see also the FGES catalog, built from a dedicated search for extended emission above 10 GeV along the Galactic plane: [77]). Overall there are 19 extended emission sources in the 3FHL catalog that are within our RoI that are also not detected as extended emission in the 3FGL catalog. These sources may very well be contributing to some of the high-energy residuals we find along the disk in our fits. However, we restrict ourselves to the 3FGL catalog for now and save a full analysis of the extended emission sources along the disk for a future study.

Another residual that remains is some extended emission close to the Galactic center. It has the peculiar feature that it is strongest in the second energy band, and appears almost as a negative residual in the third and fourth energy band. We show this feature in figure 9. The top two figures are the same as in figures 5 and 6, but more heavily smoothed with a gaussian kernel to reduce the visual effect of noise and enhance the larger residual present close to the Galactic center. In the bottom panel of figure 9, we show the ratio of the higher-energy residual map divided by the lower-energy map, after the maps were smoothed, to highlight the strength of this residual (the sharp positive-negative edges in the ratio map are the result of the smoothing). We address this residual in Run5.

Run5. This is the last and most complex model in the present analysis. On top of the components of Run4, it includes a template that accounts for the central bulge excess that we found above. It corresponds to the Fermi Galactic center excess (GCE) extensively studied in the literature [e.g., 78, 79, 80, 81, 82, 83, 13, 6, 8, 11, 84]. As template for the GCE, instead of an analytic DM-inspired model as typically done in the literature, we use here, for the first time, the profile of the 511 keV flux as measured by INTEGRAL  [e.g.  85, 86, 87, 88]. The positron annihilation signal (for a comprehensive review see ref. [89]) in the Galactic bulge shows a few common features to the GCE: both are almost spherically symmetric around the Galactic center, show a spectrum uniform in an extended RoI, and peaks in the direction of the center of the Galaxy. For a more thorough study about the connection between the two signals we refer the reader to a forthcoming publication [90]. We here report the main results of introducing such an additional template component. We only include the Galactic center and bulge component of the 511 keV flux, and neglect the disk for the present purpose. These components are modeled following ref. [88]. Furthermore, we keep the morphology of the component constrained to within $20\%$, while the spectrum remains completely unconstrained, see table 1. The corresponding residuals in figures 5 and 6, show that this component is enough to remove the remaining positive and negative residuals in the central region.

In figure 7, we show the distribution of residuals in units of standard deviations for Run1 and Run5. Comparing Run1 to Run5, we find that the introduction of nuisance parameters in the fit reduces the number of pixels with highly significant residuals (above/below $\pm 3\sigma$) by over two orders of magnitude. In fact, the most drastic reduction of residuals appears from Run1 to Run2, which can be seen in the drop of the Poisson $-2\ln\mathcal{L}_{P}$ as shown in table 1.

However, even in Run5 there are residuals above $\pm 5\sigma$ that remain. On closer inspection, we find that at energies $\lesssim 10{\rm\,GeV}$, they are usually related to regions around bright point sources which are over- or under-subtracted. A likely cause is that we do not refit the positions of the 3FGL sources in our analysis; another may be due to potential, small numerical imperfections in our treatment of the Fermi PSF. At high energies, the remaining residuals often correspond to the regions that are also visible in the residual maps in figure 6: extended regions along the Galactic disk with hard emission. The modeling of this emission component is still not optimal in the current analysis, and could be further improved by, e.g., increasing the size of the patches that we use for extended sources in our Run5 or including templates for additional known extended sources in the 3FHL.

Lastly, the distribution of residuals in figure 7 strongly deviates from a naive normal distribution, illustrated by the black line, which would be approximately expected if the residuals were due to statistical fluctuations only. Additionally, there is a strong tilt towards positive residuals, particularly for energies above 60 GeV, for both Run1 and Run5. This can be understood by realizing that the smallest possible negative residual in a pixel with an expected number of $\mu$ photons is, in our adopted schema, $-\sqrt{\mu}$. At high energies, $\mu$ can be of order one, leading to the apparent tilt. In fact, we observe the same behaviour for our mock results that are discussed in appendix A.

5 Characteristics of the reference model

In this section, we discuss in detail the aspects of Run5 that we motivated in the previous section. In particular, we are interested in the plausibility of the values and variances of the modulation parameters in the analysis.

5.1 Spatial modulation parameters

The various best-fit fluxes for the extended components of Run5 are shown in figure 10. The spatial modulation parameters of Run5 are shown in figure 11, for the four conventional diffuse emission components. Histograms for the same modulation parameters are shown in figure 12, weighted by pixel intensity. The original morphology of the four components is shown in figure 3, where the spatial modulation parameters are set to one.

We find that for the innermost gas ring I, the nuisance parameters remain rather close to one. The weighted $90\%$ central quantile covers the range [0.88, 1.04], indicating typically less than $10\%$ deviations from the original template. Weighting these parameters is important, since a large number of pixels of the template correspond to very small fluxes, which are not constrained by the data and hence not rescaled. This does not, however, necessarily reflect that the original template was very accurate, but more likely that measuring this template is difficult in the bright Galactic disk.

The gas rings II and III have modulation parameters with a broader distribution. Their weighted $90\%$ central quantiles cover [0.80, 1.46] and [0.70, 1.45], respectively. Note that there are also a few outliers, which likely correspond to isolated point sources that get incorrectly absorbed by the template. These do not affect the quantiles, which are relatively resilient towards outliers. The larger variations indicate that the templates are well constrained by the data, and give an estimate for the quality of gas maps in describing the data.

The gas ring III, which includes the position of the Sun and is responsible for the structured emission at mid-latitudes, exhibits an interesting pattern in the nuisance parameters. The spatial modulation parameters are in some regions as large as 2–3 (in the tail of the distribution); in other regions we observe a suppression of similar magnitude. As for Run3, the emerging structures are very similar to the ‘dark gas’ corrections [71] (see discussion for Run3).

In figure 11, we also show the rescaling of the ICS emission w.r.t. the initial template from figure 3. The weighted $90\%$ quantile covers here the values [0.84, 1.21], which is a smaller range than for the gas emission. However, there are some similarities between the large-scale features of the gas ring III and the ICS rescaling map. This suggests that both components are not completely separable, and probably some of the ICS emission is falsely absorbed by the gas map and vice versa. This is expected if, for instance, the spectra of the two components are not perfectly uniform throughout the entire RoI (this effect could be mitigated by further splitting up these diffuse components). On the other hand, looking at the total ICS emission shown in figure 10 suggests that the emission absorbed by our ICS template might actually have a bi-modal ring-like structure, with enhancements towards $\ell\approx\pm 25^{\circ}$. This corresponds to the edges of the molecular ring and could be potentially explained by the 3D distribution of CR sources, or yet by the asymmetry expected in the case of 3D interstellar radiation field.

5.2 Spectral modulation parameters

The best-fit spectra from Run5 are shown in figure 13, both for the full RoI considered in this work as well as for a smaller central region that covers the Galactic bulge. In both regions, the measured spectrum is well reproduced. The error bars on the ICS, gas ring II, and gas ring III components are very small, which is a consequence of the large sky region that they cover. The same is true for the combined emission of the 776 point sources in our analysis, although some $20\%$ uncertainties becomes visible at high energies in the central region. Gas ring I, however, has sizable uncertainties, as it only covers a small sky region and overlaps along the line-of-sight with the other gas rings. We remark that our results are relatively stable w.r.t. choices of the adopted smoothing parameters, as discussed in appendix A.3.

We find that the inner gas ring II has a somewhat harder spectrum at higher energies than the outer gas ring III (the effect is not visible is the poorly constrained ring I). This hardening can also be observed in figure 14, where we show the strength of the spectral modulation parameters for the diffuse emission components in Run5. There is no strong drag for gas ring I and the ICS components, which indicates self-consistency in the results. However, the spectra of the gas rings II and III harden by more than 3$\sigma$. Although (part of) the indicated spectral change could be also caused by some unmodeled extended sources, this hardening of the spectrum towards the inner Galaxy could be a real feature of the interaction of CR protons with the gas, and point towards a concentration of CR accelerators or specific propagation effects in that region [12]. Although a more systematic analysis would be required in this respect, we find results qualitatively in agreement with Refs. [66, 25].

5.3 The Galactic Center Excess

Lastly, the 511 keV template, which we use to model the Fermi Galactic center excess (GCE), exhibits a spectrum that peaks, as usual, at 1–3 GeV, and falls quickly off at lower and higher energies. The error bars are here at the level of $10-20\%$ at the peak. As can be seen in the right panel of figure 13, at energies above 10 GeV, the template does not absorb any significant emission, which shows that the high-energy emission found in previous studies [6, 91] is here absorbed by other components, in particular the Fermi bubbles. Here, the Fermi bubbles absorb the high-energy emission as a consequence of not having imposed a uniform brightness bubble template. Very marginal negative residuals remain at energies between 10 and 50 GeV, most likely because of the spectral constraint on the Fermi bubbles, which could be reduced by relaxing the constraints on that component.

The only change in our model between Run4 and Run5 is the addition of the GCE component. We can in theory determine its significance straightforwardly with the standard likelihood ratio test; the difference in the poisson likelihood between Run4 and Run5 is approximately $208$. However, the difference in the degrees of freedom between the two models cannot be read off from table 1, because the difference between Run4 and Run5 is smaller than the uncertainty in the number of effective degrees of freedom (see section 6.1).

We therefore performed a test where the difference in degrees of freedom is obvious: instead of allowing some small spatial modulation in the GCE component, we fix the spatial template completely, but allow the spectral modulation to be completely free. We kept all other parameters the same as in Run5. For this run, the difference in the degrees of freedom between Run4 and Run5 is equal to $25$, or the number of spectral bins. The Poisson likelihood decreases such that the signifcance for the GCE is about $12\sigma$. The overall fit is for this run is quantitatively very similar to our original version of Run5, as expected. We can therefore conclude that the GCE component is highly significant.

6 Discussion

6.1 The goodness-of-fit

We are interested in estimating the goodness-of-fit for our final Run5 model. The traditional Pearson’s chi-squared goodness-of-fit test cannot be applied to the current problem for two reasons. First, the data is Poisson and not normally distributed. Second, given that we introduced a number of penalization constraints on the model parameters, as well as smoothing constraints, it is not clear what the effective number of free model parameters really is.

We use an approach that fully takes into account the Poisson nature of the problem as well as the unclear effective number of degrees of freedom (DOF), and which is based on the mock results obtained in appendix A.3. First, we consider here only the Poisson part of the likelihood function, and define $\chi^{2}\equiv-2\ln\mathcal{L}_{P}$. We assume that the $\chi^{2}$ function is chi-squared distributed under repeated data realizations (this is true to very good approximation, since for the very high number of DOF the central limit theorem dominates the shape of the distribution). We estimate the effective number of DOF from fits to mock data that are based on the best-fit models for each of the five runs in table 1. It can be estimated to be $k\approx\langle-2\ln\mathcal{L}_{P}\rangle_{\rm mock}$; its associated uncertainty is simply the square root of the variance of the chi-square distribution, $\sqrt{2k}$, about $\sim 1100$ for our five runs. The number of effective DOF along with the final Poisson likelihood and regularization terms, $-2\ln\mathcal{L}_{P}$ and $-2\ln\mathcal{L}_{R}$, can be found in table 1 for each fit.

Naively, the number of model parameters is just given by the number of free parameters in the fit, as listed in table 1. The number of independent data points would normally be simply equal to the number of pixels times the number of energy bins,

[TABLE]

For the current analysis, we have $N_{\text{data}}=360\times 81\times 25=729000$. However, these estimates are only valid if we are in the Gaussian regime in all data bins, and if all model parameters are free to vary and can independently and without degeneracies improve the fit to the data. However, these conditions are not satisfied in our current setup. We therefore estimate the effective number of data bins, model parameters, and DOF in a way that makes these estimates useful for goodness-of-fit calculations.

The effective number of data bins is calculated by averaging the $-2\ln\mathcal{L}_{P}$ over mock realizations of the best fit model, without refitting model parameters. This can be written as

[TABLE]

where we explicitly show that the average can be written as a sum over all data bins. We furthermore defined a function $d(\mu)$, which can be interpreted as the effective number of “degrees of freedom” that correspond to a data bin with expectation value $\mu$, assuming Poisson noise. The sum of this function over all data bins gives the total number of effective data bins, $N_{\text{data}}^{\text{eff}}$. We show that function in figure 15. For $\mu\ll 1$, one can show that it follows $2\mu\ln(1/\mu)$; for $\mu\gg 1$ it converges to one. We use a tabulated version of the function $d(\mu)$ to calculate $N_{\text{data}}^{\text{eff}}$ for each of the five runs. The results are similar for all runs, with values around 620000, and shown in table 1.

The number of effective model parameters can now be estimated in a similar way. We define it as

[TABLE]

where $\boldsymbol{\hat{\theta}}$ are the maximum-likelihood estimators for the model parameters, and different for each data set the average is taken over. Note that $\left\langle-2\ln\mathcal{L}_{P}(\boldsymbol{\hat{\theta}})\right\rangle_{\mathcal{D}(\boldsymbol{\theta})}$ is simply equal to the DOF, $k$, which we define above. For each of the five runs, instead of performing the average, we in practice use only one mock realization to obtain $N_{\text{param}}^{\text{eff}}$ (the expected standard deviation is $\sim 1100$ and not relevant for our discussion). The resulting values are always significantly smaller than the nominal number of free model parameters, $N_{\text{param}}$, and summarized in table 1.

We estimate the fit-quality of the five runs by using some ‘model fidelity’ parameter,

[TABLE]

It can be interpreted as the expected statistical significance (in standard deviations) of any additional single-parameter component that would bring the model at hand in perfect agreement with the data. Since the expected significance is lower for a multi-parameter component, our model fidelity parameter $\mathcal{F}$ provides an upper limit on the significance at which additional unmodeled components could be identified in the data by improvement of the model. We consider this quantity to be more useful than the standard $p$-values, which remain extremely small for all models considered in this work ($p\lll 10^{-315}$ for Run1 and $p=10^{-77}$ for Run5).

The model fidelities $\mathcal{F}$ for Run1–5 are summarized in table 1. We find that the most drastic improvement is between Run1 and Run2, where the fit fidelity drops from 627 to 164. Subsequent refinements of the model appear to lead only to a relatively mild improvement of the fit quality. The still relatively large value for Run5, 144, is likely dominated by residuals around bright point sources and boundary effects related to our treatment of the PSF. We stress that the statistically significance for including additional model components with many parameters cannot be directly read off by comparing the improvement in $\mathcal{F}$ alone.

6.2 Future directions and potential applications

In this paper we present a novel method for fitting the gamma-ray sky that can help us in deepening our understanding of gamma-ray emission mechanisms. We discuss here a number of possible or planned extensions of the current framework, and some of the science questions that it can be used for to address.

On the scientific side, we are interested in full parameter scans over GALPROP/DRAGON predictions for cosmic-ray diffusion and gamma-ray emission. To this end, we would leave both the propagation model parameters as well as the relevant template modulation parameters free to vary. One could then study constraints on GALPROP/DRAGON parameters (like the halo height, the diffusion coefficient, source distributions, etc) while simultaneously accounting for, e.g., hard-to-model components like the Fermi bubbles. However, such a simultaneous fit is currently hindered by the long run-times of GALPROP/DRAGON when generating high-resolution gamma-ray predictions for large regions of the sky.

A combined analysis could be used, e.g., to study potential cosmic-ray gradient, or substructure in the inverse Compton emission. This is for instance relevant to understand whether there is really a hardening of the proton spectrum suggested by the data, as previously claimed [66, 25]. One particularly interesting aspect is the potential observability of spiral arm structures in the gamma-ray emission of the Milky Way. It is well known that the distribution of CR sources as modeled in numerical codes might not be the most realistic one. For example, we do know that CR sources are expected to be distributed in spiral arms. The impact of spiral arms modeling has been discussed in the literature in the context of CR observables [92, 93, 94]. Also, the spiral arm dynamic can have an impact of secondary-to-primary ratios [95]. However, so far, the implications for gamma rays observables remain largely unexplored, with a pioneering work [96] showing that the spiral arm dynamic might have a clear signature on gamma-ray pion spectrum. Fitting both propagation model parameters and nuisance parameters simultaneously would be likely important to make sure that a detection of spiral arm structures is not just an effect of mismodeling of gas.

Finally, an obvious application is the characterization of Fermi bubbles at low latitudes and of possible degeneracies with the Galactic center excess. As we saw above, we are able to reconstruct the bubbles template given their spectral distribution, without including a specific spatial template. A more detailed analysis of the degeneracy between the spectral and spatial characteristics of the bubbles and the GCE would enable us to characterize very large scale emission down to the center of the Galaxy and might provide a way to disentangle various emission processes.

On the technical side, a full Bayesian extension of the current analysis would be interesting and desirable. To this end, Hamiltonian Monte Carlos are particularly promising, since they make use of the available gradient information. In this framework, our sampling from the inverse Fisher information to estimate errors could be replaced by directly drawing from the proper posterior distributions of the model parameters. However, a sufficiently fine sampling of the posterior would require speeding up the code significantly and/or using variational inference for the posterior distributions. Another interesting technical extension would be the parallelization of the code, which is, as mentioned above, currently hindered by the difficulties of parallelizing sparse matrix multiplications.

7 Conclusions

We presented a new hybrid approach for studying, modeling and decomposing diffuse gamma-ray emission, which we dubbed SkyFACT (Sky Factorization with Adaptive Constrained Templates). Our approach combines methods of image reconstruction and adaptive spatio-spectral template regression in one coherent framework, based on penalized Poisson likelihood regression.

We discussed in detail the implementation of our approach, as well as various solutions to technical challenges related to the high dimensionality of the optimization problem. In particular, we showed how the L-BFGS-B algorithm can be used to find unique solutions in high-dimensional parameter spaces. The largest example in this work has $1\times 10^{5}$ parameters. We showed that the optimization problem is, for scenarios considered here (and for sufficiently soft smoothing), convex, and hence has no non-global minima. With conventional desktop computers, convergence can be reached in a few dozens of core hours. We demonstrated how Cholesky decomposition of the Fisher matrix, which we use here as approximation to the inverse covariance matrix, can be used to sample parameter errors. From these, uncertainties and correlations of various model components can be inferred. This is an important step towards going beyond the typically adopted strategy of ‘bracketing uncertainties’ with discrete model choices, although model assumptions still have to be made.

We applied our approach to the gamma-ray emission from the inner Galactic disk, $|\ell|<90^{\circ}$ and $|b|<90^{\circ}$, using 7.6 years of Fermi-LAT data binned in $0.5^{\circ}$ pixels and covering the energy range 0.34–228.65 GeV. We presented a series of five increasingly complex models that show the challenges and typical residuals for conventional fits to the data, and the reduction of residuals when realistic nuisance parameters are included in the fit. The final, most complex model includes three gas rings, an ICS template, the IGRB, a spectral template for the Fermi bubbles, a spatial template for the Fermi Galactic center excess (based on INTEGRAL measurements of the 511 keV line), as well as Fermi-LAT extended and point sources. As a result, we can remove most of the residuals over much of the RoI and across all energies in the analysis.

We also performed a series of tests with synthetic data to explore the performance of the regularization used in our fits; see appendix A. In general, we found that moderate regularization in the form of either modulation or smoothing can effectively recover model components without signficant overfitting.

The main physics results of our analysis will be presented elsewhere. However, we summarize a few findings:

• •

We do not include dark gas corrections in the gas maps, but instead recover the missing dark gas from the nuisance parameters in the fit. In general, the spatial nuisance parameters rescale the a priori gas and ICS templates that we use from DRAGON by less than $\pm 40\%$, which is enough to provide good fits to the data, and within the uncertainties of these components.

• •

We furthermore do not include a spatial template for the Fermi bubbles, but instead recover their flux at high and lower Galactic latitudes. We find that the bottom of the Fermi bubbles is not coincident with the Galactic center, but lies at around $\ell\approx-2^{\circ}$, which was also observed by [49]. However, despite being bright, the emission is not very significant, and we cannot exclude that this is an extended unmodeled component of the foreground.

• •

For the Galactic center excess we use a spatial template based on the morphology of the Galactic center and bulge emission of the 511 keV line measured by INTEGRAL. We find that the template is capable of removing the inner Fermi Galactic center excess and is detected at about the $12\sigma$ level, and that it picks up a spectrum compatible with the combined emission of (young or millisecond) pulsars.

Acknowledgements.

We are indebted to D. Gaggero for providing the package GammaSKY and helping with the GDE templates. F.C. thanks E. Charles for valuable discussion. C.W. and E.S. acknowledge funding from an NWO Vidi fellowship.

Appendix A Tests with synthetic data

A.1 Image reconstruction, MEM and smoothing

We first test both the effect of smoothing and MEM on image reconstruction for varying signal-to-background ratios and number of counts. We build a mock data set that consists of two spatial templates: a background that varies exponentially with latitude and is constant in longitude, and a set of three rings with flat, uniform intensity. No energy dependence is included in the data or fits. The ring and background templates are normalized such that the signal-to-background ratio is roughly equal to $10^{-2}$ at the top of the map ($b=20.25^{\circ}$), $10^{-1}$ at $b=10.125^{\circ}$, $1$ at $b=0^{\circ}$, $10$ at $b=-10.125^{\circ}$, and $10^{2}$ at the bottom ($b=-20.25^{\circ}$).

Poisson noise is added and the data is divided into three equal-longitude sections of $60^{\circ}$. The exposure in each section varies such that the total number of counts is different in each section. There are approximately $10^{5}$ in the left section, $10^{6}$, in the middle, and $10^{7}$ in the right section. In the real Fermi-LAT data we use in our fits, there are roughly $5\times 10^{6}$ photons in the lowest energy bin. At the top of the map, where the signal-to-background is small, the average number of background counts is $10^{2}$ in the left section, $10^{3}$ in the middle, and $10^{4}$ in the right section. Correspondingly, at the bottom where the signal-to-background is large, the average number of background counts is 0.01, 0.1, and 1 in the left, middle, and right sections. We show the counts map used for this analysis in figure 16.

To perform fits, we fix the background component completely. For the ring component, we initialize the template to be completely flat, and set the template modulation parameter to be essentially free ($\lambda=0$). We perform four fits with increasing amounts of smoothing applied to the ring template ($\eta=0$, $1$, $4$, and $100$).

We show the results of these runs in figure 17. The top-left image shows the recovered “flux” map where no smoothing was used in the template reconstruction. It is obvious even by eye that some amount of background is being recovered in the ring template; the histograms of residuals for this run also indicate this. In general, the rings are more accurately recovered in regions of high signal-to-background and large number of counts, as seen in the bottoms of right-most sections of the “flux” maps. For regions with fewer photons, and low signal-to-background, the rings are less well-recovered, seen in the tops of the middle and left sections of the flux maps. However, moderate amounts of smoothing can help to recover the rings and suppress too much fitting of the Poisson flucuations; compare, for example, the left sections of the middle two flux maps. If the smoothing parameter is too large, the ring shape is washed out, as seen in bottom flux map, especially for low photon counts. The strong positive and negative residuals seen in the histogram for this run with the strongest smoothing are mostly the result of the edge of ring getting smeared out: there are strong positive residuals on the inside edge of the ring and strong negative residuals on the outside edge.

We perform a second set of runs without smoothing, but with varying template modulation parameters, to illustrate how MEM regularization alone influences image reconstruction. Again, we fix the background component and initialize the ring template to be flat. We vary the template modulation parameter in each fit ($\lambda=1,10^{-2},10^{-4}$). If $\lambda>1$, there is insufficient variation permitted in the template parameters to recover the ring flux.

These results are shown in figure 18. In the top panel, where the template modulation is strongest, there are still significant remaining residuals. Essentially, this template modulation is almost too strong, especially in the low-photon section, to allow a flat template to completely reconstruct the ring flux. For weaker template modulation, there are very few positive residuals, indicating that the ring flux is fully recovered (this is evident from the residual maps, not shown, as well.) The negative residuals that remain are the result of the flat template not being able to go completely to zero; the rescaling for most of those pixels is at their smallest allowed value ($10^{-10}$; this is also true for the runs with smoothing). However, also evident is some amount of Poisson noise being picked up by this template, seen especially in the high-background-count upper regions of the flux maps. As seen in the previous set of runs, even a weak smoothing constraint applied here would drastically reduce the amount of fitted Poisson noise, especially in low-photon regions.

A.2 Component separation and Poisson-induced bias with two components

For the next set of mock data runs, we wish to test how fits perform at separating model components (both templates and spectra) with different values for smoothing and template modulation parameters. For these tests, we build a mock data set with two components: the gas and ICS components described in section 3. Here we also include energy dependencies for both components, using the spectra described in section 3 for the gas and ICS components. The energy range we fit over is the same as for the real LAT data: $0.34-228.65$ GeV. We choose these two components because they are roughly similar in overall brightness in our RoI, have similar spectral shapes at high energies, and share some morphological characteristics. The diffuse emission from the gas and ICS also dominate the real LAT data we consider. We show the mock data counts map and the initial spectral energy distribution in figure 19.

For all of the following fits, we allow the overall normalization of each of the templates to freely vary ($\lambda^{\prime\prime}=0$). The templates and spectra are initialized to their original components described in section 3. All runs were performed with MEM regularization. We show the hyper-parameter values for each fit in table 2.

We first show the results of a run with moderate template and spectra modulation parameters, and no smoothing ($\lambda,\lambda^{\prime}=1$, $\eta,\eta^{\prime}=0$). As seen in figure 20, the fit fails to properly recover the two components, instead separating into a “low-energy” component and a “high-energy” component. If the modulation is set to be weaker, the same behavior is seen (the error bars on the spectra become larger). While a histogram of residuals and map of residuals summed over all energies indicates no substantial deviations, there is clear structure remaining in at least the higher energy residual map, all shown in figure 21. Additionally, as shown in figure 22, there are large deviations from the original spectral and spatial components. While this is formally not a bad fit (its likelihood is similar to those of the following runs), it is obvious from the figures that it is an unphysical result.

Next, we show the results of a run with stronger template and spectra modulation parameters, still with no smoothing ($\lambda,\lambda^{\prime}=4$, $\eta,\eta^{\prime}=0$); these regularization values are near those considered in fits to real data. In this case, the original model components are recovered well, as shown in figure 23. We show the residuals for this fit in figure 24. From the rescaling parameters shown in figure 25, we see very little deviation from the original spectra, and small deviations in the template rescaling

Strong smoothing can act to effectively constrain the allowed variation in the templates, even if the modulation parameters are relatively weak. In figure 26, we show the results from a fit with the same modulation parameters as in figure 20, but with strong smoothing applied to both templates (but still no smoothing on the spectra; $\lambda,\lambda^{\prime}=1$, $\eta=0.1$, $\eta^{\prime}=0$). There are no strong residuals obvious in the maps or histograms in figure 27. As shown in figure 28, there are essentially no deviations from the original spectra or spatial templates.

A.3 Tests with Run5 mock data

For our final set of tests, we build a mock data set that is composed of the best-fit model components to run5, described in section 5. We show the counts map and spectra used for this set of tests in figure 29.

We first fit the mock data set, with the same model components as in Run5, but starting from the best-fit parameters for Run5 (see section  5). As shown in in the top left panel of figure 30, the spectral components are almost perfectly recovered, to within $\sim 30\%$ on average (except for a few obvious energy bins in the 511 keV spectrum, which is not smoothed). We show histograms of the residuals for this mock data set after fitting in the top left panel of figure 31; this can be compared to the original fit for Run5 in figure 7. The Poisson tail is clearly seen in the highest energy bin.

To test the robustness of our results for Run5 on real LAT data, we vary the template modulation and smoothing parameters for this mock data set. In our first test, we set the smoothing on all components to zero. The results of these tests are shown in the top-right panels of figures 30 and 31. The most obvious feature in the spectra plot is the deviation in the ICS spectrum at one high energy bin. The mockdata set is Poisson-dominated at this energy, and there are few photons actually being fit. Additionally, the spatial and spectral modulation parameters are somewhat weaker for the ICS component than the other components, especially the gas, which, in the absence of of smoothing, means that the ICS component is permitted more variation. This deviation is clearly unphysical, and demonstrates the need for moderate regularization in low photon regimes, as was also seen in the previous two sections. Interestingly, the overall deviation from the data is also quite high, at the $\sim 2\sigma$ level for three energy bins, indicating that these energies are not well fit. The overall likelihood is not significantly different, however, from any of the other mock data runs in this section.

We also perform fits where the spatial modulation parameters for the ICS and three gas components are set to be weaker and stronger by a factor of 4 relative to those in table 1. The other modulation parameters and all smoothing parameters remain unchanged from the original Run5 set-up. The results of these tests are shown in the bottom panels of figures 30 and 31. Again, in the case of the weaker modulation (bottom left), the ICS component is considerably suppressed while the outer gas ring III is enhanced relative to the original best-fit parameters; the 511 keV template is also somewhat enhanced. This illustrates how weakly constrained the ICS component is relative to the other components. The residual histograms in figure 31 all show reasonable overall fits.