Dissecting the inner Galaxy with γ-ray pixel count statistics
LLAPTH-008/21, TTK-21-06
Dissecting the inner Galaxy with γ -ray pixel count statistics F. Calore, ∗ F. Donato,
2, 3, † and S. Manconi ‡ Univ. Grenoble Alpes, USMB, CNRS, LAPTh, F-74940 Annecy, France Dipartimento di Fisica, Universit`a di Torino, via P. Giuria, 1, I-10125 Torino, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria, 1, I-10125 Torino, Italy Institute for Theoretical Particle Physics and Cosmology,RWTH Aachen University, Sommerfeldstr. 16, 52056 Aachen, Germany
We combine adaptive template fitting and pixel count statistics in order to assess the nature ofthe Galactic center excess in
Fermi -LAT data. We reconstruct the flux distribution of point sourcesin the inner Galaxy well below the
Fermi -LAT detection threshold, and measure their radial andlongitudinal profiles. Point sources and diffuse emission from the Galactic bulge each contributes O (10%) of the total emission therein, disclosing a sub-threshold point-source contribution to theGalactic center excess. Introduction
The Galactic center excess (GCE)shows up as an unexpected γ -ray component in thedata of the Large Area Telescope (LAT), aboard the Fermi satellite, at GeV energies, from the inner degreesof the Galaxy [1–5]. Despite the great interest raisedby the GCE discovery, its nature is still unknown. Whilethe GCE morphology is consistent with a Navarro, Frenkand White (NFW) profile [6] for annihilating particledark matter (DM) [1, 4, 7–9], it could also be due toa population of millisecond pulsars, as proposed by [10].Stellar distributions were used as tracers of point sources(PS) emitting below threshold, and turned out to matchthe morphological features of GCE photons better thanDM-inspired templates [11–13]. Complementary studiesof photon-count statistics revealed that the GCE can beentirely due to a population of PS [14, 15]. Recently, theDM interpretation was brought back by [16], althoughhampered by systematics affecting photon-count statis-tical methods [17–21]. As a conclusive probe of the PSnature of the GCE, a fully multiwavelength approach hasbeen proposed, from radio to gravitational wave observa-tions [22–24].A major limit to all these studies is the modeling of theGalactic diffuse foreground, and the impact of residualmis-modeled emission on the results’ robustness. As fortemplate fitting methods, the analysis of the diffuse emis-sion has been recently approached with the skyFACT al-gorithm, which fits the γ -ray sky by combining meth-ods of image reconstruction and adaptive spatio-spectraltemplate regression [25]. The skyFACT method has beentested in the Inner Galaxy (IG) region, and probed to beefficient in the removal of most residual emission for a ro-bust assessment of the GCE properties [11, 25]. Anothersource of uncertainty is the contribution of sub-thresholdPSs. Photon-count statistical methods can discriminatephotons from γ -ray sources based on their statisticalproperties [26]. In particular, the 1-point probability dis-tribution function method [27] ( ) fits the contribu- tion of diffuse and PS components to the γ -ray 1-pointfluctuations histogram. Employing on Fermi -LATdata, it was possible to measure the PS count distributionper unit flux, dN/dS , below the LAT detection thresholdat high latitudes [27–29], and to set competitive boundson DM [30].The scope of this
Letter is to apply the methodto
Fermi -LAT data from the IG to understand the role offaint PS to the GCE, while minimizing the mis-modellingof diffuse emission components. To this end, we adopta hybrid approach which combines, for the first time,adaptive template fitting methods as implemented in skyFACT , and techniques.
Rationale, data and methodology.
We followa two-step procedure: First, we fit γ -ray data with skyFACT in order to build a model for the emission inthe region of interest (ROI), maximally reducing residu-als found to bias photon-count statistical methods [19].Secondly, we run fits with skyFACT -optimized dif-fuse models as input, and assess the role of PS to theGCE.We analyze 639 weeks of P8R3 ULTRACLEANVETO
Fermi -LAT data [31] until 2020-08-27. For the skyFACT fit, we consider an ROI of 40 ◦ × ◦ around theGC [32], and the 0 . −
300 GeV energy range. We closelyfollow [11] and update the analysis for the increased dataset and 4FGL catalog [33]. The emission model includes γ rays from inverse Compton scattering, π decay, 4FGLpoint-like and extended sources, the Fermi bubbles, theisotropic γ -ray background (IGRB), and the GCE. Forthe latter, we consider a template for the Galactic bulgeemission as in [11], and one for a generalized NFW DMdistribution with slope 1.26 (NFW126) [3, 4]. We referto [34] for more details.We operate the analysis in the energy range2 − ◦ × ◦ , IG ROI hereafter. We cut at latitudes a r X i v : . [ a s t r o - ph . H E ] F e b | b | > . ◦ or 2 ◦ to check the stability of results.The -fit model components are: An IGRB tem-plate (free normalization), a diffuse emission template(free normalization), and an isotropic PS (IPS)[35] pop-ulation with dN/dS defined by a multiple broken powerlaw:d N d S = A S · (cid:16) SS (cid:17) − n S > S b1 ; (cid:16) S b1 S (cid:17) − n + n (cid:16) SS (cid:17) − n S b2 < S ≤ S b1 ;... (cid:16) S b1 S (cid:17) − n + n (cid:16) S b2 S (cid:17) − n + n · · · (cid:16) SS (cid:17) − n N b+1 S ≤ S b N b . (1)The free parameters are A S , the flux break positions,and the broken power-law indices, n i [34]. The IPS dN/dS measured by the fit should recover the dN/dS of Fermi -LAT detected PS in the bright regimewhile pushing the PS detection threshold down to lowerfluxes [27, 28]. Our goal being to quantify the role ofPS to the GCE within the , we add a
GCE smoothtemplate (free normalization) in the fit. As a base-line, we use the best-fit skyFACT bulge template as GCEin the fit ( ), and we define the sF-B diffusemodel as the sum of best-fit inverse Compton, π decay, Fermi bubbles, and extended sources, thus subtractingthe bulge emission.On the one hand, the use of skyFACT best-fit diffusemodel guarantees a robust characterization of GCE spec-trum and morphology against systematics related to themis-modeling of the diffuse emission [3, 19], resolvingover/under-subtraction issues by including a large num-ber of nuisance parameters. The limitations of such asystematic uncertainty are indeed also relevant for thereconstruction of faint PS with methods [34]. Onthe other hand, the skyFACT optimization procedure mit-igates possible systematics related to the mis-modelingof unaccounted components [17], by allowing spatial re-modulation in the fit templates. Also, we stress that thisis the first time the stellar distribution in the Galacticbulge as tracer of GCE photons is used in pixel countstatistical analyses.Besides the bulge, we also consider NFW126 as GCEin the analysis ( ). In this case,we construct the corresponding skyFACT -optimized dif-fuse model ( sF-NFW126 ) from the skyFACT run adoptingNFW126 as GCE, in analogy with the sF-B model. Sucha procedure guarantees maximal consistency betweenGCE and diffuse models adopted as input in the .Finally, to bracket the uncertainties related to the opti-mization of the diffuse model, we also build a skyFACT -optimized diffuse template from the skyFACT run not in-cluding any GCE additional template ( sF-noGCE ). Results.
Our updated analysis of
Fermi -LAT data with skyFACT confirms previous findings from [11–13]. Abulge distribution for GCE photons is strongly preferredby data on top of the NFW126-only model ( ∼ σ ), andthere is no evidence for an additional NFW126 contribu-tion on top of the bulge-only model ( ∼ σ ), cf. [34]. Thisimplies that the model maximally reducing the residualsis the skyFACT best-fit of the run with the bulge.We then use skyFACT -optimized diffuse and GCE tem-plates as input for fits, testing 0.5 ◦ and 2 ◦ lati-tude cuts. Our results are summarized in Fig. 1, wherewe show the best-fit dN/dS for the IPS in the IG ROIfor several fit configurations. First, we notice thatwhatever GCE template is added to the fit compo-nents (bulge or NFW126), its normalization never con-verges toward the lower bound of its prior interval, re-gardless of the skyFACT diffuse template adopted. Thesame is valid for the IPS normalization. In all fit se-tups shown, an IPS population is recovered below theLAT flux threshold. The reconstructed IPS dN/dS is sta-ble against systematics related to the choice of skyFACT -optimized diffuse template, and latitude cut. Moreover,it does not present any spurious effect at the Fermi -LATthreshold ( ∼ − ph cm − s − ), and IPS are resolveddown to ∼ − ph cm − s − for | b | > . ◦ , dependingon the GCE modeling. This holds true even when noGCE template is included neither in the skyFACT fit norin the one, contrary to what happens using non-optimized diffuse models [17, 19]. We therefore demon-strate, also in the context of methods, that reduc-ing large-scale residuals from mis-modeling of the diffuseemission improves the reconstruction of PS dN/dS (seealso [34]). The reconstructed dN/dS has a normaliza-tion decreasing with increasing latitude cuts, suggestingthat PS are more numerous towards the very GC. Whenan NFW126 template is included in the fit, theIPS dN/dS is compatible with the case.In both cases, the second break in the dN/dS – inaddition to the one set in the bright regime – is re-covered close to the LAT flux threshold. Instead, the reconstructs PS down to lower fluxes, regardlessof sF-noGCE or sF-B diffuse models. For these setups, thesecond flux break is found at ∼ · − (8 · − ) phcm − s − for | b | > . ◦ (2 ◦ ). Going from | b | > . ◦ to | b | > ◦ , the dN/dS is resolved down to even lower fluxes.A posteriori, we associate such a better sensitivity to IPSto the ability of the fitted diffuse components to furtherreduce fit residuals.We quantify now the evidence for models with an ad-ditional GCE template. To this end, we compare theglobal evidence, ln Z , for the , and setups, with different skyFACT diffusemodel inputs. For each model combination, we com-pute the Bayes factor between model i and j , B ij =exp(ln Z i − ln Z j ), and assess the strength of evidenceof model i with respect to model j . Our results are pre-sented in Tab. I. Regardless of the skyFACT -optimized − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > . ◦ − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > ◦ FIG. 1.
IPS source count distribution in the IG ROI from the fit for | b | > . ◦ (left) or 2 ◦ (right). Solid (dashed) linescorrespond to sF-noGCE ( sF-B ) diffuse template. The black line illustrates the case. The blue (red) line refers to ( ) case. The colored areas correspond to 1 σ uncertainty bands. The black (gray) points represent thecount distribution of 4FGL sources (without any analysis flag, see [33]).TABLE I. Results for the analysis of the IG
LAT data. First four columns: setup of the analysis and latitude mask ofthe IG. The ln( Z ) is the nested sampling global log-evidence extracted from Multinest [27]. Last two columns: flux percentageof different model components with respect to the total emission in the ROI (for
S < − ph cm − s − , see [34]), andnormalization of GCE template in the . Flux percentage always sum to unity within errors. Description setup skyFACT diffuse | b | cut [ ◦ ] ln( Z ) Point sources/diffuse/GCE % A B / NFW126
No GCE (both) − / / − -Bulge ( only) − / / . ± . only) − / / . +0 . − . Bulge ( skyFACT only) − / / − -Bulge (both) − / /
10 1 . ± . − / / . ± . − / / − -Bulge ( only) − / / . ± . only) − / / . ± . skyFACT only) − / / − -Bulge (both) − / / . ± . − / / . ± . diffuse template adopted, data always more strongly sup-port models which include an additional smooth templatefor the bulge with respect to models without GCE in the skyFACT and/or fits (ln B ij > and modelswith an additional smooth NFW126 component in the skyFACT and/or fits (ln B ij > sF-B , and | b | > ◦ ,the evidence for an additional bulge template ( ),with respect to is ln B ∼
95. Moreover,in this case the normalization of the bulge template is A B = 1 . ± .
1, supporting the consistency between GCEand diffuse model adopted. This evidence is as strongalso for | b | > . ◦ , ln B ∼ sF-noGCE diffuse model in the fit including the bulge ( ), we findcomparable evidence to the , sF-B setup. In-deed, skyFACT is able to re-absorb part of the photonsfrom the bulge by re-modulating (spatially) other dif-fuse templates, and so, partially reduces the residualsalso in the sF-noGCE case. This is perfectly consistentwith the fact that A B = 0 . ± .
1. Models with PS and asmooth bulge component are therefore strongly preferredby data, regardless of the optimized diffuse model em-ployed. On the contrary, the evidence for an additionalsmooth NFW126 template with respect to models with-out GCE in the skyFACT fit and/or fits dependson the choice of the skyFACT -optimized diffuse templateadopted, as well as on the latitude cut.We have tested our results against a number of sys-tematic effects, which are detailed in [34].
GLAT [deg] d N / d Ω [ d e g ] ExtragalacticOuter Galaxy4FGL1pPDF − − − − GLON [deg] d N / d Ω [ d e g ] ExtragalacticOuter Galaxy4FGL1pPDF
FIG. 2.
Latitude (left) and longitude (right) source density dN/d Ω profiles , as reconstructed by the fit using the sF-B diffuse model. We also display source density profiles for 4FGL sources (black points), and average source densities in theOG and EG ROIs. The flux percentages reported in Tab. I illustrate that fits to
Fermi -LAT data find non-null (and evencomparable) emission from both the IPS population and the smooth GCE template, in most cases each contribut-ing about 10% of the total emission in the ROI. Since4FGL sources (2 ◦ cut, without analysis flag, see Fig. 1)account for 7% (10% including flagged sources) of thetotal IG emission, the remaining flux comes from sub-threshold IPS. We have verified [34] that our results arenot driven by PS in the ultra-faint regime [18], wherethe sensitivity of the method drops (as quantifiedby the magnitude of uncertainty bands in Fig. 1), andan IPS population may become degenerate with a trulydiffuse emission.We also measure the IPS dN/dS in two control regions:The outer Galaxy (OG, | b | < ◦ , 60 ◦ < | l | < ◦ ) andthe extragalactic region (EG, | b | > ◦ , | l | > ◦ ). Thereconstructed dN/dS in both OG and EG ROIs does notpresent any spurious threshold effect and can identifyIPS down to the statistical limit of the method around ∼ − ph cm − s − [27]. We compute the source den-sity dN/d Ω in the flux interval [10 − − − ] ph cm − s − , finding ∼ . in the OG, and ∼ . in the EG, see Fig. 2.Since the spatial distribution of PS is isotropic by con-struction, we test the PS spatial behavior by dissect-ing the IG ROI into three concentric annuli, maskedfor latitudes | b | < . ◦ . We extract the dN/dS sepa-rately in each ring, and integrate it over the flux interval[10 − − − ] ph cm − s − . The result is reported inFig. 2 as a function of latitude, for our baseline , sF-B setup. We observe a decreasing trend of the dN/d Ωin the IG with latitude. Also, the dN/d
Ω in the inner-most ring is about a factor of three higher than 4FGLsources, as well as than in OG and EG. For the most external latitude bin, the source density is instead com-parable with the catalog, OG and EG ones. This cor-roborates the evidence that the IG PS population is notpurely isotropic nor extragalactic in origin, but rather itpeaks towards the GC. Similarly, we build the longitudeprofile of IG PS, Fig. 2. The dN/dS has been fitted in6 longitude slices from the GC bound at | l | = 6 ◦ , ◦ and 20 ◦ . The derived dN/d Ω shows again a distributionpeaked around the GC, and compatible with OG (andpartially with 4FGL and EG) sources only in the mostexternal longitude interval. This result adds a piece ofevidence that the GCE is contributed by faint PS clus-tered at the GC, supporting their
Galactic origin.
Conclusions.
For the first time, we analyzed theIG
Fermi -LAT sky by means of the photon-countstatistics technique in order to understand the role ofPS to the GCE. To minimize the systematic effects in-herent the modeling of the γ -ray sky, we introduced im-portant methodological novelties. First, we implementedwithin the new, optimized , models for the diffuseemission from skyFACT adaptive template fits, developinga self-consistent procedure which effectively reduces dif-fuse mis-modeling. Secondly, besides PS, in the fitwe included an additional smooth GCE template whichtraces the stellar distribution in the Galactic bulge.The updated skyFACT analysis of the IG confirms thatthe GCE is better described by a bulge template thanan NFW126 model at high significance. Moreover, wefind that the method, supplied with skyFACT dif-fuse emission templates, always recovers an IPS popu-lation well below the
Fermi -LAT flux threshold, downto ∼ − ph cm − s − for | b | > . ◦ . The recon-structed IPS dN/dS is stable against a number of sys-tematics, in particular related to the choice of skyFACT -optimized diffuse template and latitude cut. Regardlessof the skyFACT -optimized diffuse template, data always prefer models which include an additional smooth tem-plate for the bulge with respect to both models withoutGCE and models with an additional NFW126 template,in the skyFACT and/or fits.Our results show that, within the statistical validity ofthe and the setups tested, IPS and diffuse bulgeeach contributes about O (10%) to the IG γ -ray emission.In particular, within our baseline model the foundsthat PS (bulge) contribute 13% (10%) of the total emis-sion of the IG. Subtracting the contribution from cata-loged sources, a non-negligible fraction of the IG emissionis accounted by sub-threshold PS. This further corrobo-rates a possible, at least partial, stellar origin of the GCE.We also verified that this IPS population is not purelyisotropic nor extragalactic in origin, rather it peaks to-wards the very GC. Although the final confirmation ofthe PS nature of the GCE will most likely come frommultiwavelength future observations, we undoubtedly gotone step closer to the understanding of the mysteriousnature of the GCE emission. Acknowledgments.
We very kindly acknowledge thework formerly done by H.S. Zechlin on the code.We warmly thank P. D. Serpico for inspiring discussion.We also thank F. Kahlhoefer, M. Kraemer, P. D. Serpico,and C. Weniger for a careful reading of the manuscriptand for insightful comments. The work of F.D. hasbeen supported by the ”Departments of Excellence 2018- 2022” Grant awarded by the Italian Ministry of Edu-cation, University and Research (MIUR) (L. 232/2016).F.C. acknowledges support by the Programme NationalHautes Energies (PNHE) through the AO INSU 2019,grant “DMSubG”, and the Agence Nationale de laRecherche AAPG2019, project “GECO”. S.M. acknowl-edges computing resources granted by RWTH AachenUniversity under project rwth0578. ∗ [email protected] † [email protected] ‡ [email protected][1] K. N. Abazajian and M. Kaplinghat, Phys. Rev. D86 ,083511 (2012), arXiv:1207.6047 [astro-ph.HE].[2] C. Gordon and O. Macias, Phys. Rev.
D88 , 083521(2013), [Erratum: Phys. Rev.D89,no.4,049901(2014)],arXiv:1306.5725 [astro-ph.HE].[3] F. Calore, I. Cholis, and C. Weniger, JCAP , 038(2015), arXiv:1409.0042 [astro-ph.CO].[4] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden,S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Phys.Dark Univ. , 1 (2016), arXiv:1402.6703 [astro-ph.HE].[5] M. Ajello et al. (Fermi-LAT), Astrophys. J. , 44(2016), arXiv:1511.02938 [astro-ph.HE].[6] J. F. Navarro, C. S. Frenk, and S. D. M. White, Astro-phys. J. , 493 (1997), arXiv:astro-ph/9611107. [7] F. Calore, I. Cholis, C. McCabe, and C. Weniger, Phys.Rev. D91 , 063003 (2015), arXiv:1411.4647 [hep-ph].[8] P. Agrawal, B. Batell, P. J. Fox, and R. Harnik, JCAP , 011 (2015), arXiv:1411.2592 [hep-ph].[9] S. Murgia, Ann. Rev. Nucl. Part. Sci. , 455 (2020).[10] K. N. Abazajian, Journal of Cosmology and Astropar-ticle Physics , 010 (2011), arXiv:1011.4275 [astro-ph.HE].[11] R. Bartels, E. Storm, C. Weniger, and F. Calore, Na-ture Astronomy , 819 (2018), arXiv:1711.04778 [astro-ph.HE].[12] O. Macias, C. Gordon, R. M. Crocker, B. Coleman,D. Paterson, S. Horiuchi, and M. Pohl, Nature Astron-omy , 387 (2018), arXiv:1611.06644 [astro-ph.HE].[13] O. Macias, S. Horiuchi, M. Kaplinghat, C. Gordon,R. M. Crocker, and D. M. Nataf, JCAP , 042 (2019),arXiv:1901.03822 [astro-ph.HE].[14] R. Bartels, S. Krishnamurthy, and C. Weniger, Phys.Rev. Lett. , 051102 (2016), arXiv:1506.05104 [astro-ph.HE].[15] S. K. Lee, M. Lisanti, B. R. Safdi, T. R. Slatyer,and W. Xue, Phys. Rev. Lett. , 051103 (2016),arXiv:1506.05124 [astro-ph.HE].[16] R. K. Leane and T. R. Slatyer, Phys. Rev. Lett. ,241101 (2019), arXiv:1904.08430 [astro-ph.HE].[17] R. K. Leane and T. R. Slatyer, Phys. Rev. Lett. ,121105 (2020), arXiv:2002.12370 [astro-ph.HE].[18] L. J. Chang, S. Mishra-Sharma, M. Lisanti,M. Buschmann, N. L. Rodd, and B. R. Safdi,Phys. Rev. D , 023014 (2020), arXiv:1908.10874[astro-ph.CO].[19] M. Buschmann, N. L. Rodd, B. R. Safdi, L. J. Chang,S. Mishra-Sharma, M. Lisanti, and O. Macias, Phys.Rev. D , 023023 (2020), arXiv:2002.12373 [astro-ph.HE].[20] Y.-M. Zhong, S. D. McDermott, I. Cholis, and P. J. Fox,Phys. Rev. Lett. , 231103 (2020), arXiv:1911.12369[astro-ph.HE].[21] R. K. Leane and T. R. Slatyer, Phys. Rev. D , 063019(2020), arXiv:2002.12371 [astro-ph.HE].[22] F. Calore, M. Di Mauro, F. Donato, J. W. T. Hessels,and C. Weniger, ApJ , 143 (2016), arXiv:1512.06825[astro-ph.HE].[23] F. Calore, T. Regimbau, and P. D. Serpico, PRL ,081103 (2019), arXiv:1812.05094 [astro-ph.HE].[24] J. Berteaud, F. Calore, M. Clavel, P. D. Serpico,G. Dubus, and P.-O. Petrucci, (2020), arXiv:2012.03580[astro-ph.HE].[25] E. Storm, C. Weniger, and F. Calore, JCAP , 022(2017), arXiv:1705.04065 [astro-ph.HE].[26] D. Malyshev and D. W. Hogg, ApJ , 181 (2011),arXiv:1104.0010 [astro-ph.CO].[27] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo,and A. Vittino, Astrophys. J. Suppl. , 18 (2016),arXiv:1512.07190 [astro-ph.HE].[28] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo,and M. Regis, Astrophys. J. Lett. , L31 (2016),arXiv:1605.04256 [astro-ph.HE].[29] S. Manconi, M. Korsmeier, F. Donato, N. Fornengo,M. Regis, and H. Zechlin, Phys. Rev. D , 103026(2020), arXiv:1912.01622 [astro-ph.HE].[30] H.-S. Zechlin, S. Manconi, and F. Donato, Phys. Rev. D , 083022 (2018), arXiv:1710.01506 [astro-ph.HE].[31] Publicly available at https://heasarc.gsfc.nasa.gov/ \FTP/fermi/data/lat/weekly/photon/ . More detailsare found at https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Data/LAT_DP.html .[32] This ROI still allows to discriminate the GCE morphol-ogy without suffering from systematics induced by select-ing too narrow ROIs [13].[33] S. Abdollahi et al. (Fermi-LAT), Astrophys. J. Suppl. , 33 (2020), arXiv:1902.10045 [astro-ph.HE].[34] See Supplemental Material at [URL will be inserted bypublisher]. Sec. I for details on the SkyFACT analysis,Sec. II for details on the 1pPDF analysis, Sec. III-IV forsystematic checks on the diffuse emission template anddN/dS modeling within the 1pPDF, and Sec. V for thedN/dS of the outer and extragalactic regions.[35] As currently implemented, the method does notallow to test spatially-dependent dN/dS .[36] S. Mishra-Sharma, N. L. Rodd, and B. R. Safdi, Astron.J. , 253 (2017), arXiv:1612.03173 [astro-ph.HE]. [37] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt,F. K. Hansen, M. Reinecke, and M. Bartelmann, ApJ , 759 (2005), astro-ph/0409513.[38] F. Feroz, M. P. Hobson, and M. Bridges, MNRAS ,1601 (2009), arXiv:0809.3437.[39] W. A. Rolke, A. M. L´opez, and J. Conrad, Nuclear In-struments and Methods in Physics Research A , 493(2005), physics/0403059.[40] F. Acero et al., ApJS , 26 (2016), arXiv:1602.07246[astro-ph.HE].[41] M. Ackermann et al., ApJ , 86 (2015),arXiv:1410.3696 [astro-ph.HE].[42] L. Maccione, C. Evoli, D. Gaggero, and D. Grasso,“DRAGON: Galactic Cosmic Ray Diffusion Code,”(2011), ascl:1106.011.[43] M. Ackermann et al. (Fermi-LAT), Astrophys. J. , 64(2014), arXiv:1407.7905 [astro-ph.HE].[44] M. Ackermann and others, ApJ , 3 (2012),arXiv:1202.4039 [astro-ph.HE]. Supplemental Material:
Dissecting the inner Galaxy with γ -ray pixel count statistics F.Calore, F.Donato, and S.Manconi
THE
SKYFACT
ANALYSIS
In this section, we discuss in more detail the analysis and results of the
Fermi -LAT γ -ray fit with skyFACT [25].For the skyFACT analysis, we consider all FRONT+BACK events (evtype=3) to maximize the statistics against thelarge number of free parameters in the fit. The data are binned in energy into 30 logarithmically-spaced bins from 0.2to 500 GeV, and spatially into cartesian pixels of size 0.5 ◦ . The fit is performed in the energy range 0 . −
300 GeV,and the main ROI restricted to the inner 40 ◦ × ◦ . This skyFACT ROI is larger than our inner Galaxy (IG)ROI since, for the purpose of template fitting, the ROI must be large enough to be able to correctly assess the GCEmorphology [13]. The model of the γ -ray sky diffuse components and the statistical analysis closely follow [11]. Themain novelties here are: (i) The increased data set; (ii) the restricted ROI to ensure stability of the fit with theincreased data set; and (iii) the use of the 4FGL catalog to model Fermi -LAT point-like and extended sources.Every model component is characterized by an input spectrum and morphology, which are fitted to γ -ray datawith the adaptive template fitting algorithm implemented in skyFACT , and based on penalized maximum likelihoodregression. Hyperparameters in the regularization term of the likelihood control the allowed variation of spectral andspatial free parameters, preventing overfitting. Spectral and spatial modulation parameters (i.e. nuisance parameters)are allowed to vary to account for mis-modelling of the input templates. The minimization is performed by theL-BFGS-B (Limited memory BFGS with Bound constraints) algorithm. We refer to [25] for more details about thetechnical implementation. The diffuse γ -ray (spectral and spatial) model components we input are: (i) An isotropicspatial component with the best-fit IGRB spectrum from [41]; (ii) an inverse Compton spectral and spatial componentcomputed for a typical scenario of cosmic-ray sources and propagation parameters with the DRAGON code [42]; (iii) threerings for the spatial distribution of photons from π decay, as traced by the sum of atomic and molecular hydrogendistribution and available within the GALPROP public release (the π decay input spectrum is taken from [44]),(iv) the Fermi bubbles with spectrum from [43] and a uniform geometrical template as input morphology, and (v)the GCE. For the latter, we test different spatial models: NFW templates with slopes equal to 1 (NFW100) and 1.26(NFW126) and a bulge template, composed by a boxy-bulge and a nuclear bulge as modeled in [11]. Additionally,we refit all point-like and extended sources at the position of 4FGL cataloged sources. We refer to [11, 25] for detailsabout the γ -ray model. Spectral and spatial uncertainties on the model components are set by regularization terms.We allow variations as in run5 of [25] for all components, except for the additional GCE template. For the GCE, we https://galprop.stanford.edu/ dN/dS , see Sec. , we run fits for different values of the spatial smoothing hyperparameter, η , forthe gas and Fermi bubble templates. As defined in [25], the spatial smoothing hyperparameter is η = 1 /x , where x is the admitted variation between neighboring pixels. The reference values are η g = 25 for the gas and η b = 4 for the Fermi bubbles. By varying the smoothing scale, we therefore check the bias induced by mis-modeling at small scales.This updated analysis confirms previous findings about the preference for a bulge morphology of the GCE on topof DM-only templates. In Tab. SII, we report the log-likelihood values of skyFACT fits with different GCE templates.From the likelihood values, using the δ − χ statistics [11], one can compute the evidence for any additional templatein nested models. In this case, we find that: Adding an NFW template (no matter the slope) on top of a model withthe bulge already included does not significantly improve the fit (3 . σ for NFW100, and 4 . σ for NFW126), whileadding a bulge template on top of an NFW-only model does (13 . σ for NFW100, and 11 . σ for NFW126). In general,NFW100 provides larger residuals than NFW126.Finally, we also run skyFACT fits in the OG control region ( | b | < ◦ , 60 ◦ < | l | < ◦ ), where we do not findevidence for the additional bulge component, see Tab. SII. ROI skyFACT run − L IG r5 noGCE r5 NFW126 r5 NFW100 r5 bulge r5 bulge NFW126 r5 bulge NFW100 r5 noGCE r5 bulge Log-likelihood values for skyFACT fits with various GCE templates. Results for the 40 ◦ × ◦ IG ROI and for theOG ROI, for an unconstrained GCE spectrum.
THE
ANALYSIS
In this section, we provide further details on the analysis. We start by a small review of the methodand implementation as presented in [27, 28, 30], with some details specifically important for the analysis of the IG.The
Fermi -LAT dataset used for the analysis is then described, before focusing on the model parameters andpriors for the analysis.
The method
Different implementations of photon-count statistical methods applied on
Fermi -LAT data have been presented [27](also called ), and [36] (also called NPTF). The method have been extensively presented in [27, 28, 30],to which we refer for any further detail. We here briefly highlight the variations introduced in this work for theapplication to the IG with skyFACT templates.In the actual implementation of the , the probability p ( p ) k of finding k photons is pixel-dependent, p denotingthe evaluated map pixel. Different photon sources will contribute to the p ( p ) k with different statistics. Truly diffuse,isotropic emissions will contribute to p ( p ) k with counts following a Poissonian distribution. The presence of non-Poissonian sources, such as PS, and more complex diffuse structures alters the shape of p k , which permits to investigatethese components by means of the of the data (see [27] for details).The total diffuse contribution is given by the sum of the Galactic diffuse emission (as derived by skyFACT or takenfrom existing models, see next), a diffuse component describing the GCE as derived by skyFACT following a bulge orDM morphology, and a truly isotropic diffuse emission. Being x ( p )diff the number of diffuse photon counts expected ina map pixel p we have: (see also [30]) x ( p )diff = A gal x ( p )gal + A GCE x ( p )GCE + x ( p )iso F iso F iso . (S2)3 counts/px counts/px counts/px FIG. S3. Left panel: Region of interest for the analysis in the
Inner Galaxy . Middle panel:
Outer Galaxy . Right panel: extragalactic . The
Fermi -LAT counts per pixel in the energy bin [2, 5] GeV are reported. F iso is the integral flux of the isotropic diffuse emission, and the expression in Eq.(S2) permits to use directly F iso as asampling parameter, in order to have physical units of flux. The skyFACT templates for the Galactic diffuse emissionand GCE enter the fit with the best-fit normalization as found within the corresponding skyFACT analysis,and we allow A gal , GCE as additional normalization factors. When dealing with real γ -ray data, one has to take intoaccount source-smearing effects due to a finite detector point-spread function (PSF), which cause the photon fluxdetected from a given point source to be spread over a certain area of the sky (i.e. adjacent pixels, when skymaps aredivided in pixels). We correct for the PSF effect by statistical means as detailed in [27]. Fermi -LAT dataset for
We restrict the analysis of
Fermi -LAT data to events in the quartile with the best angular reconstruction, i.e.evtype=PSF3. The spatial binning of photon events is performed using the HEALPix equal-area pixelation scheme[37] with resoluton parameter κ = 7 [37] (the total number of pixels covering the entire sky is N pix = 12 N and N side = 2 κ ). The analysis is restricted to the 2 − analysis in smaller energy bins might allow to further investigate the nature of IPS, whichis left to future work.A benchmark upper flux cut of 10 − ph cm − s − is applied. We do not mask resolved sources, but we model themtogether with fainter PS to reduce possible systematics connected to the source masking or to additional Poissoniantemplates. A lower flux cut is only effectively applied within the hybrid method, see next section.To dissect the source-count distribution of PS in the IG we restrict to specific ROI. The fiducial IG ROI is definedby a 20 ◦ × ◦ square around the GC. As discussed in [19], the choice of the ROI can influence the results of photon-count statistical analysis, and needs to be chosen carefully. We found the 20 ◦ × ◦ square optimal to contain enoughstatistics to constrain the IPS parameters, and not too large, to avoid possible over-subtraction of the Galactic diffuseemission [19]. Nevertheless, we tested that our main conclusions are unchanged when going to 30 ◦ × ◦ . Althoughwe cannot test directly in the fit the preference for different PS spatial distributions, we can nevertheless verify ifthese PS are truly isotropic over the ten-degree scale of the IG ROI, and if they are mainly Galactic or extragalactic.To this end, the IG ROI is partitioned in radial and longitudinal slices to build the profiles shown in Fig. 2. For thesake of definiteness, the corresponding portions of the IG are illustrated in Fig. S4.Furthermore, two additional ROIs are used to investigate the nature of the IPS found in the IG ROI: The outerGalaxy (OG, | b | < ◦ , 60 ◦ < | l | < ◦ ) and the extragalactic ((EG, | b | > ◦ , | l | > ◦ ) ROIs. The first is meantto represent Galactic IPS, but away from the GC, while the second to provide the truly extragalactic population ofIPS. We define the OG such that, from the skyFACT runs, we do not have any longer evidence for a GCE emission.Since there is no evidence for the GCE, the diffuse models obtained with and without this additional component areconsistent. For the IG and OG ROIs we also cut the innermost region along the Galactic plane. We test different cuts | b | > . ◦ , ◦ , ◦ , finding consistent results. We stress that a latitude cut of 2 ◦ was used in past analyses [15, 17, 19],while we demonstrate the consistency of our results down to 0 . ◦ for the first time.An illustration of these three ROIs is provided in Fig. S3. For each ROI, the counts per pixel of Fermi -LAT datain the energy bin 2 − − − −
10 0 10 20 30
GLON [deg] − − − G L A T [ d e g ] − − −
10 0 10 20 30
GLON [deg] − − − G L A T [ d e g ] FIG. S4.
Radial and longitude slices.
The Inner Galaxy division in radial (longitude) slices in the left (right) panel used forcomputing the source number density depicted in Fig. 2, with consistent color coding.
Model parameters, priors, fitting procedure
We recall that the fits to
Fermi -LAT data are performed with the following components: Isotropic diffuseemission template, diffuse emission template (optimized or not through skyFACT fit), IPS population with source-countdistribution per unit flux dN/dS and, if included, a diffuse template describing the GCE, following a Galactic bulgeor DM morphology. The IPS population is described by a unique multiple broken power law (MBPL, see Eq. (1))with two or three free breaks, while the Poissonian components have one free normalization each. The normalizationconstant in Eq. (1) is fixed to S = 5 · − ph cm − s − . The results presented in the main text have been obtainedusing a Galactic diffusion emission template optimized by skyFACT , to reduce background systematics. Results usingalternative diffuse emission templates are discussed in the following. Both the diffuse emission templates and theGCE templates used as inputs for the fits are the optimized outputs of skyFACT , and enter the fits withan additional free normalization, A gal and A B / NFW126 , respectively. The pixel-dependent likelihood function L ( Θ ) isdefined following the L2 method in [27]. In this way, the spatial morphology of the skyFACT diffuse templates is takeninto account in the fits. In fact, while below the sensitivity of the method the PS in the ultra-faint regimeare in principle degenerate with a diffuse emission, our pixel dependent likelihood is sensitive to the morphology ofthe diffuse spatial templates. The full list of free parameters Θ, along with their prior intervals are summarized inTab. SIII.To sample the posterior distribution P ( Θ ) = L ( Θ ) π ( Θ ) / Z (where π ( Θ ) is the prior and Z is the Bayesian evidence)the MultiNest framework [38] was used in its standard configuration, setting 1000 live points with a tolerance criterionof 0.2. One-dimensional profile likelihood functions [39] for each parameter are built from the final posterior sample inorder to get prior-independent frequentist parameter estimates. If not differently stated, best-fit parameter estimatesrefer to the obtained maximum likelihood parameter values. Consistent results are found for the Bayesian parameterestimation. The nested sampling global log-evidence ln( Z ) is used to build Bayes factors for model comparison, seemain text.Different fitting techniques have been introduced within the framework [27]. We here use the MBPL approachas benchmark, where the parameters of the MBPL in Eq. (1) are sampled directly. To check for possible systematicsintroduced in the ultra-faint regime, we also perform our analysis using the Hybrid approach (see [27] for details). Thisis characterized by fixing a grid number of nodes (i.e. fixed flux positions) for the dN/dS fit, around the sensitivitythreshold of the analysis. The Hybrid approach was introduced to address possible underestimation or bias in thereconstructed source-count distribution fit and its uncertainty bands at the lower end of the faint-source regime, asdemonstrated using Montecarlo simulations in [27]. This has been shown to alleviate possible bias when measuringthe dN/dS in the ultra-faint regime, which can cause a loss of sensitivity of the method.5 TABLE SIII. Prior types and ranges for the analysis of the IG. The first two parameters are in common with all thesetups, while the third is present only in the , analyses. The second block refers to the parameters forthe IPS when using the MBPL fit approach, while the last one to the Hybrid approach, where the MBPL was extended with anode. The normalizations A S , A nd1 are given in units of s cm sr − . The break positions S bn , snd1 are in units of ph cm − s − .The F iso is given in units of cm s − sr − . All other parameters are dimensionless.Method Parameter Prior Range N b = 2 (=3) A gal log-flat [0.1,10] F iso log-flat [10 − , 10 − ] A B / NFW126 log-flat [10 − , 10 ]MBPL A S log-flat [10 , 5 · ] S b1 log-flat [5 · − , 10 − ] S b2 log-flat [10 − (5 · − ), 5 · − ] S b3 log-flat [. . . (10 − ), . . . (5 · − )] n flat [-1, 3] n flat [-1, 3.5] n flat [-2 (1), 2 (3)] n flat [. . . (-2),. . . (2)]Hybrid A S log-flat [10 , 5 · ] S b1 log-flat [10 − , 10 − ] S b2 log-flat [10 − (2 · − ), 10 − ] S b3 log-flat [. . . (10 − ), . . . (2 · − )] n flat [2.5 ,4.3] n flat [1.3,2.3] n flat [1.3,3] n flat [. . . (1.3),. . . (3)] A nd1 log-flat [10 , 10 ] S nd1 fixed 5 · − (3 · − ) n f fixed -10 DIFFUSE EMISSION TEMPLATE SYSTEMATICS
One of the main novelties of this work is the fact to employ consistently diffuse emission models optimized using skyFACT within the . We here apply the to the IG using other widely used diffuse emission templates:The official spatial and spectral template released by the
Fermi -LAT Collaboration for
Pass 8 data (Official P8)( gll iem v06.fits , see Ref. [40]), and the models labeled A (modA) and B (modB), optimized for the study of theIGRB in [41]. We note that within modA and modB the
Fermi bubbles are not modeled. The results for the dN/dS ofthe IG are illustrated in the left (right) panel of Fig. S5 when cutting the innermost 2 ◦ (4 ◦ ).By using standard diffuse models (modA, modB and Official P8), we reconstruct spurious sources at ∼ × − ph cm − s − , well above the sensitivity of the . Such a peak of the IPS dN/dS disappears instead if we usediffuse emission templates as optimized with skyFACT . Large scale residuals are indeed reduced when allowing thespatial diffuse templates to be remodulated in the fit. Even in the absence of an additional GCE template, the skyFACT fit remodulates the diffuse components such to partially absorbs GCE photons, therefore reducing residualsand improving the fit with respect to standard diffuse models. Also, in this ROI, all the diffuse models, except the skyFACT one, do not properly reproduce the 4FGL catalog bright sources. We notice that the spurious IPS peak ofthe dN/dS corresponds to a peak of flagged 4FGL sources, further corroborating the conclusion that it is indeed aspurious reconstruction effect. We stress that any comparison with 4FGL cataloged sources is purely illustrative, andserves for cross-checking our results at high fluxes.We therefore confirm previous findings [19] that large residuals due to mis-modelling of diffuse emission induce abias in the reconstruction of PS in the inner Galaxy.We also identify spatially critical regions within the IG where this mis-modeling effect is more pronounced, notablythe Northern hemisphere (both West and East quadrants). This might be connected to the North/South asymmetryfound within the NPTF analysis of the GCE discussed in [17]. As shown in Fig. S6 the spurious IPS peak of the dN/dS reconstructed with the using the Official P8 template is found to be strongly pronounced in the NorthIG ROI in the same flux region as found in Fig. S5, while it is not present in the analysis of the South IG ROI.When using the diffuse emission templates as obtained with skyFACT , we find instead a smoother dN/dS which iscompatible within 1 σ uncertainty with the 4FGL unflagged sources. Moreover, the results for the dN/dS using6 − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > ◦ SkyFACTOfficial P8 modAmodB 4FGL w/o flag4FGL − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > ◦ SkyFACTOfficial P8 modAmodB 4FGL w/o flag4FGL
FIG. S5.
Diffuse emission systematics.
Source count distribution in the IG obtained from the analysis cutting the inner2 ◦ (left panel) and 4 ◦ (right panel). The black line is obtained from the when using the model for the Galactic diffuseemission obtained from skyFACT (without any component modeling the GCE, sF-noGCE ). The colored lines are instead obtainedfrom the using the official Fermi -LAT model for
Pass 8 (cyan line), or modA and modB (orange and indaco lines). Theblack (gray) points represent the count distribution of 4FGL sources (without any analysis flag, intended as a cautionary indexfor the reality of a source or the magnitude of its systematic uncertainties [33]). − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] North Inner Galaxy, | b | > ◦ , Official P8 4FGL w/o flag4FGL − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] South Inner Galaxy, | b | > ◦ , Official P8 4FGL w/o flag4FGL FIG. S6.
North and South Inner Galaxy.
Left (right) panel: Source-count distribution of the North (South) region of the IGobtained from the analysis. Results are here reported using the Official P8 and the sF-B models for the diffuse emission.Points as in Fig. S5. the Official P8 and the skyFACT diffuse templates in the South IG ROI are compatible within the obtained 1 σ bands. dN/dS MODELING SYSTEMATICS
The stability of the dN/dS results in the IG from the combined - skyFACT analysis of Fermi -LAT data wastested against a number of systematics.The dN/dS in the IG and OG is well described by a MBPL with two free breaks. We verified that an additionalfree break is not preferred by data, and that the MBPL obtained with three free breaks (see prior intervals in TableSIII) is compatible, within the uncertainties, with the case of two free breaks. This is illustrated in the left panel ofFig. S7.Fluxes from bright sources with
S > − ph cm − s − are not considered in the analysis. We verified that7 − − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > ◦ , MBPL N b = 2 , MBPL N b = 3 , Hybrid N b = 2 , Hybrid N b = 3 − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Inner Galaxy, | b | > ◦ η g = 25 , η b = 4 η g = 25 , η b = 25 η g = 400 , η b = 4 η g = 100 , η b = 4 η g = 11 , η b = 4 η g = 4 , η b = 4 FIG. S7.
Systematic for dN/dS reconstruction in the IG.
Source count distribution of the IG obtained from the 1pPDF analysiscutting the inner 2 ◦ . Left panel: Effect of the number of free breaks N b (dashed lines) and of the Hybrid fit approach withdifferent number of breaks and varying the node position. The dotted line illustrates the position of S nd1 for the correspondingHybrid fit. Right panel: Effect of the smoothing scales η g , b used to obtain the skyFACT diffuse emission. Points as in Fig. S5. shifting the upper flux cut down to S > − ph cm − s − does not change our main results. The only differencewe observe is a change in the position of the first flux break, driven by few, bright sources, always compatible with4FGL points and uncertainties within 1 − σ C.L. As for the low-flux regime, in our benchmark setup we donot include any lower flux cut. To test for possible effects connected to the faint end of the source-count distribution,we repeated the main analysis using the hybrid approach introduced in [27]. We set a fixed node at S nd1 , with theindex of the power-law component below the last node, n f = −
10, thus effectively suppressing possible contributionsin the ultra-faint regime below the fixed node. Note that a fixed node S nb0 at the lower limit of the prior for the lastfree break is technically imposed, since the first free node S nb1 is continuously connected to the MBPL componentwith a power law at higher fluxes. We tested different values for the position of S nd1 in the faint source regime,together with a MBPL with two or three free breaks. Results are summarized in the left panel of Fig. S7. To theextent we have tested, a node in the faint source regime at 3 − · − ph cm − s − does not affect the reconstructed dN/dS of the IG, which is well compatible, within 1 σ uncertainty bands, with the benchmark results discussed inFig. 1. In particular, the dN/dS is well compatible in the flux interval 10 − − − ph cm − s − , where the radialand longitude profiles are computed.We also tested the effect on the results of changing the smoothing scale η of the Galactic diffuse componentsin the skyFACT fit. The benchmark values used for the results illustrated so far are η g = 25 and η b = 4 for the gasand the Fermi bubbles components, respectively. Variations for η b = 25 and different values for η g = 400 , , , η corresponds to a highersmoothing in the skyFACT template used for analysis. A different smoothing scale in the skyFACT diffusetemplate could affect the reconstructed dN/dS in the , as this could leave residuals at small angular scales thatcould be wrongly attributed to PS by the . In Fig. S7 the corresponding results for the dN/dS are reported,compared to the benchmark results also shown in Fig. 1. This test is performed using the sF-B diffuse emission withinthe fit. The best fit and 1 σ band of the reconstructed source-count distribution are consistent for all theexplored variations in η b , g . The same is true for the best-fit parameters of the IPS. As for the other parameters,variations of η b do not affect significantly the best fit of A B and the ln( Z ). We instead observe a slight decrease(increase) of ln( Z ) for decreasing (increasing) η g . This is expected, as for smaller η g in the skyFACT fit the Galacticdiffuse emission gas template can more easily account for small variations between adjacent pixels, decreasing theoverall residuals. OUTER GALAXY AND EXTRAGALACTIC dN/dS
RESULTS
We finally report extended results on the OG and EG ROI, which are illustrated in Fig. S8. We note that in theseROI there are no 4FGL flagged sources, and the dN/dS is well described by a single power law from S ∼ · − ph8 − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] Outer Galaxy, | b | > ◦ − − − S [ph cm − s − ] − − − − S d N / d S [ ph c m − s − d e g − ] ExtragalacticOfficial P8modA4FGL w/o flag4FGL
FIG. S8.
Outer Galaxy and extragalactic dN/dS . Source count distribution of the outer Galaxy (left panel) and of theextragalactic ROI (right panel) obtained from the analysis. Outer Galaxy: Results for the dN/dS are reported for twocases: (i) using the skyFACT
Galactic diffuse emission model without any component modeling the GCE (black line), and for thebenchmark sF-B model, where an additional component modeling the bulge is added. Extragalactic: Results are here reportedusing the Official P8 model and modA. Points as in Fig. S5. cm − s − down to S ∼ − ph cm − s − . The IPS in the OG and EG are well consistent with 4FGL source countsfor bright sources. In both cases, we obtain compatible results among different Galactic diffuse emission models, bothfor skyFACTskyFACT