Gamma-ray image reconstruction of the Andromeda galaxy
LLAPTH-007/21
Gamma-ray image reconstruction of the Andromeda galaxy
C´eline Armand
1, 2, 3, ∗ and Francesca Calore Univ. Grenoble Alpes, USMB, CNRS, LAPTh, F-74000 Annecy, France Univ. Grenoble Alpes, USMB, CNRS, LAPP, F-74000 Annecy, France Astronomy Department, University of Geneva, Chemin d’Ecogia 16, 1290 Versoix, Switzerland (Dated: February 15, 2021)We analyze about 12 years of
Fermi -LAT data in the direction of the Andromeda galaxy (M31).We robustly characterize its spectral and morphological properties against systematic uncertaintiesrelated to the modeling of the Galactic diffuse emission. We perform this work by adapting andexploiting the potential of the
SkyFACT adaptive template fitting algorithm. We reconstruct the γ -rayimage of M31 in a template-independent way, and we show that flat spatial models are preferred bydata, indicating an extension of the γ -ray emission of about 0 . − . ◦ for the bulge of M31. Thisstudy also suggests that a second component, extending to at least 1 ◦ , contributes to the observedtotal emission. We quantify systematic uncertainties related to mis-modeling of Galactic foregroundemission at the level of 2.9%. I. INTRODUCTION
Distant galaxies constitute promising targets to un-derstand astrophysical processes occurring inside suchsystems. Over the past decade, several star-forming galax-ies have been detected in γ rays by Fermi -LAT, amongwhich are the Large Magellanic Cloud, the Small Magel-lanic Cloud, Andromeda galaxy, M82, NGC 253, NGC2146, and Arp 220. The γ -ray emission of these distantgalaxies is produced by the interactions of the large-scalepopulation of cosmic rays with the interstellar mediumand by high-energy sources such as supernova remnantsand pulsars therein [1]. The study of external MilkyWay-like galaxies can provide an independent and outsideperspective of the γ -ray astrophysical processes whichare in place in our own Galaxy as well [2]. Therefore, itwould bring out additional and complementary knowledgeto our understanding of the Milky Way, especially if wecan discriminate individual contributions to the γ -rayemission of such extragalactic systems. Necessary (butnot sufficient) condition to this end is that an extendedemission signal is significantly detected in γ rays fromthe external galaxy. With this ultimate goal in mind, inthe present work, we focus on our nearest galaxy neigh-bor, the Andromeda spiral galaxy (M31), located at adistance of approximately 785 ±
25 kpc [10] at Galacticcoordinates ( l = 121 . ◦ , b = − . ◦ ) [11].Its close proximity allows us to optically resolve itsstellar disk and bulge as two separate components. Thisdistinction is not possible in our Galaxy as the bulge ∗ [email protected] Several measurements of the distance have been performed rang-ing from 765 ±
28 kpc [3] to 790 ±
45 kpc [4]. Most of themare found around 775-785 kpc [5–9]. In our study, we adopt adistance of 785 ±
25 kpc [10] to be consistent with previous workson the subject that refer to this value. is obscured by the bright emission of the disk [1]. M31spans 3 . ◦ × ◦ on the sky [2, 12]. This spiral galaxyhas a total mass of (0 . − . × M (cid:12) [13–20]. Thestellar component accounts for (10 − × M (cid:12) ofwhich 30% is in the bulge and 56% is in the disk [19].Furthermore, disk galaxies are typically surrounded bya large cosmic-ray halo extending up to a few hundredsof kpc [21], as well as by a circumgalactic medium mademostly by ionised hydrogen which can extend up to thevirial radius [22].As for γ -ray studies, M31 has the advantage of lyingat high Galactic latitudes, away from the Galactic planewhich makes it less polluted by the diffuse γ -ray fore-ground emission of the Milky Way [23]. M31 has beenthe object of several dedicated γ -ray analyses aimed atits spectral and morphological characterization. The firststudy in this direction was performed by the Fermi -LATcollaboration who analyzed about two years of
Fermi -LAT data, detecting M31 as a point-like source with a5 . σ significance and founding marginal evidence (1 . σ )for its spatial extension [24]. A seven-year data analysisshowed evidence for point-like source detection at ∼ σ ,and a 4 σ preference for extended emission, compatiblewith either a Gaussian distribution of width 0 . ◦ ± . . ◦ ± .
05. Several studiestested the extension of M31 in γ rays using various modelsfor its morphology [2, 12, 23, 25, 26], mostly confirmingthe preference for a centrally concentrated emission. Ad-ditionally, [25] found weak evidence for the presence of“Fermi bubble”-like structures perpendicularly to M31’sgalactic plane, while [26, 27] showed the existence of anextended excess up to 120 −
200 kpc (8 − ◦ ) away fromthe center of M31.Interestingly, no evidence for γ -ray emission correlatedwith M31’s gas-rich regions or star formation activitywas found. On the contrary, a mild correlation with in-frared stellar templates emerged, hinting to a possibleorigin from old stars [1]. Such a centrally concentrated a r X i v : . [ a s t r o - ph . H E ] F e b signal, perhaps associated with an old stellar populationin M31, was suggestive of another longstanding excessin γ -ray astrophysics, the so-called Fermi
Galactic cen-ter GeV excess [28–38]. A possible connection betweenthe two signals was explored by [39], which discussedthe possibility of a common origin from primordial anddynamically formed millisecond pulsars. Possibly, M31extended signal can also be compatible with dark matterannihilation from an adiabatically contracted dark matterdensity profile [2].Despite the growing evidence for M31 γ -ray spatialextension, the main limitation of current analyses remainsthe evaluation of the systematic uncertainties related tothe contamination of the Galactic diffuse emission, asdiscussed also in [2, 26].In what follows, we use the SkyFACT adaptive templatefitting algorithm [40] to characterize the γ -ray emissionin the direction of M31. Our goal is twofold: first, weaim to robustly detect M31 as an extended source againstforeground model systematic uncertainties, to character-ize its spectral and spatial distributions, and to test thepresence of multiple γ -ray emission components. To thisend, we will for the first time include in the fit M31 mor-phology templates tracing its stellar distribution from [19].We willl perform spectral and spatial fits of several mod-els to the data and make the first proper comparison ofour non-nested models. Secondly, we seek to reconstructM31’s morphology in a fully data-driven way based onthe image reconstruction feature of the SkyFACT code.This work represents the first study where a “template-independent” approach is applied to the analysis of anextended γ -ray signal. This image reconstruction featureallows us to build the galaxy intensity profile and deriveits parametrization.The paper is organized as follows: In Sec. II, we presentthe data selection and the fitting algorithm. Section IIIbriefly describes how we model the different emission tem-plates required in the analysis. Section IV illustrates theresults obtained in the case of standard template fittingtechniques and it is meant to reproduce and strengthenliterature results. In Sec. V, we exploit the innovativefeatures of SkyFACT and demonstrate that the evidencefor extension is robust against foreground contamination.In Sec. VI, we model-independently reconstruct the mor-phology of M31 and provide its intensity profile. Wefinally discuss the results and conclude on this work inSec. VII.
II. DATA SELECTION AND FITTINGALGORITHMA. Data selection
We use 624 weeks (12 years) of
Fermi -LAT data fromAugust 4th, 2008 to July 16th, 2020, collected from a10 ◦ × ◦ region of interest (ROI) centered at ( l = 121 ◦ , b = − ◦ ). We select Pass 8 SOURCE class events with Point Spread Function (PSF) PSF0+PSF1+PSF2+PSF3type and use the corresponding instrument response func-tions (IRFs) P8R3 SOURCE V2. We apply a cut onmaximum apparent zenith angle ( z max = 105 ◦ ), andrecommended data-quality filters (DATA QUAL >
0) &(LAT CONFIG==1). We analyze data in the energyrange 300 MeV to 100 GeV, logarithmically spaced into22 bins. Table I summarizes the data selection criteria.This work makes use of the
Fermi
ScienceTools v11r5p3software package . We spatially split the data of the ROIinto 200 ×
200 angular bins of size 0 . ◦ in cartesian skyprojection. Event selection criteriaEvent class Pass 8 SOURCEEvent type PSF0+PSF1+PSF2+PSF3IRFs P8R3 SOURCE V2 z max ◦ ROI size 10 ◦ × ◦ ROI center l = 121 ◦ , b = − ◦ Pixel resolution 0.05 ◦ Pixel binning 200 ×
200 pixelsSky projection Cartesian “CAR”Energy range 300 MeV to 100 GeVEnergy binning 22 binsFilters (DATA QUAL >
0) & (LAT CONFIG==1)
TABLE I. Summary of the selection criteria applied to thedata using the
Fermi
ScienceTools . B. SkyFACT
To perform the γ -ray analysis of M31 we use SkyFACT [40], a code which combines a hybrid approachtemplate fitting and image reconstruction for studyingand decomposing the γ -ray emission. SkyFACT offers theadvantage of accounting for expected spatial and spectraluncertainties in each model emission component of thefit, contrary to standard template fitting methods wherethe spatial distributions of model components are fixed.The algorithm is based on penalized maximum likelihoodregression, with additional nuisance parameters that ac-count for uncertainties modeling the imperfections of themodel templates. More details on the analysis code canbe found in Appendix A.The detection of the source and the evidence for itsextension are probed through a likelihood ratio (LR) sta-tistical test, denoted TS, between our various modelsincluding M31 and our model without the source. The TSresults are then interpreted in terms of standard devia-tions σ using the Chernoff theorem, describing a mixturebetween a χ and a Dirac function. A detailed descrip-tion of the statistical analysis framework can be found inAppendix B. http://fermi.gsfc.nasa.gov/ssc/data/analysis/ III. MODELING THE γ -RAY EMISSION In our model the γ -ray emission in the ROI is made upby the following contributions: the Galactic interstellaremission (IEM), the Isotropic Gamma Ray Background(IGRB), 13 Point-like Sources (PS) as reported in the4FGL [11], and M31. No extended source is presentin the chosen ROI. We describe below how we modelthe individual templates (spectral and spatial), used by SkyFACT in the fit of the data.
A. IEM, IGRB and PS modeling
The Galactic IEM represents the foreground emissionfrom the Milky Way, and includes the contribution fromInverse Compton, π decay, and Bremstrahlung γ -rayproduction mechanisms. We adapt the Galactic IEM gll_iem_v06.fits provided by [41] to our energy andspatial binnings, and normalize the spatial template by thesum of all pixels for each energy bin. The normalizationof each energy bin corresponds to the spectral profile ofthis diffuse component.As spectral template of the IGRB component,we use P8R3_SOURCE_V2.txt [41] associated to gll_iem_v06.fits . We produce a flat templateto model the spatial part, as the IGRB is assumed to beisotropic after exposure correction of the LAT [42].PS best-fit spectra are extracted from the 4FGL catalog, gll psc v22.fits , of
Fermi -LAT sources [11], and areused as input spectra for the
SkyFACT fit.
B. M31 modeling
As for M31 galaxy, we model its γ -ray emission usingvarious templates motivated either by previous works [2,12, 23] or by studies of the stellar distribution of thegalaxy [19], with the final goal of understanding whichM31 morphology is preferred by data.As for the morphology ( i.e. spatial component) of M31,we test the following distributions: a. Gaussian spatial profile. We build a first set ofspatial templates for M31 using a Gaussian function whosewidth σ defines the 68% containment extension angle. Weproduce Gaussian templates of width ranging from 0 . ◦ (PS limit) to 1 ◦ around the position of M31. b. Uniform spatial profile. A second set of spatialtemplates is constructed with uniform disks of radii ex-tending from 0 . ◦ (PS limit) to 1 ◦ . We voluntarily use an older version of the IEM template model( v06 ) since the latest one ( v07 ) may include some emission in thedirection of M31, and due to M31 gas not subtracted from thatIEM model. c. Einasto spatial profile.
Finally, we model M31 asthree distinct stellar components for the nucleus, bulge,and disk, using three ellipses inclined at 77 . ◦ [43]. Thedensity distribution of each of them follows the Einastolaw whose parameters are defined in [19] and reads: ρ Ein ( r ) = ρ c exp (cid:40) − d N (cid:34)(cid:18) rr c (cid:19) /N − (cid:35)(cid:41) , (1)where d N is a function of the coefficient N defined as d N ≈ N − / . /N for N ≥ .
5. The parameter ρ c corresponds to the density at a distance r c whichencloses half of the total mass of the component. Thevariable r is the distance of a given pixel from the centerof M31 expressed in terms of line of sight s as: r ( s, b, l ) = ( d + s − ds cos b cos l ) / (2)where ( b, l ) are the relative galactic coordinates to thecenter of M31 and d is the distance between the observerand M31.We integrate the Einasto profile over the line of sighttowards each pixel p of the ROI in order to obtain aspatial template. The integration is performed accordingto: T p ( b, l ) = (cid:90) s ρ Ein ( r ( s, b, l )) ds , (3)whose upper and lower limits are derived by solving Eq. 2for s with r = R vir where R vir is the virial radius, set to213 kpc [19]. The two solutions obtained for s are: s max / min = d cos b cos l ± (cid:113) d [(cos b cos l ) −
1] + R (4)We consider three sub-components to model M31, i.e. the nucleus, the bulge, and the disk, each defined by anellipse of ratio q and an Einasto parametrization. Wenote that 68% of the total emission expected for eachcomponent is contained in an ellipse of major (minor)axis 2.15 ◦ (0.37 ◦ ) in the case of the disk, and 0.26 ◦ (0.18 ◦ )in the case of the bulge. For the nucleus template, in-stead, about 95% of the total emission comes from thecentral pixel. Table II presents the parameters ρ c , r c , N , d N , and q characterizing the spatial template of eachsub-component. In what follows, we will consider thesetemplates individually or combined in order to model theemission of M31.As for the M31 spectral template, we use the M31spectral energy distribution as provided in the 4FGLcatalog, and modeled by a power law with exponentialcutoff. IV. STANDARD TEMPLATE FITTING
Although
SkyFACT was developed for simultaneous spa-tial and spectral template fitting, it can also be used to
Ellipse parameters and Einasto templates
Component ρ c r c N d N q Nucleus 1.713 0.0234 4.0 11.67 0.99Bulge 0 .
920 1.155 2.7 7.77 0.72Disk 0 .
013 10.67 1.2 3.27 0.17
TABLE II. Parameters of the Einasto templates where ρ c [ M (cid:12) . pc − ] is the density, r c [kpc] is the radius enclosing halfof the total mass of the sub-component, d N is a function of thecomponent N , and q the ratio between the major and minoraxes of the ellipse. These values are taken from [19]. perform traditional or standard template fits (STF). Theanalysis within this setup allows us to directly comparewith the previous literature, as well as with the other SkyFACT runs we will explore in the following sections.For this purpose, we set all spatial hyper-parameters sothat the spatial templates are kept fixed to the input mod-els in the fit, for all components. By contrast, the spectralparameters are left completely free to vary, i.e. bin-by-binenergy independent fit.
A. Detection of M31 as PS
We first perform STF to probe the evidence for M31 in
Fermi -LAT data in the energy range between 300 MeVand 100 GeV. We compare the fit of a model that onlycontains IEM, IGRB, and the 13 PS, i.e. our
STD baselinemodel , to the fit of a second model including an additionalcomponent modeling M31. For our PS-limit Gaussianspatial template (0 . ◦ width) and input spectrum takenfrom the 4FGL, we compute the significance of an ad-ditional source in a nested model with bin-by-bin freespectral parameters according to the Chernoff theorem(see Appendix B). M31 is significantly detected in this γ -ray energy range with a 7 . σ significance. Figures 1and 2 show respectively the flux of each component of themodel fitted to the data, and the residuals in the absenceof M31 compared to the fit including M31. We found areconstructed spectral flux for M31 consistent with theone reported in the 4FGL catalog and parametrized bya power law with exponential cut-off. We note that thelack of statistics generates very large uncertainties on thereconstruction of the M31 flux in the highest energy bins.From the residuals map, we can see that the addition of asource at the M31 position captures some of the residualsleft out by the model without M31. B. Evidence for the extension of M31
In order to (i) demonstrate the preference for an ex-tended emission morphology and (ii) measure the ex-tension of the signal, we perform a scan over the 68% extension angle θ of the Gaussian spatial profile . Wetherefore obtain the log-likelihood ( − L ) profile as afunction of θ , as shown in Fig. 3. The log-likelihood pro-file is minimized by non-zero extension angles, indicatinga preference for an extended γ -ray emission for the M31model component. The highest detection significance ofthe source is σ = 9 .
2, and it is associated to the Gaussiantemplate of width θ = 0 . ◦ .In order to estimate the best-fit extension angle, we fitthe log-likelihood profile with a cubic spline (green linein Fig. 3). The minimum of the curve corresponds to θ = 0 . ◦ ± .
07. The uncertainties on the extension arecomputed by determining the two values of θ at which: − L ( θ ± ∆ θ ) = − L (ˆ θ ) + ( −
2∆ ln L ( θ )) (5)where 2∆ ln L ( θ ) = 1 .
00 for a coverage probability of68 .
27% in the case of one parameter [44]. The 1 σ statisti-cal uncertainties of each side of the distribution are thenderived by the following subtractions: σ + θ = θ + − ˆ θ ; σ − θ = θ − − ˆ θ. (6)This result is compatible with the value of [2] who founda size of 0 . ◦ ± .
10 using a Gaussian template.
C. Systematics due to PSF mis-modeling
We quantify the impact of systematic uncertainties onthe angular extension θ due to the possible mis-modelingof the PSF. To this end, we build two alternative PSFkernels based on the prescription of [45], where the averagePSF, P ( θ, E ), is modified according to: P min ( θ ; E ) = P ( θ × [1 + f ( E )]; E )(1 + f ( E )) P max ( θ ; E ) = P ( θ × [1 + f ( E )] − ; E )(1 + f ( E )) − , (7)where θ is the reconstruction angle and E is the energy.The function f ( E ) corresponds to the scaling functionfor the relative PSF uncertainty depending on the energyand is defined by the following: f ( E ) = (cid:26) . E ≤
10 GeV0 .
05 + 0 . ( E/
10 GeV)
E >
10 GeV , (8)where the uncertainty is constant at 5% below 10 GeV andincreases up to 25% at an energy of 1 TeV. Using these twobracketing models of the PSF, we scan over the angularextension θ and we find the corresponding best-fit angularextensions: θ P min = 0 . ◦ +0 . − . and θ P max = 0 . ◦ +0 . − . . We will discuss below the extension measurement for other typesof template. E [GeV]10 E d N / d E [ G e V / c m / s / s r ] Fit without M31 (STF)
TotalIEM IGRBPoint Sources Data 10 E [GeV]10 E d N / d E [ G e V / c m / s / s r ] Fit with M31 - Gaussian template of width 0.001° (STF)
TotalIEM IGRBM31 Point SourcesData
FIG. 1. Spectra of each component for the fit without M31 (left) and with M31 modeled as a 2D Gaussian of width 0 . ◦ (PS-limit) (right) . b [ ° ] R e s i d u a l c o un t s b [ ° ] R e s i d u a l c o un t s b [ ° ] R e s i d u a l c o un t s b [ ° ] R e s i d u a l c o un t s FIG. 2. Residual maps for the fits without (left column) and with (right column) the M31 component. The residual mapsshow the difference between the data and the theoretical model. The top row shows the residuals of the whole ROI wherethe center of M31 is marked by a green cross, while the bottom row is a zoom-in of the area within the green square. Thezoom-in residual maps show that the model including the additional M31 component (bottom right) captures some of the γ -ray residuals clearly present when M31 is not included in the model (circled in green) (bottom left) . Extension angle [°] L i k e li h oo d l n L STF Gaussian spatial template
InterpolationSTF results
FIG. 3. Log-likelihood, − L , profile from the scan over theextension angle θ (dots) of a Gaussian spatial template. Thesolid green line corresponds to the interpolation from whichwe estimate the best-fit extension angle and its statisticaluncertainties. To estimate the systematic uncertainty on the quantityof interest, θ , we follow Eq. 4 of [45], and compute thedispersion between the nominal value of θ and the valuesobtained with the two bracketing PSF kernels P min ( θ ; E )and P max ( θ ; E ). We find the systematic uncertainty on θ due to PSF mis-modeling to be 1.5%. V. SEMI-ADAPTATIVE TEMPLATE FITTING
After having obtained results for the STF setup, weuse
SkyFACT to test systematic uncertainties related tothe uncertainties on the foreground emission template.In particular, by allowing some more freedom on thespatial hyper-parameters of the IEM diffuse component(in addition to the full freedom on the spectral templatealready tested by the STF approach), we can test (i)whether or not the extension of the source is still preferredagainst uncertainties on the IEM template, and (ii) therobustness of the reconstructed spectra of M31 and IEMcomponents against the freedom given to the baselineparameters.In what follows, we refer to this setup as Semi-Adaptative Template Fitting (sATF). We define threeconfigurations of baseline models allowing a 30%, 50%,and 100% variation of the IEM spatial parameters. Wetest the evidence for M31 and measure its extension ontop of each of these baseline models and for differentM31 spatial profiles: Gaussian, uniform, and Einastotemplates.
A. sATF with Gaussian spatial templates
We perform three scans over θ , one for each new base-line model, and we derive the best-fit angular extension from the cubic spline interpolation of the log-likelihoodprofile. The three IEM configurations give compatiblevalues with best-fit extension angles of θ = 0 . ◦ ± . θ = 0 . ◦ +0 . − . (IEM 50%), and θ = 0 . ◦ +0 . − . (IEM 100%). These results are also in agreement withthe STF results of Sec. IV. We evaluate the systematicuncertainty on θ due to the IEM mis-modeling to be2.9%, by computing the dispersion among the nominalvalue from the STF scan (no freedom allowed on the IEMcomponent) and the best-fit extensions obtained with thethree sATF configurations, with Eq. 4 of [45].For all three IEM configurations, the reconstructedspectra of the IEM and M31 model components result tobe compatible with each other. In Fig. 4, we compare theIEM (top panel) and M31 (bottom panel) best-fit spectraobtained with the three sATF configurations and withthe STF where no freedom on IEM spatial parametersis allowed. All spectra are extracted from the fits whereM31 is described by a Gaussian template of 0 . ◦ width.The reconstruction of the spectra for both the IEM andM31 model components is therefore robust against thevariation allowed on the IEM spatial parameters, andcompatible with STF results.In Fig. 5, we show the remodulation map of the IEMcomponent for the three IEM configurations, for a Gaus-sian template of 0 . ◦ width. The remodulation mapindicates how much a component contributes to the totalflux observed in each pixel, and what is its variation w.r.t.the input template. We notice that the IEM remodulationmap for the 100% configuration shows some significantre-modulation around the position of M31. This over-absorption is instead absent if we decrease the freedomon the IEM spatial parameters to 30% and 50%. Wefind a negligible relative difference when comparing theoutputs of the 50% and 30% configurations (0.02% onaverage over the ROI, with 1% variation for the pixelcorresponding to M31 center). Instead, the 100% configu-ration remodulation map show variations at the % levelwith respect to the 50% and 30% configurations (about3% on average over the ROI, with 5-6% variation for thepixel corresponding to M31 center). We note that STFand sATF fits provide a reconstruction of the IEM compo-nent with a statistical uncertainty of 4-5% at low energies, i.e. the most constrained energies, that increases up to28% at the highest energies where there is less statistics.Given that the best-fit extension angles and the recon-structed spectra are well compatible for the three IEMconfigurations, we adopt the IEM 30% configuration forthe rest of this study. B. sATF with uniform spatial templates
We perform a scan where M31 is modeled by a flatdisk of radius extending from 0 . ◦ (PS-limit) to 1 ◦ .From the log-likelihood profile interpolation, we find twominima in this case: the first minimum gives an extensionof θ = 0 . ◦ ± .
04 which is compatible with the value E [GeV]10 E d N / d E [ G e V / c m / s / s r ] Diffuse Emission Flux
STF - IEM 0%sATF - IEM 30%sATF - IEM 50%sATF - IEM 100%10 E [GeV]10 E d N / d E [ G e V / c m / s / s r ] M31 Flux
STF - IEM 0%sATF - IEM 30%sATF - IEM 50%sATF - IEM 100%
FIG. 4. Best-fit spectra of the IEM (top) and M31 (bottom) model components using four different baseline configurations,where we vary the freedom on the IEM spatial parameters:0% (STF), 30%, 50%, and 100%. For more clarity, each pointis slightly shifted on purpose on the energy scale with respectto the green points, placed at the true energies. of 0 . ◦ ± .
04 found by [2] for a uniform disk template.The second minimum is located at θ = 0 . ◦ +0 . − . . Bothextensions have been detected with significance of 8.5 σ .The first minimum is compatible with the bulge size.The second one can be instead explained by the presenceof a more extended structure such as the stellar disk [19],or a bubble-like component [25], while it does not seem tobe consistent with the extent of the M31 disk as observedin optical light (3 . ◦ ) [12].Motivated by previous evidence for bubble-like struc-tures lying above and below M31 galactic plane [25], wetest a specific model for these structures. We performthe fits of nested models where we add an M31 bubblecomponent to our previous best-fit models, the 0 . ◦ widthGaussian bulge and the 0 . ◦ radius uniform disk. Wemodel the bubbles spectrum with a power law of index − . NASA/IPAC EXTRAGALACTIC DATABASE - https://http://ned.ipac.caltech.edu/ . disks of 0 . ◦ radius lying perpendicularly to the M31galactic plane. The center of each bubble is offset by0 . ◦ from the center of M31. We find only marginalevidence (2 . σ at best) for the bubble component, a bitless significant than what found by [25] who reported a5 . σ evidence. C. sATF with Einasto templates
We finally test Einasto spatial templates as motivatedby recent stellar mass models of M31. In this case, we donot perform any scan over the Einasto template param-eters since we want to test specific models tracing M31stellar components: nucleus, bulge, and disk. We keepthe same spectral template of M31 extracted from the4FGL for the nucleus and the bulge, while we use a powerlaw of index -2.4 to model the spectrum of the disk [25].We test the evidence for each of these three components,separately, on top of our baseline model with sATF (IEM30% configuration). As for the results in Tab. III, thenucleus, bulge, and disk are all, individually, detectedat about 7 . − . σ . In Fig. 6, we display the best-fitspectra of these three templates when fitted as a singleadditional component. The reconstructed spectral fluxesof the nucleus and the bulge components follow the sametrend with a bump between 1-2 GeV and a cut-off athigher energies. In the case of the disk model, a similarbump feature is also observed at low energies. In addition,the spectral reconstruction shows a high energy tail inthe spectral flux. The authors of [39] report a similarobservation which could potentially be confirmed by thenext-generation of γ -ray experiments such as CTA [46]and LHAASO [47]. sATF results for Einasto spatial templates Model − L σ det no M31 311107.65 -N 311015.52 7.5B 311010.85 7.7D 310997.55 8.4N+B 311010.06 7.8N+D 310986.28 9.0B+D 310986.16 9.0N+B+D 310985.30 9.0TABLE III. Log-likelihood value, − L , and detection sig-nificance, σ det , of the M31 model over the model without M31.We consider Einasto spatial templates for nucleus (N), bulge(B), and disk (D) and their combinations. We also run fits for all possible pair combinations andfor the setup with all three templates in the fit. We com-pute the test statistic and significance of all combinationsof nested models to determine the significance of eachcomponent (Tab. IV): • The addition of a disk to the nucleus (bulge) modelimproves the fit with a significance for the disk at b [ ° ] R e m o d u l a t i o n I E M % b [ ° ] R e m o d u l a t i o n I E M % b [ ° ] R e m o d u l a t i o n I E M % FIG. 5. Remodulation map of the IEM component in the case of sATF for the Gaussian template (0 . ◦ width) for the IEMconfigurations 30% (top left) , 50% (top right) , and 100% (bottom) . The black cross indicates the center of M31 given bythe 4FGL. E [GeV]10 E d N / d E [ G e V / c m / s / s r ] Reconstructed fluxes of single Einasto components
NucleusBulgeDisk
FIG. 6. Comparison of M31 reconstructed flux for the individ-ual spatial Einasto models with either a nucleus, a bulge, or adisk modeling M31 emission. . σ (2 . σ ). Adding instead a nucleus or a bulgecomponent to the disk model shows a very smallimprovement of the fit with a significance of lessthan 1 σ . • The addition of a bulge to the nucleus model doesnot improve significantly the fit (0 . σ ). The addi-tion of a nucleus to the bulge does not show anyimprovement either. • The case of the fit with all three sub-components(nucleus, bulge, and disk) does not provide anysignificant improvement over the nucleus+disk andbulge+disk models.The preference for a disk-like extended structure isin agreement with our previous findings for the uniformspatial disk template. The extension of the Einasto disk(1.08 ◦ semi-major axis) is indeed compatible with thesecond minimum found in the uniform disk scan (0 . ◦ ). Einasto spatial templates: nested models comparison
Model σ Model σ N+B vs N 0.2 N+D vs D 0.8N+D vs N 2.9 B+D vs D 0.8N+B+D vs N 3.0 N+B+D vs D 0.9N+B vs B ∼ ∼ ∼ σ , of the nested models consideringcombinations of nucleus (N), bulge (B), and disk (D) templates. D. Akaike Information Criterion (AIC) modelcomparison
In order to properly determine which model best de-scribes the spatial distribution of M31, we compare ourdifferent non-nested templates using the AIC score, asdefined in Appendix C.In Tab. V, we quote the AIC differences between thevarious template combinations. Based on the ∆AIC val-ues, we can rank for the very first time the differentmodels used to describe the morphology of M31: thelower the AIC score, the better the model is. In general,models for which the ∆AIC relative to the model witha lower AIC is less than 2 are typically considered tohave substantial support from data [48]. The uniformtemplates with θ = 0 . ◦ ± .
04 and θ = 0 . ◦ +0 . − . turnout to be the best models we tested with an almost iden-tical score (AIC min ), followed by the nucleus+disk, thebulge+disk and disk-only Einasto templates, and finallyby the Gaussian template of width θ = 0 . ± . ◦ .We note that, from the log-likelihood results only ( i.e. nested model comparison), the addition of a bulge or anucleus to the disk model improves the likelihood butnot enough (0 . σ ) to claim the presence of an additionalcomponent. The AIC score, on the other hand, factorsin the number of effective parameters of the fits. In thecase of the bulge+disk model, the improvement in the log-likelihood is not enough to counterbalance the increaseof effective parameters due to the addition of the bulgecomponent. The nucleus+disk model, on the other hand,has about the same number of effective parameters – dueto the small size of the nucleus –, while the log-likelihoodimproves. Thus, for the case nucleus+disk, we get asignificantly better model. VI. IMAGE RECONSTRUCTION OF M31MORPHOLOGY
We finally take full advantage of
SkyFACT ’s flexibilityby performing an adaptive template fit allowing a 100%freedom on the spatial and spectral parameters of theM31 component. We keep the same input M31 spectralshape from the 4FGL as above. As for the spatial part,we do not assume any spatial input morphology, but werather define a large circular region of 1 . ◦ radius aroundthe source position, within which the M31 component isadjusted to the data in each pixel and each energy bin.The reconstructed morphology of M31 is shown in Fig. 7.From the reconstructed M31 image, we build the in-tensity profile of M31 shown in Fig. 8 by summing thecontributions of all pixels in concentric annuli spaced of0 . ◦ from 0 ◦ to 1 . ◦ . The intensity of each annulus isthen averaged over the solid angle ∆Ω annulus = N ∆Ω pixel ,where N is the number of pixels in a given annulus. Theannulii are centered around M31 4FGL position.The intensity profile appears rather flat and does notfollow a simple known function. We fit the intensity using0 ∆ AIC model comparison
Model Gaussian 0 . ◦ Uniform disk 0 . ◦ Uniform disk 0 . ◦ Einasto D Einasto B+DUniform 0 . ◦ -56.51 - - - -Uniform 0 . ◦ -55.93 0.58 - - -Einasto D -44.79 11.71 11.13 - -Einasto B+D -44.67 11.26 21.65 0.12 -Einasto N+D -55.33 1.18 0.60 -10.53 -10.66TABLE V. ∆AIC computed for (Row - Column). A positive (negative) ∆AIC indicates that the model listed in the correspondingcolumn gives a better (worse) fit to data than the model in the corresponding row. l o g I n t e n s i t y [ / c m / s / s r ] FIG. 7. Reconstruction of the γ -ray emission of M31 usingthe full adaptive template fitting (ATF) technique. The blackmarker indicates the position of the source given by the 4FGLat ( l = 121 . ◦ , b = − . ◦ ) [11]. F l u x w i t h i n e a c h a nnu l u s [ / c m / s / s r ]
1e 6
FIG. 8. Intensity profile of M31 built using concentric annuliiof width 0 . ◦ (dots). The green line shows the total fit of theintensity profile which corresponds to the combination of aGaussian function and a power law. a combination of a Gaussian function and a power lawgiven by: Ξ( x ) = a exp (cid:26) − ( x − µ ) σ (cid:27) + bx c , (9) where the best-fit coefficients are found to be a = 1 . × − , µ = − . × − , σ = 4 . b = 4 . × − , and c = 2 . R regressionsum of squares method, defined in Appendix D. We finda value of R = 0 . x ).These two functions in Ξ( x ) can be interpreted as asso-ciated to the two distinct components of M31 emissionalso hinted by our previous studies: the Gaussian functionwould model M31 central emission as it is peaked at thecenter of the source and decreases at larger radii, while thepower law could correspond to the extended structuresaround the galaxy, such as galactic bubbles, the outerpart of the disk, and/or a larger cosmic-ray halo. We findan inflection point of the intensity profile at 0 . ◦ by takingits second derivative. This corresponds to the angle atwhich the power law starts kicking in and may delimitthe apparent size of the bulge in γ rays. This value iscompatible with the result of [2] who found an extensionangle of 0 . ◦ ± .
04 using a uniform disk template tomodel M31’s bulge, and also marginally consistent withour first minimum of the uniform disk scan 0 . ◦ ± . VII. DISCUSSION AND CONCLUSIONS
In the present work, we analyzed the
Fermi -LAT γ -rayemission at GeV energies from the direction of M31. Were-assessed the evidence for an additional extended γ -raysource at the position of M31, and carefully characterizedits morphology by quantifying systematic uncertaintiesfrom the IEM foreground emission. To this end, weintroduced the following methodological novelties: • For the first time, we adapted and applied the
SkyFACT code to the analysis of extended γ -raysources, generalizing its scope (which was previouslyonly limited to the inner Galaxy analysis). • We fully exploited the power of
SkyFACT of account-ing for a large number of free parameters by definingthree setups which allowed us (a) to recover resultsfrom the literature and estimate the systematicuncertainties due to the PSF mis-modeling (STF,1Sec. IV), (b) to properly compare non-nested mod-els with various spatial templates, (c) to quantifysystematic uncertainties from IEM mis-modeling(sATF, Sec. V), and (d) to reconstruct the image ofM31 in a template-independent way (ATF, Sec. VI). • We tested a large number of M31 spatial templatesof three kinds, the Gaussian and uniform modelsmotivated by previous literature, and the Einastospatial templates which are directly motivated bythe different M31 stellar components (nucleus, bulge,and disk). • For the first time, we ran a proper model comparisonof our different spatial models for the emission ofM31, and ranked them based on their AIC score.This allowed us to be able to state what morphologyis preferred by the data, among the ones we tested.We summarize below our main findings: • We confirm the evidence for an extended γ -ray emis-sion from the direction of M31 ( ∼ σ detectionsignificance). • For a Gaussian spatial template, the best-fit an-gular extension is θ = 0 . ◦ ± .
07 (68.27% sta-tistical uncertainty). This determination is veryrobust against systematic uncertainties related tomis-modeling of the PSF and of the IEM. We quan-tify PSF systematic uncertainties at the level of1.5%. To estimate IEM-related uncertainties, weset up sATF runs with increasing freedom in thevariation of the IEM spatial parameters (30%, 50%and 100% allowed variation). We quantify the IEMsystematic error to be 2.9%. • When modeling M31’s morphology with a uniformspatial template, we interestingly find two minimaof the log-likelihood profile: the first one at an angu-lar extension of 0 . ◦ ± .
04 associated to M31 bulge,and the second at 0 . +0 . − . ◦ indicating the presenceof a more extended structure around M31, suchas the galactic bubbles [25]. We explicitly testedthe presence of bubbles-like structures lying aboveand below the M31 galactic plane on top of a cen-trally concentrated emission. We only find marginalevidence ( < σ ) for this additional emission compo-nent. • The use of spatial templates tracing the stellar dis-tribution of M31 nucleus, bulge, and disk compo-nents provides some more insight on the physicalnature of the γ -ray emission. We observe that thedata prefer to complement centrally concentratedcomponents (nucleus- or bulge-only models) withan extended disk-like component ( ∼ σ evidence),which has an angular extension of about 1 ◦ (semi-major axis). This is a further indication of thefact that the M31 γ -ray emission is indeed moreextended than what was previously found in most of the literature showing an evidence for the bulgeonly component [2, 12, 24]. On the other hand,we cannot claim evidence for a centrally concen-trated component (nucleus or bulge) on top of adisk-only model, as we found a marginal improve-ment of 0 . σ only. The case of a fit with all threesub-components (nucleus, bulge, and disk) insteaddoes not provide any significant improvement eitherover the nucleus+disk or bulge+disk models. • We properly compare spatial template modelsthrough the evaluation of the AIC score of eachmodel. Based on ∆AIC values, we find the datarather prefer flat templates to describe the γ -rayemission of M31. As a matter of fact, the uniformdisk templates come first in our ranking, followed bythe Einasto nucleus+disk model, then the Einastodisk and bulge+disk models, and finally the Gaus-sian templates. Again, this highlights the preferencefor a more extended emission, which goes beyondM31’s bulge. • We find that the M31’s reconstructed spectrum isextremely robust against all variations of fit setupsand spatial templates we tested. The emission fromM31 is significant up to, at most, 10 GeV in energy. • We reconstruct the morphology of M31 in a fullytemplate-independent way. To this end, we do notimpose any a priori spatial distribution for M31,and we leave the
SkyFACT spatial parameters tofreely adjust to the data. From the reconstructedangular intensity profile, we show that the M31intensity profile can be fitted by a two-componentparametrization made up by a Gaussian function,dominant within ∼ . ◦ , and a power law at largerradii. The latter component could correspond to asuperposition of several contributions such as theouter part of the disk, galactic bubbles similar tothe ones of the Milky Way, and/or a large comic-rayhalo surrounding M31 [27].Our results therefore indicate that the emission fromM31 extends beyond its nucleus and bulge, with significantemission up to at least 1 ◦ . This can be compatible eitherwith a flat disk-only emission, or with a two-componentmodel such as a nucleus+disk. We stress that thesecomponents should trace the stellar distribution in M31,rather than gas regions given the lack of evidence for γ -ray gas-correlated emission [49]. A stellar origin ofM31’s extended emission is also strongly supported bythe reconstructed spectral energy distribution, which isnot compatible with the Milky Way IEM spectrum andclearly distinguishable from it. A plausible explanationfor the emission can be, for example, a large populationof millisecond pulsars in the disk and bulge of M31 [39].This finding is of relevance especially for models tryingto interpret, in a common framework, the γ -ray emissionfrom M31 and from the Milky Way inner region, a.k.a. theGeV excess [12, 39].2With its better resolution and wide energy coverage,the upcoming γ -ray Cherenkov telescopes CTA [46] andLHAASO [47] will contribute to demystifying the natureof M31’s emission, by, for instance, testing the high-energytail of its spectrum and possibly detecting an interstellaremission component. ACKNOWLEDGMENTS
We would like to thank P. D. Serpico, C. Eckner, F. Do-nato, and M. Di Mauro for useful and enlightening dis-cussions. We also thank C. Eckner for comments on themanuscript, and C. Weniger for initial support with the
SkyFACT code. We acknowledge support from the ANRJCJC 2019 grant “GECO” (PI: F. Calore). This work hasbeen done thanks to the interactive servers of the Univ.Savoie Mont Blanc - CNRS/IN2P3 MUST computingcenter. We would also like to thank M. Gauthier-Lafayefor his work on the servers upon our multiple requests.
Appendix A: Description of the fitting algorithm
SkyFACT
For each emission model component of the fit, the γ -rayflux is modeled for a pixel p and an energy bin b by thefollowing equation:Φ pb = (cid:88) k T ( k ) p τ ( k ) p · S ( k ) b σ ( k ) b · ν ( k ) (A1)where T ( k ) p and S ( k ) b are the spatial and spectral templatesrespectively of the input component k and τ ( k ) p and σ ( k ) b are their associated spatial and spectral modulation (ornuisance) parameters. The quantity ν ( k ) represents anoverall normalization factor. More details can be foundin [40].The modeled γ -ray flux is then divided by the Fermi -LAT exposure and convolved with the
Fermi -LAT PSF toderive the expected number of γ rays per energy bin andper pixel. This model is then fitted to the Fermi -LATcount map, minimizing the following total log-likelihood:ln L = ln L P + ln L R , (A2)where L P is the standard Poisson likelihood, and L R repre-sents the regularization term which controls the variationof the modulation parameters.The optimization is performed with the L-BFGS-B(Limited memory BFGS with Bound constraints) al-gorithm [50–52], which is similar to the BFGS quasi-Newtonian method but uses less memory. These algo-rithms are developed to find the local extrema of functions As for point-like sources, they are described only by a spectraltemplate. based on Newton’s method where a second degree approx-imation is used to find the minimum function. They bothuse the Hessian inverse matrix estimation to search forthe parameters maximizing the function. The advantageof the L-BFGS-B is that the algorithm does not store adense n × n approximation to the inverse Hessian matrixbut instead stores a few vectors representing the approxi-mation implicitly. Therefore, the L-BFGS-B algorithm iswell suited for very large-scale problems, i.e. with a largenumber of parameters. In addition, the L-BFGS-B algo-rithm supports boundary conditions which is necessaryto impose non-negativity constraints on the parameters,see discussion in [40].The uncertainties of individual model parameters andcomponent fluxes are computed using a sampling method.This technique circumvents the computation of the co-variance matrix which would require a significant compu-tational time due to the large number of parameters. Itrelies on the Fisher information matrix given by: I ( θ ) ij = − (cid:28) ∂ ∂θ i ∂θ j ln L (cid:29) D ( θ ) , (A3)where θ are the model parameters and i, j define a matrixelement. The average is calculated using mock data gen-erated from the best-fit model D ( θ ). The Fisher matrixis then decomposed into triangular and diagonal matri-ces using sparse Cholesky decomposition [53] from whichsample model parameter vectors δθ are computed. Meanvalues of component fluxes and model parameters arederived from the best-fit value of θ while standard devi-ations are derived by computing model predictions for θ + δθ which are averaged over many samples. Appendix B: Statistical analysis framework
The total likelihood given by
SkyFACT is the sum ofa Poisson likelihood L P and a regularization term L R that controls the modulation parameters. To probe theextension of M31, we use the likelihood ratio (LR) teststatistic where we state two hypotheses. The null hy-pothesis H corresponds to the absence of M31, whereonly the diffuse components (IEM, IGRB) and PS areconsidered, while the alternative hypothesis H includesan emission component for M31. The LR test statistic(TS) is defined as:TS = − L H ( ˆ θ Diffuse , ˆ θ IGRB , ˆ θ PS ) L H ( ˆ θ M31 , ˆ θ Diffuse , ˆ θ IGRB , ˆ θ PS ) , (B1)where θ M31 , θ Diffuse , θ IGRB , and θ PS are the parameterscorresponding to M31, IEM, IGRB, and PS componentsrespectively. The hat refers to the values maximizingthe likelihood functions L H and L H . The parameters θ i are all constrained in R + . Typically, the LR TS isdistributed as a χ q distribution if the null hypothesis H is true, according to Wilk’s theorem [54], where q is the3number of free parameters. However, in this case, theparameters are not allowed to take on negative values. Asa result, the Wilk’s theorem does not hold anymore sincethe values of the parameters of L H are on the boundaryof the allowed parameter space. For this theorem to bevalid, the maximum likelihood estimators must find a non-boundary value where the function is differentiable [55].This problem can be solved using the Chernoff theorem[56] for nested models which combines a χ and a Diracfunction δ at 0: P ( T S ) = 2 − n (cid:34) δ (0) + n (cid:88) i =1 (cid:18) ni (cid:19) χ i ( T S ) (cid:35) , (B2)where the 2 − n term is the number of distinct ways that n energy bins can take on non-negative values. The δ function ensures that all energy bins fulfil the nonnegativity condition while the binomial coefficient (cid:0) ni (cid:1) describes the number of possible configurations of non-negative amplitudes where each configuration has a χ i distribution [57].The TS can then be interpreted in terms of standarddeviations σ . First, one needs to compute the p -value,also known as the survival function , p = P (TS > Λ) =1 − CDF, where CDF is the cumulative distribution. Thesurvival function corresponds to the probability of obtain-ing a TS value larger than the given value Λ. In otherwords, it represents the proportion of data “surviving”above a certain value. The number of σ is then derived byperforming the square root of the inverse survival function“InverseCDF” of a χ q distribution with argument p andreads as: σ = (cid:113) InverseCDF( χ , CDF[p(TS) , (cid:99) TS]) , (B3)where (cid:99) TS is the observed TS value [57, 58]. In our hy-pothesis testing ( i.e. presence of M31), when the M31spectral parameters are fixed, the number of free parame-ters q is one, corresponding to the overall normalization ν M31 . However, when the spectral parameters are free tovary, the number of free parameters equals the number ofenergy bins [58].
Appendix C: Non nested model comparison
Comparing complex, non-nested, models is quite subtlesince the number of effective parameters vary from a modelto another. The larger the number of effective parameters,the easier it is for a fit to converge. Therefore, one needsto take into account the number of effective parametersin order to assess the goodness of the models from one toanother. This is what is required to meaningfully comparedifferent spatial models for M31.In the case of non nested models, the δ - χ q mixturecannot be applied. As an alternative, model comparisoncan be performed using the Akaike Information Criterion(AIC) defined as [59]:AIC = 2 N effParam − L Data ) , (C1) where N effParam is the number of estimated effective param-eters in the model and − L P Data ) is the minimizedPoisson likelihood function on the data.The effective number of parameters of each model isderived using [40]: N effParam = N effData − N effDOF , (C2)where N effData is the effective number of data bins and N effDOF is the effective degree of freedom (DOF). The quantity N effData is derived by averaging Poisson likelihood values ofseveral sets of mock data obtained without refitting thedata: N effData = (cid:104)− L P ( θ ) (cid:105) Mock . (C3)The second quantity N effDOF , however, is defined by the av-erage of the Poisson likelihood values obtained by refittingthe mock data: N effDOF = (cid:104)− L P (ˆ θ ) (cid:105) Mock , (C4)where ˆ θ are the values of the model parameters thatminimize the likelihood function. In what follows, we usetwenty sets of mock data to derive N effData and N effDOF . Appendix D: R regression In Sec. VI, we use the R regression sum of squaresmethod to estimate the goodness of the intensity profilefit. The R regression is defined by [60]: R = SS Regression SS Total = 1 − SS Residual SS Total , (D1)where SS Regression = N (cid:88) i =1 ( ˆ Y i − ¯ Y ) (D2) SS Residual = N (cid:88) i =1 ( Y i − ˆ Y i ) (D3) SS Total = N (cid:88) i =1 ( Y i − ¯ Y ) (D4)where Y i and ˆ Y i are the observed and best-fit intensityvalues in the annulus i respectively and ¯ Y is the meanintensity. The term SS Regression quantifies how much thebest fit deviates from ¯ Y , while SS Residual estimates thediscrepancy between the model and the actual intensity,and SS Total measures the variability of the data. Thecloser SS Regression is to SS Total , the better the model,in contrast with SS Residual where the closer the valueto SS Total , the worse the model. The R value is givenbetween 0 and 1, 0 indicating a bad fit and 1 a perfect fit.4 [1] M. Ackermann et al. (Fermi-LAT), Astrophys. J. ,208 (2017), arXiv:1702.08602 [astro-ph.HE].[2] M. Di Mauro, X. Hou, C. Eckner, G. Zaharijas,and E. Charles, Phys. Rev. D99 , 123027 (2019),arXiv:1904.10977 [astro-ph.HE].[3] A. G. Riess, J. Fliri, and D. Valls-Gabaud, The Astro-physical Journal , 156 (2012).[4] Y. Joshi, A. Pandey, D. Narasimha, R. Sagar, andY. Giraud-Hiraud, Astron. Astrophys. , 113 (2003),arXiv:astro-ph/0210373.[5] K. Stanek and P. Garnavich, Astrophys. J. Lett. ,L131 (1998), arXiv:astro-ph/9802121.[6] A. R. Conn, R. A. Ibata, G. F. Lewis, Q. A. Parker, D. B.Zucker, N. F. Martin, A. W. McConnachie, M. J. Irwin,N. Tanvir, M. A. Fardal, and et al., The AstrophysicalJournal , 11 (2012).[7] P. R. Durrell, W. E. Harris, and C. J. Pritchet, Astron.J. , 2557 (2001), arXiv:astro-ph/0101436.[8] S. Holland, Astron. J. , 1916 (1998), arXiv:astro-ph/9802088.[9] I. Ribas, C. Jordi, F. Vilardell, E. L. Fitzpatrick, R. W.Hilditch, and E. F. Guinan, Astrophys. J. Lett. , L37(2005), arXiv:astro-ph/0511045.[10] A. W. McConnachie, M. Irwin, A. Ferguson, R. Ibata,G. Lewis, and N. Tanvir, Mon. Not. Roy. Astron. Soc. , 979 (2005), arXiv:astro-ph/0410489.[11] S. Abdollahi et al. (Fermi-LAT), Astrophys. J. Suppl. , 33 (2020), arXiv:1902.10045 [astro-ph.HE].[12] L. Feng, Z. Li, M. Su, P.-H. T. Tam, and Y. Chen,Res. Astron. Astrophys. , 046 (2019), arXiv:1810.10721[astro-ph.HE].[13] E. Tempel, A. Tamm, and P. Tenjes, (2007),arXiv:0707.4374 [astro-ph].[14] M. A. Fardal et al. , Mon. Not. Roy. Astron. Soc. ,2779 (2013), arXiv:1307.3219 [astro-ph.CO].[15] L. M. Widrow, K. M. Perrett, and S. H. Suyu, Astrophys.J. , 311 (2003), arXiv:astro-ph/0301144.[16] N. Evans and M. Wilkinson, Mon. Not. Roy. Astron. Soc. , 929 (2000), arXiv:astro-ph/0004187.[17] N. Evans, M. Wilkinson, P. Guhathakurta, E. Grebel, andS. Vogt, Astrophys. J. Lett. , L9 (2000), arXiv:astro-ph/0008155.[18] M. S. Seigar, A. J. Barth, and J. S. Bullock, Mon.Not. Roy. Astron. Soc. , 1911 (2008), arXiv:astro-ph/0612228.[19] A. Tamm, E. Tempel, P. Tenjes, O. Tihhonova, andT. Tuvikene, Astronomy & Astrophysics , A4 (2012).[20] E. Corbelli, S. Lorenzoni, R. Walterbos, R. Braun, andD. Thilker, Astronomy and Astrophysics , A89 (2010).[21] R. Feldmann, D. Hooper, and N. Y. Gnedin, The Astro-physical Journal , 21 (2012).[22] N. Lehner, J. C. Howk, and B. P. Wakker, The Astro-physical Journal , 79 (2015).[23] Z. Li, X. Huang, Q. Yuan, and Y. Xu, JCAP , 028(2016), arXiv:1312.7609.[24] Abdo, A. A. et al. , .[25] M. Pshirkov, V. Vasiliev, and K. Postnov, Mon. Not. Roy.Astron. Soc. , L76 (2016), arXiv:1603.07245 [astro-ph.HE].[26] C. Karwin, S. Murgia, S. Campbell, and I. Moskalenko,(2019), arXiv:1903.10533 [astro-ph.HE].[27] A. Do, M. Duong, A. McDaniel, C. O’Connor, S. Pro- fumo, J. Rafael, C. Sweeney, and W. Vera, (2020),arXiv:2012.14507 [astro-ph.HE].[28] L. Goodenough and D. Hooper, arXiv e-printsarXiv:0910.2998.[29] D. Hooper and L. Goodenough, arXiv:1010.2752.[30] D. Hooper and T. Linden, arXiv:1110.0006.[31] D. Hooper and T. R. Slatyer, arXiv:1302.6589.[32] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden,S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer,arXiv:1402.6703.[33] M. Ajello et al. (Fermi-LAT), arXiv:1511.02938.[34] F. Calore, I. Cholis, and C. Weniger, arXiv:1409.0042.[35] V. Vitale and A. Morselli, arXiv e-prints (2009),arXiv:0912.3828.[36] K. N. Abazajian and M. Kaplinghat, arXiv:1207.6047.[37] A. Boyarsky, D. Malyshev, and O. Ruchayskiy,arXiv:1012.5839.[38] W.-C. Huang, A. Urbano, and W. Xue, arXiv e-printsarXiv:1307.6862.[39] C. Eckner et al. , Astrophys. J. , 79 (2018),arXiv:1711.05127 [astro-ph.HE].[40] E. Storm, C. Weniger, and F. Calore, JCAP , 022(2017), arXiv:1705.04065 [astro-ph.HE].[41] .[42] M. Di Mauro and F. Donato, Phys. Rev. D , 123001(2015), arXiv:1501.05316 [astro-ph.HE].[43] J. Ma, Chinese Physics Letters , 1420-1422 (2001).[44] P. D. Group et al. , .[45] M. Ackermann et al. (Fermi-LAT), Astrophys. J. Suppl. , 32 (2018), arXiv:1804.08035 [astro-ph.HE].[46] B. S. Acharya et al. (CTA Consortium), Science with theCherenkov Telescope Array (WSP, 2018) arXiv:1709.07997[astro-ph.IM].[47] X. Bai et al. , (2019), arXiv:1905.02773 [astro-ph.HE].[48] K. P. Burnham and D. R. Anderson,