Dark Matter Interpretation of the Fermi-LAT Observations Toward the Outer Halo of M31
Chris Karwin, Simona Murgia, Igor Moskalenko, Sean Fillingham, Anne-Katherine Burns, Max Fieg
DDark Matter Interpretation of the
Fermi -LAT ObservationsToward the Outer Halo of M31
Christopher M. Karwin,
1, 2, ∗ Simona Murgia, † Igor V. Moskalenko, ‡ Sean P. Fillingham,
4, 2
Anne-Katherine Burns, and Max Fieg Department of Physics and Astronomy, Clemson University, Clemson, SC, USA Department of Physics and Astronomy, University of California, Irvine, CA, USA Hansen Experimental Physics Laboratory and Kavli Institute for ParticleAstrophysics and Cosmology, Stanford University, Stanford, CA, USA Department of Astronomy, University of Washington, Seattle, WA, USA (Dated: October 20, 2020)An excess γ -ray signal toward the outer halo of M31 has recently been reported. Although otherexplanations are plausible, the possibility that it arises from dark matter (DM) is valid. In this workwe interpret the excess in the framework of DM annihilation, using as our representative case WIMPDM annihilating to bottom quarks, and we perform a detailed study of the systematic uncertaintyin the J -factor for the M31 field. We find that the signal favors a DM particle with a mass of ∼ ∼ × − − × − cm s − . This high uncertainty is due totwo main factors, namely, an uncertainty in the substructure nature and geometry of the DM halosfor both M31 and the Milky Way (MW), and correspondingly, an uncertainty in the contributionto the signal from the MW’s DM halo along the line of sight. However, under the conditions thatthe minimum subhalo mass is < ∼ − M (cid:12) and the actual contribution from the MW’s DM haloalong the line of sight is at least ∼
30% of its total value, we show that there is a large overlap withthe DM interpretations of both the Galactic center (GC) excess and the antiproton excess, whilealso being compatible with the limits for the MW dwarf spheroidals. More generally, we summarizethe results from numerous complementary DM searches in the energy range 10 GeV −
300 GeVcorresponding to the GC excess and identify a region in parameter space that still remains viablefor discovery of the DM particle.
I. INTRODUCTION
Observational evidence for dark matter (DM) in M31comes from measurements of its rotational velocitycurve [1–5]. These observations provide coarse-grainedproperties of the DM distribution near the central re-gions of the halo where the galaxy resides. With theexisting data, the fine-grained structure of DM and itsdistribution outside of the galaxy is primarily inferredfrom simulated halos. Within the standard cosmologi-cal paradigm, M31’s DM halo is expected to extend wellbeyond the galactic disk, and it is also expected to con-tain a large amount of substructure. However, there iscurrently a high level of uncertainty regarding the exactnature of the halo properties, i.e. the geometry, extent,and substructure content, especially on galactic scales [6–32].Due to its mass and proximity, the detection sensitiv-ity of M31 to DM searches with γ -rays is competitivewith the Milky Way (MW) dwarf spheroidal galaxies,particularly if the signal is sufficiently boosted by sub-structures [33–38]. Moreover, M31 is predicted to be thebrightest extragalactic source of DM annihilation [39, 40]. ∗ [email protected] † [email protected] ‡ [email protected] A detailed study of the γ -ray emission observed to-wards M31’s outer halo has recently been made inRef. [32]. In that study evidence is found for an excesssignal that appears to be distinct from the conventionalMW foreground, having a total radial extension upwardsof ∼ γ -ray emission from M31’s inner galaxy has also beendetected, but the exact nature of the emission still re-mains an open question, as the morphology of the signaldoesn’t appear to trace regions rich in gas and star forma-tion [32, 41–47]. On the other hand, the total γ -ray lumi-nosity is found to be in general agreement with the well-known scaling relationship between the γ -ray luminos-ity and infrared luminosity (8–1000 µ m) for star-forminggalaxies [48]. Ultimately, a better determination of the γ -ray signal from M31’s inner region is still needed, whichwill require a refinement of the underlying gas maps (H i )used to model the Galactic foreground emission, as thecurrent maps may be holding a fraction of gas that actu-ally resides in the M31 system [32]. The Doppler-shiftedvelocity of the gas, together with the Galactic rotation a r X i v : . [ a s t r o - ph . H E ] O c t curve, is used to separate the MW and M31 gas. The un-certainty arises from two main conditions. First, there isa partial overlap of the rotational velocities for M31 andthe MW. Second, M31 is at a fairly high latitude wherethere is an increased uncertainty in the rotational speedof the MW gas, which is measured in the Galactic disk.In this work we interpret the excess γ -ray emission ob-served towards M31’s outer halo in the framework of DMannihilation. We consider WIMP (i.e. weakly interactingmassive particle) DM, and focus the analysis on the un-certainties associated with the properties of the DM halo.Moreover, we consider a realistic observational perspec-tive, in which the line of sight towards M31’s outer DMhalo naturally extends through a similar DM halo aroundthe MW. In general, this is not directly accounted forwhen modeling the MW foreground γ -ray emission, andcan significantly impact the results.The paper is organized as follows. In Section II wegive a qualitative description of M31’s outer halo. InSection III we present the M31 data, DM fit, and an-alytical J -factor calculations. In Section IV we presentresults for our best-fit models, and we consider these re-sults in the context of the Galactic center (GC) excess,and more generally, in the context of the current statusof DM indirect detection. In Section V we conclude. Ad-ditional details for the complementary DM searches weconsider are given in Appendix A. II. M31’S OUTER HALO
For observations of γ -ray emission arising from DMannihilation towards M31’s outer halo, the total signalwould ostensibly contain emission from the MW’s DMhalo along the line of sight, emission from the local fil-amentary structure connecting the MW and M31, andemission from the entire DM halo of M31, plus any sec-ondary emission (from M31 and the MW). For the MWhalo, a DM signal should be pretty bright, but since theobservation occurs from within the halo, the emissioncan be easily confused with the isotropic component (andother components of the MW interstellar emission model(IEM)). For M31, we observe the entire halo from theoutside, and therefore we see the total integral signal.Thus M31 is advantageous for halo searches with γ -raysbecause it breaks the observational degeneracy.Figure 1 provides a qualitative description of M31’souter halo, including an accounting of some notablestructures along the line of sight that may provide hintsof the DM distribution. The γ -ray counts map (shown inblack and white) is from Ref. [32]. The bright emissionalong zero degree latitude is the plane of the MW. Thesize of M31’s DM halo is indicated with a dashed cyan cir-cle, which corresponds to a projected radius of 300 kpc,for an M31-MW distance of 785 kpc. The dash-dot lime-green circle shows the outer boundary of the sphericalhalo (SH) region, which we use for the DM fit, as dis-cussed in Section III. M31’s satellite population is shown FIG. 1. The line of sight looking towards M31’s outer halo.The size of M31’s DM halo is indicated with a dashed cyancircle, which corresponds to a projected radius of 300 kpc, foran M31-MW distance of 785 kpc. The dash-dot lime-greencircle shows the outer boundary of the SH region ( r tan = 117kpc), which we use for the DM fit. M31’s population of satel-lite galaxies is shown with red open circles. M33 can be seenin the lower left corner. Also plotted are some notable gasclouds in the region, namely, the M31 cloud (orange regionsurrounding the M31 disk), Wright’s cloud (WC), and Com-plex H. See text for more details. with open red circles. A subset of the satellites in M31(which are thought to reside within DM substructures)are known to be positioned within a large thin plane (TheGreat Plane of Andromeda, GPoA); and likewise, a sub-set of the MW satellites are known to be part of a largeplanar structure as well (The Vast Polar Structure of theMilky Way) [49–55]. In addition, the satellite system ofM31 is highly lopsided, as about 80% of its satellites lieon the side closest to the MW [51, 56]. For members ofthe GPoA, those to the north of M31 recede from us, andthose to the south of M31 move toward us, in the planeof rotation.Also shown in Figure 1 are two notable, highly ex-tended gas clouds in the direction of M31, namely, Com-plex H [8, 57–59] and the M31 cloud [8, 60]. The gascontours show H i emission from the HI4PI all-sky sur-vey (based on EBHIS and GASS) [61]. The M31 cloud isa highly extended lopsided gas cloud centered in projec-tion on M31, originally reported in Ref. [8]. It remainsuncertain whether the M31 cloud resides in M31 or theMW. Most recently Ref. [60] has argued that M31’s diskis physically connected to the M31 cloud. If at the dis-tance of M31 ( ∼
785 kpc) the total gas mass is estimatedto be ∼ –10 M (cid:12) . Complex H can be seen towardthe top of M31’s DM halo. The distance of Complex Hfrom the MW is uncertain, although its likely distancehas been estimated to be ∼
30 kpc from the GC, whichcorresponds to the cloud having a diameter of about ∼ i mass of ∼ M (cid:12) [8, 58, 59]. ComplexH does not appear to contain any stars, and it has beenpostulated to be either a dark galaxy of the Local Groupor an example of a cold accretion flow [59].Figure 1 also shows H i emission contours correspond-ing to M33. γ -ray emission from M33 has recently beendetected [32, 62, 63], making it the only extragalacticsatellite galaxy to be detected in γ rays. The total H i mass of the M33 disk is ∼ M (cid:12) . The hook-shapedgas cloud to the right of M33 is Wright’s cloud, first re-ported in Ref. [64]. The distance of Wright’s cloud re-mains uncertain [11]. The H i mass of Wright’s cloud atthe distance of M33 is ∼ . × M (cid:12) [65]. Although nocontours are shown, we note that below M33 is “the darkcompanion to M33”, which is another highly extendedgas cloud originally reported in Ref. [66], and labeled asa compact high-velocity cloud. If at the distance of M33,Ref. [65] estimates the H i mass to be ∼ M (cid:12) , and thesize to be ∼ . × . γ rays (aside from M33). Forthe gas clouds, any γ -ray emission would depend on theactual location of the cloud, along with the CR densityin the region. To investigate this in depth would requirea detailed modeling which is beyond the scope of thisanalysis. III. ANALYSISA. Gamma-Ray Data for M31
To determine whether the excess γ -ray emission ob-served towards M31’s outer halo is consistent with a DMinterpretation, we employ the best-fit γ -ray spectra fromRef. [32]. In that study M31’s halo is characterized usingthree symmetric components centered at M31 labeled as:inner galaxy (IG; r ≤ . ◦ ), spherical halo (SH; 0 . ◦ > r ≤ . ◦ ), and far outer halo (FOH; r > . ◦ ). For anM31-MW distance of 785 kpc, the IG, SH, and FOH cor-respond to projected radii of 5.5 kpc, 117 kpc, and ∼ γ -ray emission from standardastrophysical processes. The FOH component overlapswith the MW plane at the top of the field, which sig-nificantly complicates the interpretation of the emissionfrom this region. In addition, properly modeling the FOH will require a thorough treatment of secondary emissionfrom DM, which we leave for a future study.Two different fit variations were performed in Ref. [32]to determine the spectrum of the SH component. In themain variation (full) the entire template was used. In analternative variation (north and south) the template wasseparated into north and south components. In this casethe spectral parameters for the two halves are allowedto vary independently, although they are fit simultane-ously. This results in three different determinations of thespectrum, which we label as spherical halo (SH), spher-ical halo north (SHN), and spherical halo south (SHS).We use these variations to quantify the systematic un-certainty of the signal related to modeling the MW fore-ground, which differs in the two regions.It is important to emphasize that the line of sight to-wards M31 extends through the MW DM halo, in addi-tion to the M31 DM halo. However, the potential γ -raycontribution from the MW component is not explicitlyaccounted for when determining the M31 contribution.Some of the MW halo component would likely be at-tributed to the isotropic component, as well as to theother components of the IEM; however, it is unclear theextent to which this would occur. This is partly dueto the fact that the absorption of a MW DM halo sig-nal by other MW components in large part depends onthe actual halo geometry and substructure content in thedirection of the M31 field. Thus the spectra for the M31-related components from Ref. [32] contain the total excessemission along the line of sight, which may also includesome significant contribution from the MW’s extendedDM halo. This is taken into account in our J -factor cal-culations. B. Dark Matter Fit
As our representative DM model we consider annihila-tion into bottom quarks. This channel has been shown toprovide a good fit to the γ -ray GC excess. The DM spec-tra are obtained from PPCC 4 DM ID [67, 68], and theyinclude electroweak corrections. We scan DM massesfrom 10 GeV to 280 GeV, using a 10 GeV spacing.The γ -ray flux for DM annihilation is given by d Φ dE = < σ f v > πηm χ dN fγ dE J, (1)where < σ f v > is the velocity averaged annihilationcross-section for final state f , m χ is the DM mass, η = 2(4) for self-conjugate (non-self-conjugate) DM, dN fγ /dE is the number of γ -ray photons for annihilation into finalstate f , and J is the astrophysical J -factor, which will be available at DM Mass
DM fit to SHDM fit to SHNDM fit to SHS Energy [MeV] N × (E × dN f /dE) [MeV cm s ] bb DM Flux
DM fit to SHDM fit to SHNDM fit to SHSSH dataSHN dataSHS data
FIG. 2.
Left panel : ∆ χ profile for the three different fit variations: spherical halo (SH): solid black curve; spherical halonorth (SHN): dash-dot turquoise curve; spherical halo south (SHS): dashed grey curve. The light grey dotted lines show the1, 2, and 3 sigma contour levels, for 1 degree of freedom. Right panel : Best-fit spectra overlaid to the corresponding data.Arrows give the 1 σ upper limits. discussed in Section III C. In general Eq. (1) is summedover all final states f . In this analysis we use η = 2.By multiply each side of Eq. (1) by the energy squaredwe obtain units of MeV cm − s − : E d Φ dE = < σ f v > πηm χ ( E dN fγ dE ) J. To fit to the γ -ray data we freely scale the quantity inparentheses by a normalization factor N , using a χ fit.This then implies: N = < σ f v > πηm χ J. (2)The M31 data contains upper limits which need to beaccounted for in the fit procedure. For n measurementsof x i with uncertainties σ i and m upper limits with x j For the best-fit models the corresponding annihilationcross-section is calculated using Eq. (2). This requiresknowledge of the J -factor, which is the greatest uncer-tainty in the analysis. The J -factor characterizes thespatial distribution of the DM, and is given by the inte-gral of the mass density squared, over the line of sight.When describing the DM distribution as an ensemble ofdisjoint DM halos, the J -factor is: J = (cid:88) i (cid:90) ∆Ω d Ω (cid:90) LoS dsρ i ( r i ( s, n )) , (7)summed over all halos in the line of sight (LoS), where ρ i ( r ) is the density distribution of halo i , and r i ( s, n ) isthe position within that halo at LoS direction n and LoSdistance s . J -factors determined from these spherically-averagedprofiles are an underestimate of the total J -factor becauseof the effect of the non-spherical structure. This under-estimate is typically encoded with a boost factor ( B ).To calculate J -factors we use the CLUMPY code [71–73]. For a detailed discussion of the boost factor cal-culation see the CLUMPY papers/website, as well asRefs. [9, 19, 25, 26, 30] and references therein. Here wesummarize the key points. The main parameters for theboost factor are the following: ◦ minimum subhalo mass ◦ mass-concentration relationship ◦ subhalo mass function (index and normalization),i.e. the number of subhalos per volume in a givenmass range ◦ mass distribution of subhalos ◦ distribution of subhalos in the main haloSince the γ -ray flux from DM annihilation scales asthe square of the DM density, the effect of substructureis very important for indirect detection, as it provides aboost to the total flux. The flux enhancement is mostsignificant for larger halos, since they enclose more levelsof hierarchical formation. The size of the smallest DMsubhalo is determined by the free streaming scale of theDM particles [26, 74, 75]. This depends on the specificparticle physics and cosmological models, and in generalit is highly uncertain. In this study we consider minimumsubhalo masses in the range M min = 10 − − M (cid:12) .The lower limit is typically expected for thermal WIMPDM with a mass of ∼ 100 GeV [74], and the upper limitreflects the typical resolution power of DM simulations.The concentration parameter c ∆ , at a given character-istic overdensity ∆, can be defined as c ∆ = R ∆ r − , (8)where R ∆ is the radius of the DM halo correspondingto the overdensity ∆, and r − is the position where theslope of the DM density profile reaches − 2. The boostfactor is highly sensitive to the concentration parameter,as it scales as the concentration to the third power [71–73]. In general the concentration is a function of halomass and redshift. In the top panel of Figure 3 we plotdifferent determinations of the concentration-mass rela-tion at z=0. The solid lines (black, purple, magenta, andred) are from Ref. [30], which is based on two N-bodycosmological simulations of MW-sized haloes: VL-II [20] available at https://clumpy.gitlab.io/CLUMPY/ M [M ] c c -M Relations Wang+19Wang+19 (with free-streaming cutoff)Ishiyama+19Moline+17 (0 < x sub < 0.1)Moline+17 (0.1 < x sub < 0.3) Moline+17 (0.3 < x sub < 1.0)Moline+17 (1.0 < x sub < 1.5)Sanchez-Conde+14Bullock+01 r [kpc] ( r ) [M kpc ] DM Density Profile Karwin+19, NFWTamm+12, NFWTamm+12, BurkertTamm+12, Einasto M vir [M ] Boost Boost Factor Ishiyama+19 high (-6,2.0,0.35)Ishiyama+19 low (semianalytic model)Moline+17 (-6,2.0,0.35)Moline+17 (-6,1.9,0.19)Sanchez-Conde+14 (-6,2.0,0.35)Sanchez-Conde+14 (-6,1.9,0.19)Gao+12 (PL extrapolation)Bullock+01 (-6,2.0,0.35)Bullock+01 (-6,1.9,0.19)Bullock+01 (6,1.9,0.12) FIG. 3. Top panel: Concentration-mass relations fromRefs. [30, solid black, purple, magenta, and red], [75, dot-ted yellow], [76, dash-dot grey], [9, dashed green], and [77,solid and dashed cyan]. Middle panel: Different DM den-sity profiles for M31. The region bounded by the red dashedlines corresponds to the SH. Bottom panel: Mass depen-dence of the boost factor for different parameters. The namein the legend specifies the model of the concentration-massrelation, and in parentheses the numbers give (in order) thepower of the minimum subhalo mass, the PL index of thesubhalo mass function, and the fraction of the halo resolvedin substructure. The red dashed lines correspond to the massrange for M31 and the MW. and ELVIS [27]. These results summarize some of themain properties of the concentration parameter; namely,for a given halo the concentration decreases with increas-ing radius, and the concentration of subhalos is higherthan that of field halos. In particular, the solid lines inFigure 3 are for different radial bins defined in terms of x sub ≡ R sub /R ∆ . The solid black line is calculated out-side of the virial radius, and it gives an approximation forfield halos (see [30] for further details). For simplicity, inour benchmark model we use the relation from Ref. [9],plotted with a dashed green line in the top panel of Fig-ure 3. As can be seen, this serves as a good intermediatemodel between the different estimates. Note that we havealso tested the model from Ref. [76] and the results arequalitatively consistent.The boost factor also depends on the subhalo massfunction, which specifies the number of subhalos at agiven mass. This function is given by a simple powerlaw (PL), having an index of ∼ − . − . J -factor forboth M31 and the MW, we vary the index of the subhalomass function ( α ) and the fraction of the halo resolved insubstructure ( f sub ) in the ranges 1 . − . . − . γ -ray data is performed. The solid black curve is fromRef. [32], and the other curves are from Ref. [79]. Forour J -factor calculations we test two profiles. We usethe NFW profile from Ref. [32], which has correspondinghalo properties of R vir = 210 kpc, R s = 18 . ρ s = 2 × M (cid:12) kpc − . In CLUMPY this correspondsto the kZHAO profile with parameters α, β, γ = 1,3,1.We also use the Einasto profile from Ref. [79], which hasthe corresponding halo properties of R vir = 210 kpc, R s = 178 kpc, and ρ s = 8 . × M (cid:12) kpc − . InCLUMPY this corresponds to the kEINASTO N profilewith the parameter n=6. The overdensity factor is set to∆ = 200. We use an M31-MW distance of 785 kpc.Other major uncertainties in the boost factor calcu-lation are the spatial distribution of subhaloes in themain halo, as well as the mass distribution of the sub-haloes themselves. We assume that the density profileand the spatial distribution of the subhaloes are the sameas the density profile of the main halo for both the NFWand Einasto distributions. Note that both the spatialdistribution of subhaloes and their density profiles havebeen found to prefer an Einasto distribution comparedto an NFW, although both profiles provide a good fit(see [19] and references therein). Additionally, it’s foundthat within ∼ 25 kpc from the center of MW-sized halosthere is a depletion of the subhalo population due to tidaldisruption from the galactic disk [31].In principle each DM halo of a given mass is a hier-archical structure, so that even subhalos have subhalosthemselves. For simplicity we set the number of substruc-ture levels to 2. We have also tested including higher substructure levels, but we find that they do not make asignificant difference for our J -factor calculations, as hasbeen previously found [30].The bottom panel of Figure 3 shows the mass de-pendence of the boost factor for different choices of theminimum subhalo mass, the subhalo mass function, andthe fraction of the halo mass resolved in substructure.Within the uncertainties we have considered, the overallboost factor ranges from ∼ D. Halo Geometry Another important systematic uncertainty for deter-mining the J -factor for the M31 field is the halo geom-etry, for both M31 and the MW. Indirect DM searchestypically assume spherical symmetry for the halo shape,however, in the standard DM paradigm (ΛCDM), DMhalos are expected to be very non-spherical, and in fact,spherical halos are rare (see [14] and references therein).For the MW, numerous studies have been done to in-fer the DM halo geometry, but differing conclusions havebeen reached. The halo has been found to be spheri-cal [80], prolate [24, 81, 82], oblate [83], triaxial (includ-ing the so-called “Gaia sausage”) [23, 84, 85], and evenlopsided [22]. Further complicating the matter is that thehalo geometry may have a radial dependence [86, 87].Moreover, it’s found in both simulations and observa-tions that for galaxy pairs (similar to M31 and the MW)the halos tend to bulge toward their respective part-ners [51, 56].In general the halo geometry can be described with anellipsoid, with the axes a, b, and c. The shape is charac-terized by the axis ratios, with the normalization condi-tion abc = 1 (see the CLUMPY code for more details).For describing the MW halo, the a-axis corresponds tothe Galactic x-axis (connecting the Sun to the Galacticcenter), the b-axis corresponds to the Galactic y-axis,and the c-axis corresponds to the Galactic z-axis (per-pendicular to the Galactic plane). We use the referencescited above to calculate J -factors for different MW halogeometries. Note that we also consider a triaxial halogeometry modeled after the Gaia sausage. Although theevidence indicates that this structure may be a subdomi-nant component of the halo, for simplicity we test a moreextreme scenario where the entire halo follows this geom-etry. Figure 4 shows the three main halo shapes that wetest, and the specific axis ratios for all geometries aresummarized in Table I.In the top panel of Figure 5 we show the J -factor ratio( J/J Sph ) for the Einasto high DM model, where J is forthe alternative geometry, and J Sph is for the sphericalhalo. The ratio range for all DM models is given in Ta-ble II. We find that at most the halo shape may increaseor decrease the MW J -factor (with respect to spherical Spherical Halo Gev cm sr Prolate Halo Gev cm sr Oblate Halo Gev cm sr FIG. 4. MW J -factors for four different geometries, as indicated above each map. Maps are shown in Galactic coordinateswith a Mollweide projection. The corresponding axis ratios are given in Table I. For the prolate halo q=1.67, and for theoblate halo q=0.6. The color scale ranges from the minimum halo value to 1/10 the maximum halo value. The DM model is”Einasto high” from Table II. Note that these particular maps don’t show individually resolved substructures, although theyare included in the analytical model.TABLE I. MW Halo GeometryHalo Geometry Axes (a,b,c)Spherical 1, 1, 1Prolate (q=1.67) 0.84, 0.84, 1.41Prolate (q=1.25) 0.93, 0.93, 1.16Oblate (q=0.4) 1.36, 1.36, 0.54Oblate (q=0.6) 1.19, 1.19, 0.71Oblate (q=0.8) 1.08, 1.08, 0.86Triaxial 0.67, 1.34, 1.113Triaxial (Gaia Sausage, α = 70 ◦ ) 1.38, 1.06, 0.69Note: The axes are normalized so that abc=1. In general,prolate halos have a=b < c, and oblate halos have a=b > c. Forconvenience we also give the ratio q=c/a. The specific axisratios come from the literature, as discussed in the text. Forvisualization purposes, the different geometries are plotted inFigure 4. geometry) by factors of 2.29 and 0.34, respectively.To test how the MW J -factor varies with Galactic lat-itude we repeat the calculations with the line of sightcentered at latitudes of − ◦ and 0 ◦ , with longitude =121 ◦ . Note that b = − ◦ corresponds to the region usedin Ref. [32] for tuning the isotropic spectrum, which werefer to as the tuning region (TR). Results for this test areshown in the middle panel of Figure 5 (for the Einastohigh model), where we plot the J -factor ratio with re-spect to the value obtained in the TR. In all cases agradient can be seen, with the amplitude of the variationdependent on the halo geometry. This is even true for aspherical halo, due to our position in the Galaxy at ∼ J -factors for the spher-ical and prolate halos decrease by a minimum factor of0.77. Alternatively, the J -factors for the oblate and tri-axial (Gaia sausage) halos increase by a maximum factorof 1.38. Since Ref. [32] tunes the isotropic spectrum ina region below the M31 field (consistent with l = − ◦ ),these results show that it is not necessarily the case thatthe MW DM halo component would be fully absorbed by the isotropic template. Moreover, even a gradient of ∼ − 40% (as is found in the gradient calculation) wouldbe a significant contribution to the total J -factor for theM31 field.We also test how the J -factor depends on the M31 halogeometry, with the main goal of estimating the full uncer-tainty range. For simplicity we test two different geome-tries. In each case the minor-to-major axis ratio is 0.4(with a > b=c). This represents a highly flattened halo,but it has also been found for M31 in particular [18]. Wetest two different orientations, one with the major axispointing along the line of sight connecting M31 and theMW (x-axis), and the other with the major axis pointingperpendicular to the line of sight (y-axis), running fromleft to right in the field of view. Note that results for thez-axis orientation are similar to those of the y-axis orien-tation. The bottom panel of Figure 5 shows the ratio ofthe J -factor for these different geometries compared toa spherical geometry (for the Einasto high model). Theuncertainty range for all DM models is given in Table II.The M31 halo geometry introduces an uncertainty in therange 0 . − . 32, where the increase is seen for the majoraxis aligned with the x-axis and the decrease is seen forthe major axis aligned along the perpendicular axes. E. J -Factor Uncertainty from the Milky WayForeground In the context of the J -factor uncertainty from the MWforeground, we consider two extreme cases. For case I weassume that none of the MW halo signal along the lineof sight has been absorbed by the isotropic component(and other components of the IEM), and thus the total J -factor is the sum of the J -factors for the MW and M31.For case II we assume that the MW halo signal along theline of sight has been completely absorbed, and so thetotal J -factor is due only to M31. In actuality, if theobserved excess is in fact related to DM then the truecase is likely somewhere between the two extremes. TABLE II. J-Factors and Cross SectionsModel α sub f sub M min [M (cid:12) ] J MW ( × )[GeV cm − ] J M31 ( × )[GeV cm − ] J/J Sph (MW) J/J Sph (M31) J MW /J TR < σv > I ( × − )[cm s − ] < σv > II ( × − )[cm s − ] Einasto high 2.0 0.35 10 − − − − J -factors are integrated over the spherical halo component (0 . ◦ to 8 . ◦ ). The largest subhalo mass is taken to be 10%the mass of the host halo. The calculations include 2 levels of substructure. For the M31 NFW profile R vir = 210 kpc, R s = 18 . ρ s = 2 . × M (cid:12) kpc − . For the M31 Einasto profile R vir = 210 kpc, R s = 178 kpc, and ρ s = 8 . × M (cid:12) kpc − .The MW profiles have the same parameters except we use the local DM density ρ (cid:12) = 0 . cm − , with a solar distance R (cid:12) = 8 . J -factor due to the halo geometry,with respect to a spherical halo ( J Sph ). Column 9 shows the J -factor gradient (low, high) with respect to the tuning region(TR) used in Ref. [32], which is centered at b = − ◦ . There are two values for each cross-section; the first value results fromfitting to the SH data (SHN data is very similar), and the second value results from fitting to the SHS data. Subscript I onthe cross-section indicates case I, where J total = J M31 + J MW , and subscript II on the cross-section indicates case II, where J total = J M31 . Corresponding curves are plotted in Figure 6. F. Total J-Factor Uncertainty Figure 6 shows the different J -factors as a function ofradial distance from the center of M31. The grey bandis the J -factor uncertainty for M31 from this work. Thepurple band is the J -factor uncertainty for the MW fromthis work. The markers are the M31 calculations forthe NFW (squares) and Einasto (circles) profiles, withthe boost factor, corresponding to the values in Table II.The dash-dot lines towards the bottom show the smoothM31 profiles corresponding to the markers. As can beseen, the smooth profiles are anti-correlated to the to-tal profiles, i.e. as the boost factor increases, the fractionof DM resolved in substructure also increases, and thefraction of the smooth DM component decreases. Thesolid curves are independent calculations for M31 fromRef. [32] (extending to 14 degrees) and Ref. [46] (extend-ing to 10 deg). Likewise the dashed lines are independentcalculations for the MW. As can be seen, there is goodconsistency between the different estimates. Our result-ing models are summarized in Table II.We note that Ref. [39] reports an M31 J -factor(integrated within the scale radius) of (6 . +7 . − . ) × GeV cm − , corresponding to a boost factor of 2.64and a scale radius of 2 . ◦ . The uncertainty in their cal-culation comes from the uncertainty in M vir and c vir .Their boost factor is comparable to our low and midmodels (with an NFW profile). When integrating overthe same scale radius, we obtain J -factor values in therange 2 . × − . × GeV cm − , in agreementwith the values reported in Ref. [39]. IV. RESULTS We calculate annihilation cross-sections using Eq. (2)with the values obtained from following the proceduredescribed in Sec. III, and results are given in Table II.There are two values for each cross-section; the first valueresults from fitting to the SH data (SHN data is very sim-ilar), and the second value results from fitting to the SHSdata. In Figure 7 we plot the corresponding best-fit DMparameters. The red data points correspond to case I,for which J = J MW + J M31 . The coral data points arefor case II, for which J = J M31 . The filled squares are forthe SH fit, and the open squares are for the SHS fit. Theresults for the SHN fit are very similar to the SH, and sowe do not include them here. Note that the error bars inthe cross-section assume that the minimum subhalo massis 10 − M (cid:12) , and they include the uncertainty due to thehalo geometry outlined in Sec. III. We compare the datapoints from M31’s outer halo to numerous complemen-tary targets for indirect DM searches. Details for all ofthe overlays are given in Appendix A.Broadly speaking, contours for the GC excess areshown in black, and contours for the antiproton excessare shown in teal. As can be seen, there is a ratherlarge range in the different determinations. This is dueto the different assumptions that are made in each anal-ysis. Generally speaking, these results can be interpretedcollectively as defining the currently explored systematicuncertainties in the respective signals. In the case ofthe GC excess, the uncertainty range in the cross-sectionspans roughly 1.5 orders of magnitude. This is because 50 40 30 20 10 0Galactic Latitude [ ]0.000.250.500.751.001.251.501.752.00 J / J T R Milky Way sphericalprolate (q=1.67)prolate (q=1.25) oblate (q=0.4)oblate (q=0.6)oblate (q=0.8) triaxialtriaxial (Gaia sausage) FIG. 5. Top: Ratio of the J -factor ( J ) for different MWhalo geometries compared to a spherical halo ( J Sph ), for anEinasto density profile. Middle: Gradient ratio for the J -factor calculated with the line of sight centered at three dif-ferent Galactic latitudes (with l = 121 ◦ ). The ratio is calcu-lated with respect to a latitude of b = − ◦ ( J TR ), which iscomparable to the region used for tuning the isotropic spec-trum in Ref. [32]. The middle data points at l = − − . ◦ correspond to the M31 field. In all cases the J -factors areintegrated over the region 0 . ◦ to 8 . ◦ , using the Einasto highmodel from Table II. Bottom: Ratio of the J -factor ( J ) fordifferent M31 halo geometries compared to a spherical halo( J Sph ), for an Einasto density profile. the GC excess is only a small fraction of the total emis-sion in the region, and thus it has a strong dependenceon the treatment of the IEM, which in general is difficultto accurately model due to the complexity of the GC re-gion. Moreover, the inferred DM parameters also havea strong dependence on the halo assumptions, such asthe local DM density, which may span between ∼ J-Factor [GeV cm sr ] Angle [deg] J-Factor [GeV cm sr ] Angle [deg] SHIG FOH FM31 J-Factor K+19, M31 NFWK+19, M31, NFW + Substructure (Mid)K+19, M31, NFW + Substructure (High)K+19, MW, NFW + Substructure (Mid)M+19 (M31 Adiabatic, Max)M+19 (M31 Einasto, Med) M+19 (M31 Burkert, Min)M+19 (MW NFW, High)M+19 (MW NFW, Med)M+19 (MW NFW, Min)M+19 (MW NFW, Smooth)M31 Einasto High M31 Einasto MidM31 Einasto LowM31 NFW HighM31 NFW High SmoothM31 NFW Mid M31 NFW Mid SmoothM31 NFW LowM31 NFW Smooth LowM31 Band (this work)MW Band (this work) FIG. 6. J -factors for M31 and the MW. The grey band is the J -factor uncertainty for M31 from this work. The blue bandis the J -factor uncertainty for the MW from this work. Themarkers are the M31 calculations for the NFW (squares) andEinasto (circles) profiles, with the boost factor. Parametersfor the different variations are given in Table II. The solidcurves are independent calculations for M31 from Ref. [32](extending to 14 degrees) and Ref. [46] (extending to 10 deg).Likewise the dashed lines are independent calculations for theMW. The dash-dot lines towards the bottom show the smoothM31 profiles corresponding to the markers. The vertical dot-ted red lines show the boundaries of M31’s IG, SH, and FOH(the fit is performed over the SH). GeV / cm [85, 88]. In the case of the antiproton excess,Refs. [89, 90] report detection contours, whereas Ref. [91]takes a less optimistic view, reporting upper limits (al-though the limits still clearly show an anomaly aroundthe signal region).Another important constraint is the upper limits fromthe MW dwarfs. Here too there is a fairly large un-certainty range. Compared to the limits reported inRef. [92], the latest limits from Ref. [93] are less con-straining. These limits of course have a strong depen-dence on the assumptions made for the J -factors, andby employing semi-analytic models of DM subhalos toderive realistic satellite priors on the J -factor (for theultrafaint dwarfs), Ref. [94] has recently shown that thelimits may be even weaker, by a factor of ∼ . ◦ disk, as determined from the emission it-self. Upper limits for a DM signal are then calculated inaddition to the disk. While this is definetly a very con-0 FIG. 7. DM parameter space. The red and coral data points are for M31’s outer halo. The red data points correspond to caseI, for which J = J MW + J M31 . The coral data points are for case 2, for which J = J M31 . The filled squares are for the SHfit, and the open squares are for the SHS fit. The results for the SHN are very similar to the SH, and so we don’t include themhere. Note that the error bars in the cross-section assume that the minimum subhalo mass is 10 − M (cid:12) , and they include theuncertainty due to the halo geometry. Contours for the GC excess are shown in black, and contours for the antiproton excessare shown in teal. Numerous limits from other targets are also overlaid, including the MW satellites shown with purple curves,and M31’s inner galaxy shown with a red curve. See Section IV for more details, as well as Appendix A. servative choice to make, it is by no means preferred, asthe γ -ray emission from M31’s inner galaxy has actuallybeen found to not correlate with regions rich in gas andstar formation.The data points for M31’s outer halo have a large over-lap with the DM interpretations of both the GC excessand the antiproton excess, while also being compatiblewith the limits from the MW dwarfs. However, this re-quires that the J -factor be towards the higher end of theuncertainty range. Correspondingly, this has two mainimplications. First, the minimum subhalo mass must be < ∼ − M (cid:12) . Second, the signal must have some contri-bution from the MW’s DM halo along the line of sight,i.e. the J -factor must correspond to case I, as it cannotbe due to M31 alone. V. SUMMARY, DISCUSSION, ANDCONCLUSION An excess γ -ray signal towards the outer halo of M31has recently been reported [32]. In this work we inter- pret the excess in the framework of DM annihilation. Asour representative case we use WIMP DM annihilatingto bottom quarks, and we fit the DM mass and annihi-lation cross-section to the observed γ -ray spectra fromRef. [32]. In that study M31’s halo is characterized usingthree symmetric components centered at M31, namely,the IG (r ≤ . ◦ ), SH (0 . ◦ > r ≤ . ◦ ), and FOH (r > . ◦ ). Here we fit just to the SH component. The IGand FOH components are difficult to disentangle fromstandard astrophysical processes and are not consideredin this study.The greatest uncertainty in our analysis is the deter-mination of the J -factor, which we calculate using theCLUMPY code. This uncertainty arises from two mainfactors. First, there is a high uncertainty in the sub-structure nature of the DM halo’s for both M31 and theMW, as well as an uncertainty in the halo geometries. Tobracket the substructure uncertainty we vary the subhalomass function, the fraction of the halo resolved in sub-structure, and the minimum subhalo mass in the ranges1 . − . 0, 0 . − . 35, and 10 − − M (cid:12) , respectively. Forthe concentration-mass relation we adopt the model from1Ref. [9]. The largest subhalo mass is taken to be 10% themass of the host halo. The calculations include 2 levels ofsubstructure. For the underlying smooth density profileswe test both an NFW profile and an Einasto profile. Thespatial distribution of subhaloes and the density profileof the subhaloes are assumed to be the same as the den-sity profile of the main halo. All calculations are madeself-consistently for M31 and the MW (i.e. they have thesame halo paramaters). Our calculated total boost fac-tor ranges from ∼ − J -factor for the M31 field. To do this we haveused the range of different halo shapes found in the lit-erature. For the MW we find that the halo shape maychange the J -factor in the range J/J Sph = 0 . − . J -factor for theM31 field is the contribution from the MW’s DM haloalong the line of sight. In Ref. [32] a detailed modelingof the foreground emission was performed, as well as anin-depth analysis of the corresponding systematic uncer-tainties. However, the model does not explicitly accountfor a potential contribution from the MW’s extended DMhalo. It is likely that such a signal could be (partially)absorbed by the isotropic component. The magnitudeof this effect, however, depends on the specific halo ge-ometry and substructure properties of the MW DM haloin the M31 field, which are not well constrained. In or-der to help control this, Ref. [32] used a region belowthe M31 field to tune the isotropic normalization. Here,we improve on this determination by considering varia-tions of the MW DM component in the M31 field and inthe tuning region due to different halo geometries. Wefind that the ratio is significant and, more specifically,in the range of J MW /J TR = 0 . − . 4. Thus even inthe ideal case where the isotropic component is able toperfectly absorb the emission from the MW’s DM halo,there could still be a gradient in the M31 field that isnot included in the foreground model and is likely to bea significant component in this region. Since the uncer-tainty in the J -factor due to the contribution from theMW’s DM halo along the line of sight is significant butcannot be precisely constrained, here we consider the twoextreme cases: one where none of the MW halo compo-nent has been absorbed by the isotropic component, andso J total = J M31 + J MW (case I); the other where theMW component has been completely absorbed so that J total = J M31 (case II).When these uncertainties are taken into account, wefind that the observed excess in the outer halo of M31favors a DM particle with a mass of ∼ ∼ × − − × − cm s − . We compare the best-fit DM param-eters for M31’s outer halo to numerous complementarytargets. We conclude that for the DM interpretation ofthe M31 outer halo excess to be compatible with the GCexcess, anti-proton excess, and current indirect detectionconstraints, it requires the J -factor to be towards thehigher end of the uncertainty range. This in turn hastwo main implications. First, the minimum subhalo massmust be < ∼ − M (cid:12) . And in fact this is expected inthe standard DM paradigm (ΛCDM). Second, the signalmust have a significant contribution from the MW’s DMhalo along the line of sight, i.e. it is too bright to be orig-inating from M31 alone. This condition cannot be ruledout, and it is in fact likely that some fraction of the MWDM halo emission is embedded in the signal toward M31.This is a feature of the methodology employed to tunethe MW foreground, as discussed in this paper. Giventhese conditions hold, we find that there is a large over-lap with the DM interpretations of both the GC excessand the antiproton excess, while also being compatiblewith the limits from the MW dwarfs. Although the un-certainty in the current measurements is clearly far toolarge to make any robust conclusions (either positive ornegative), this region in parameter space still remainsviable for discovery of the DM particle.Future prospects to confirm the excess toward theouter halo of M31, and to better understand its nature,crucially rely on improvements in modeling the interstel-lar emission towards M31. Furthermore, observations ofthe halos of other galaxies, e.g. M33, could provide a con-firmation of this type of signal, provided sufficient data isavailable since the signal is predicted to be fainter there.Other prospects may include a study of the distributionof properties of the isotropic background around the di-rection to M31 and further out with a goal to see the dis-tortions in the MW DM halo. Alternatively, constraintson the subhalo population by other astrophysical probesand, in turn, on their contribution to the M31 signal,might also provide a further test of the viability of theDM interpretation. ACKNOWLEDGEMENTS C.K. is pleased to acknowledge conversations withDaniel McKeown and James Bullock. I.M. acknowledgessupport from NASA grant No. NNX17AB48G. [1] H. W. Babcock, Lick Observatory Bulletin , 41(1939). [2] V. C. Rubin and W. K. Ford, Jr., Astrophys. J. ,379 (1970). [3] M. S. Roberts and R. N. Whitehurst, Astrophys. J. ,327 (1975).[4] C. Carignan, L. Chemin, W. K. Huchtmeier, and F. J.Lockman, Astrophys. J. Lett. , L109 (2006).[5] E. Corbelli, S. Lorenzoni, R. Walterbos, R. Braun, andD. Thilker, Astron. Astrophys. , A89 (2010).[6] M. Kamionkowski and A. Kinkhabwala, Phys. Rev. D , 3256 (1998).[7] R. Braun and W. B. Burton, Astron. Astrophys. ,437 (1999).[8] L. Blitz, D. N. Spergel, P. J. Teuben, D. Hartmann, andW. B. Burton, Astrophys. J. , 818 (1999).[9] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville,A. V. Kravtsov, A. A. Klypin, J. R. Primack, andA. Dekel, Mon. Not. Roy. Astron. Soc. , 559 (2001).[10] V. de Heij, R. Braun, and W. B. Burton, Astron. As-trophys. , 417 (2002).[11] R. Braun and D. Thilker, Astron. Astrophys. , 421(2004).[12] A. Helmi, Mon. Not. Roy. Astron. Soc. , 643 (2004).[13] J. Bailin and M. Steinmetz, Astrophys. J. , 647(2005).[14] B. Allgood, R. A. Flores, J. R. Primack, A. V. Kravtsov,R. H. Wechsler, A. Faltenbacher, and J. S. Bullock,Mon. Not. Roy. Astron. Soc. , 1781 (2006).[15] P. Bett, V. Eke, C. S. Frenk, A. Jenkins, J. Helly,and J. Navarro, Mon. Not. Roy. Astron. Soc. , 215(2007).[16] E. Hayashi, J. F. Navarro, and V. Springel, Mon. Not.Roy. Astron. Soc. , 50 (2007).[17] M. Kuhlen, J. Diemand, and P. Madau, Astrophys. J. , 1135 (2007).[18] A. Banerjee and C. J. Jog, Astrophys. J. , 254(2008).[19] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow,A. Jenkins, A. Helmi, J. F. Navarro, C. S. Frenk, andS. D. M. White, Mon. Not. Roy. Astron. Soc. , 1685(2008).[20] J. Diemand, M. Kuhlen, P. Madau, M. Zemp, B. Moore,D. Potter, and J. Stadel, Nature , 735 (2008).[21] M. Zemp, J. Diemand, M. Kuhlen, P. Madau, B. Moore,D. Potter, J. Stadel, and L. Widrow, Mon. Not. Roy.Astron. Soc. , 641 (2009).[22] K. Saha, E. S. Levine, C. J. Jog, and L. Blitz, Astrophys.J. , 2015 (2009).[23] D. R. Law, S. R. Majewski, and K. V. Johnston, Astro-phys. J. , L67 (2009).[24] A. Banerjee and C. J. Jog, Astrophys. J. , L8 (2011).[25] L. Gao, C. S. Frenk, A. Jenkins, V. Springel, andS. D. M. White, Mon. Not. Roy. Astron. Soc. , 1721(2012).[26] T. Ishiyama, Astrophys. J. , 27 (2014).[27] S. Garrison-Kimmel, M. Boylan-Kolchin, J. S. Bullock,and K. Lee, Mon. Not. Roy. Astron. Soc. , 2578(2014).[28] M. Velliscig, M. Cacciato, J. Schaye, R. A. Crain, R. G.Bower, M. P. van Daalen, C. D. Vecchia, C. S. Frenk,M. Furlong, I. G. McCarthy, et al. , Mon. Not. Roy. As-tron. Soc. , 721 (2015).[29] N. Bernal, L. Necib, and T. R. Slatyer, J. Cosmol. As-tropart. Phys. (12), 030.[30] ´A. Molin´e, M. A. S´anchez-Conde, S. Palomares-Ruiz,and F. Prada, Mon. Not. Roy. Astron. Soc. , 4974 (2017).[31] S. Garrison-Kimmel, A. Wetzel, J. S. Bullock, P. F.Hopkins, M. Boylan-Kolchin, C.-A. Faucher-Gigu`ere,D. Kereˇs, E. Quataert, R. E. Sanderson, A. S. Graus, et al. , Mon. Not. Roy. Astron. Soc. , 1709 (2017).[32] C. M. Karwin, S. Murgia, S. Campbell, and I. V.Moskalenko, Astrophys. J. , 95 (2019).[33] A. Falvard, E. Giraud, A. Jacholkowska, J. Lavalle,E. Nuss, F. Piron, M. Sapinski, P. Salati, R. Taillet,K. Jedamzik, and G. Moultaka, Astropart. Phys. ,467 (2004).[34] N. Fornengo, L. Pieri, and S. Scopel, Phys. Rev. D ,103529 (2004).[35] G. D. Mack, T. D. Jacques, J. F. Beacom, N. F. Bell,and H. Yuksel, Phys. Rev. D , 063542 (2008).[36] L. Dugger, T. E. Jeltema, and S. Profumo, J. Cosmol.Astropart. Phys. (12), 015.[37] J. Conrad, J. Cohen-Tanugi, and L. E. Strigari, J. Exp.Theor. Phys. , 1104 (2015), [Zh. Eksp. Teor. Fiz.148, no.6, 1257 (2015)].[38] J. M. Gaskins, Con. Phys. , 496 (2016).[39] M. Lisanti, S. Mishra-Sharma, N. L. Rodd, and B. R.Safdi, Phys. Rev. Lett. , 101101 (2018).[40] M. Lisanti, S. Mishra-Sharma, N. L. Rodd, B. R. Safdi,and R. H. Wechsler, Phys. Rev. D , 063005 (2018).[41] A. Abdo et al. (Fermi-LAT), Astron. Astrophys. ,L2 (2010).[42] M. S. Pshirkov, V. V. Vasiliev, and K. A. Postnov, Mon.Not. Roy. Astron. Soc. , L76 (2016).[43] M. Ackermann et al. (Fermi-LAT), Astrophys. J. ,208 (2017).[44] C. Eckner, X. Hou, P. D. Serpico, M. Winter, G. Zahar-ijas, P. Martin, M. di Mauro, N. Mirabal, J. Petrovic,T. Prodanovic, et al. , Astrophys. J. , 79 (2018).[45] A. McDaniel, T. Jeltema, and S. Profumo, Phys. Rev.D , 103021 (2018).[46] M. Di Mauro, X. Hou, C. Eckner, G. Zaharijas, andE. Charles, Phys. Rev. D99 , 123027 (2019).[47] A. McDaniel, T. Jeltema, and S. Profumo, Phys. Rev. D100 , 023014 (2019).[48] M. Ajello, M. D. Mauro, V. S. Paliya, and S. Garrappa,Astrophys. J. , 88 (2020).[49] P. Kroupa, C. Theis, and C. M. Boily, Astron. Astro-phys. , 517 (2005).[50] M. S. Pawlowski, J. Pflamm-Altenburg, and P. Kroupa,Mon. Not. Roy. Astron. Soc. , 1109 (2012).[51] A. R. Conn, G. F. Lewis, R. A. Ibata, Q. A. Parker,D. B. Zucker, A. W. McConnachie, N. F. Martin,D. Valls-Gabaud, N. Tanvir, M. J. Irwin, A. M. N.Ferguson, and S. C. Chapman, Astrophys. J. , 120(2013).[52] R. A. Ibata, G. F. Lewis, A. R. Conn, M. J. Irwin,A. W. McConnachie, S. C. Chapman, M. L. Collins,M. Fardal, A. M. Ferguson, N. G. Ibata, et al. , Nature , 62 (2013).[53] M. S. Pawlowski, P. Kroupa, and H. Jerjen, Mon. Not.Roy. Astron. Soc. , 1928 (2013).[54] F. Hammer, Y. Yang, S. Fouquet, M. S. Pawlowski,P. Kroupa, M. Puech, H. Flores, and J. Wang, Mon.Not. Roy. Astron. Soc. , 3543 (2013).[55] M. S. Pawlowski, Mod. Phys. Lett. A A33 , 1830004(2018).[56] M. S. Pawlowski, R. A. Ibata, and J. S. Bullock, Astro-phys. J. , 132 (2017). [57] A. N. Hulsbosch, Astron. Astrophys. , 1 (1975).[58] F. J. Lockman, Astrophys. J. Lett. , L33 (2003).[59] J. D. Simon, L. Blitz, A. A. Cole, M. D. Weinberg, andM. Cohen, Astrophys. J. , 270 (2006).[60] J. Kerp, P. Kalberla, N. B. Bekhti, L. Fl¨oer, D. Lenz,and B. Winkel, Astron. Astrophys. , A120 (2016).[61] N. B. Bekhti, L. Fl¨oer, R. Keller, J. Kerp, D. Lenz,B. Winkel, J. Bailin, M. Calabretta, L. Dedes, H. Ford, et al. , Astron. Astrophys. , A116 (2016).[62] M. Ajello, M. D. Mauro, V. S. Paliya, and S. Garrappa,Astrophys. J. , 88 (2020).[63] S.-Q. Xi, H.-M. Zhang, R.-Y. Liu, and X.-Y. Wang,arXiv e-prints , arXiv:2003.07830 (2020).[64] M. Wright, Astrophys. J. , 35 (1979).[65] O. Keenan, J. I. Davies, R. Taylor, and R. Minchin,Mon. Not. Roy. Astron. Soc. , 951 (2015).[66] D. A. Thilker, R. Braun, and R. Walterbos, in SeeingThrough the Dust: The Detection of HI and the Explo-ration of the ISM in Galaxies , Vol. 276 (2002) p. 370.[67] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi,M. Kadastik, P. Panci, M. Raidal, F. Sala, and A. Stru-mia, J. Cosmol. Astropart. Phys. (3), 051, [Erra-tum: J. Cosmol. Astropart. Phys. 1210, E01 (2012)].[68] P. Ciafaloni, D. Comelli, A. Riotto, F. Sala, A. Strumia,and A. Urbano, J. Cosmol. Astropart. Phys. (3),019.[69] J. Lyu, G. H. Rieke, and S. Alberts, Astrophys. J. ,85 (2016).[70] T. Isobe, E. D. Feigelson, and P. I. Nelson, Astrophys.J. , 490 (1986).[71] A. Charbonnier, C. Combet, and D. Maurin, Comput.Phys. Commun. , 656 (2012).[72] V. Bonnivard, M. H¨utten, E. Nezri, A. Charbonnier,C. Combet, and D. Maurin, Comput. Phys. Commun. , 336 (2016).[73] M. H¨utten, C. Combet, and D. Maurin, Comput. Phys.Commun. , 336 (2019).[74] A. M. Green, S. Hofmann, and D. J. Schwarz, Mon. Not.Roy. Astron. Soc. , L23 (2004).[75] T. Ishiyama and S. Ando, arXiv preprintarXiv:1907.03642 (2019).[76] M. A. S´anchez-Conde and F. Prada, Mon. Not. Roy.Astron. Soc. , 2271 (2014).[77] J. Wang, S. Bose, C. S. Frenk, L. Gao, A. Jenk-ins, V. Springel, and S. D. M. White, arXiv e-prints, arXiv:1911.09720 (2019).[78] J. Diemand, M. Kuhlen, and P. Madau, Astrophys. J. , 262 (2007).[79] A. Tamm, E. Tempel, P. Tenjes, O. Tihhonova, andT. Tuvikene, Astron. Astrophys. , A4 (2012).[80] C. Wegg, O. Gerhard, and M. Bieth, Mon. Not. R. As-tron. Soc. , 3296 (2019).[81] A. Bowden, N. W. Evans, and A. A. Williams, Mon.Not. R. Astron. Soc. , 329–337 (2016).[82] L. Posti and A. Helmi, Astron. Astrophys. , A56(2019).[83] S. R. Loebman, ˇZ. Ivezi´c, T. R. Quinn, J. Bovy, C. R.Christensen, M. Juri´c, R. Roˇskar, A. M. Brooks, andF. Governato, Astrophys. J. , 151 (2014).[84] G. Iorio and V. Belokurov, Mon. Not. R. Astron. Soc. , 3868 (2019).[85] N. W. Evans, C. A. O’Hare, and C. McCabe, Phys. Rev.D , 023012 (2019). [86] R. Emami, S. Genel, L. Hernquist, C. Alcock,S. Bose, R. Weinberger, M. Vogelsberger, F. Marinacci,A. Loeb, P. Torrey, and J. C. Forbes, arXiv e-prints ,arXiv:2009.09220 (2020).[87] E. Vasiliev, V. Belokurov, and D. Erkal, arXiv e-prints, arXiv:2009.10726 (2020).[88] K. N. Abazajian and R. E. Keeley, Phys. Rev. D ,083514 (2016).[89] A. Cuoco, M. Kr¨amer, and M. Korsmeier, Phys. Rev.Lett. , 191102 (2017).[90] I. Cholis, T. Linden, and D. Hooper, Phys. Rev. D ,103026 (2019).[91] A. Reinert and M. W. Winkler, J. Cosmol. Astropart.Phys. (01), 055.[92] M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. , 231301 (2015).[93] A. Albert et al. (Fermi-LAT, DES), Astrophys. J. ,110 (2017).[94] S. Ando, A. Geringer-Sameth, N. Hiroshima, S. Hoof,R. Trotta, and M. G. Walker, arXiv e-prints ,arXiv:2002.11956 (2020).[95] V. Bonnivard, C. Combet, D. Maurin, and M. Walker,Mon. Not. Roy. Astron. Soc. , 3002 (2015).[96] N. Klop, F. Zandanel, K. Hayashi, and S. Ando, Phys.Rev. D , 123012 (2017).[97] C. Karwin, S. Murgia, T. M. P. Tait, T. A. Porter, andP. Tanedo, Phys. Rev. D , 103005 (2017).[98] K. N. Abazajian, N. Canac, S. Horiuchi, andM. Kaplinghat, Phys. Rev. D , 023526 (2014).[99] F. Calore, I. Cholis, and C. Weniger, J. Cosmol. As-tropart. Phys. (03), 038.[100] T. Daylan, D. P. Finkbeiner, D. Hooper, T. Linden,S. K. N. Portillo, N. L. Rodd, and T. R. Slatyer, Phys.Dark Universe , 1 (2016).[101] C. Gordon and O. Macias, Phys. Rev. D , 083521(2013), [Erratum: Phys. Rev. D 89, no.4, 049901(2014)].[102] M. R. Buckley, E. Charles, J. M. Gaskins, A. M. Brooks,A. Drlica-Wagner, P. Martin, and G. Zhao, Phys. Rev. D91 , 102001 (2015).[103] R. Caputo, M. R. Buckley, P. Martin, E. Charles,A. M. Brooks, A. Drlica-Wagner, J. M. Gaskins, andM. Wood, Phys. Rev. D93 , 062004 (2016).[104] M. Ajello et al. , Astrophys. J. , L27 (2015).[105] M. Ackermann et al. (Fermi-LAT), Astrophys. J. ,91 (2012).[106] D. Hooper and S. J. Witte, J. Cosmol. Astropart. Phys. (04), 018.[107] I. Cholis, D. Hooper, and T. Linden, Phys. Rev. D ,083507 (2015).[108] A. E. Egorov and E. Pierpaoli, Phys. Rev. D88 , 023504(2013). Appendix A: DM Parameter Space Here we summarize all of the results overlaid in Fig-ure 7. The black data points (furthest four to the right)are for a DM interpretation of the GC excess, as pre-sented in Ref. [97]. The two points at lower energy arefor two of the models employed for the fore/background γ -ray emission from the MW, OB stars index-scaled , and4the points at higher energy are for the other two models, pulsars index-scaled . The NFW profile has γ = 1 . γ = 1 . s = 20 kpc and ρ (cid:12) = 0 . − . Note that theannihilation final state preferred by the fit to the datafavors mostly bottom-type quarks in the final state, witha small fraction of leptonic final states. Thus this modelis not directly comparable to the other overlays whichgenerally assume annihilation into a single final state.The black contour that is highly elongated in the y-direction is for the GC excess from Ref. [88]. The con-tour represents the total uncertainty (3 σ statistical +systematic). The uncertainty is dominated by the sys-tematics, and in particular, the value of the local DMdensity (this study also considers uncertainties due tothe index and scale radius of the DM profile, γ and R s ). The upper region of the contour corresponds to ρ (cid:12) = 0 . 28 GeV cm − (which is taken as the benchmarkvalue), and the lower region of the contour correspondsto ρ (cid:12) = 0 . 49 GeV cm − . The shift occurs at a crosssection value of ∼ × − cm s − . See Ref. [88] fordetails. Also plotted in Figure 7 is the best-fit point fromRef. [98] (the black data point to the far left).Other contours for the GC excess are also shown withdifferent shades of grey. The lowest and darkest contour(2 σ ) is from Ref. [99], then above that is the contour(2 σ ) from Ref. [100], and above that is the contour fromRef. [101]. The NFW profiles for all of these contourshave γ = 1 . R s = 20 kpc, and ρ (cid:12) = 0 . − .The two lowest purple curves show limits for the MWsatellite galaxies. The dashed curve is from Ref. [92]and results from the combined analysis of 15 dwarfspheroidal galaxies using Pass-8 data. The solid curveis from Ref. [93] and results from the combined analy-sis of 45 stellar systems, including 28 kinematically con-firmed dark-matter-dominated dwarf spheroidal galaxies,and 17 recently discovered systems that are dwarf can-didates. Note that the dwarf limits are obtained by as-suming spherical symmetry of the DM halos; however, ifthe halos are non-spherical then the limits may be weak-ened, as discussed in Refs. [95, 96]. We also plot the lim-its from Ref. [94] ( V = 10 . − ), which employssemi-analytic models of DM subhalos to derive realisticsatellite priors on the J -factor (for the ultrafaint dwarfs).This result explicitly exemplifies the uncertainty rangeassociated with limits from the MW dwarfs.The two highest purple curves are for the LMC andSMC. The dash-dot curve shows 2 σ limits from the LMCfrom Ref. [102], based on Pass-7 data. The dotted curveshows 2 σ limits from the SMC from Ref. [103].The tan band shows the 2 σ upper-limit from the extra-galactic γ -ray background (EGB) from Ref. [104]. Theband reflects the uncertainties related to the modeling ofDM subhaloes. This analysis shows that blazars, star-forming galaxies, and radio galaxies can naturally ac-count for the amplitude and spectral shape of the EGBover the energy range 0.1–820 GeV, leaving only modestroom for other contributions. The blue curve shows γ -ray limits (3 σ ) from the MWhalo from Ref. [105]. This is the limit obtained withmodeling the MW diffuse emission using GALPROP, foran NFW profile, with γ = 1 and a local DM density of0.43 GeV cm − . The limits are generally weaker withoutmodeling the diffuse emission, and they have a strongdependence on the local DM density.The light purple curve is for DM subhaloes fromRef. [106]. These limits are based on DM subhalo candi-dates from the unassociated point sources detected by Fermi -LAT. In total there are 19 subhalo candidates.The minimum subhalo mass for the upper limit calcu-lation is assumed to be 10 − M (cid:12) .The upper gray band in Figure 7 shows radio con-straints for the GC from Ref. [107]. The limits are derivedusing VLA observations at 330 MHz of the central 0 . ◦ around Sgr A*. An NFW profile is used with γ = 1 . s = 20 kpc, a local DM density of 0.3 GeV cm − , anda flat density core of 2 pc. The limits include energylosses due to IC and convection. The lower limit is for V C = 0 km s − , and the upper limit (not shown) is for V C = 1000 km s − . The limits can be much stronger (upto 3 or 4 orders of magnitude) when not including IC andconvection, or for a core radius closer to zero. There isalso a high uncertainty of the magnetic field strength inthe innermost region of the GC.The lower gray band shows radio limits from the cen-tral region of M31 ( ∼ B = 5 µ G and DM concentra-tion of c = 12, the middle region is for B = 50 µ G andDM concentration of c = 20, and the lowest region isfor B = 300 µ G and DM concentration of c = 28. AnNFW profile is used for the DM density, with γ = 1, and aflat core for r < 50 pc. The limits have a large uncertaintydue to the uncertainties in the DM profile and magneticfield strength in the inner regions of M31. The mag-netic field is modeled with an exponential dependence ingalactocentric radius and height above the galactic plane.The analysis accounts for leptonic energy losses due toIC emission, synchrotron emission, Bremsstrahlung, andCoulomb scattering, with synchrotron emission being thedominant loss mechanism over most of the energy range.We note, however, that uncertainties in the astrophysicalmodeling of these processes may weaken the limits evenfurther. In particular, the limits have a strong depen-dence on the relative strength of the inverse Comptonlosses compared to the synchrotron losses, which in turndepends on the energy density of M31’s interstellar radi-ation field.Also shown are contours for a recently reported ex-cesses in the flux of antiprotons. The upper light teal con-tour (2 σ ) is from Ref. [89]. The lower dark contour (2 σ )is from Ref. [90]. The NFW profiles for these contourshave γ = 1 . R s = 20 kpc, and ρ (cid:12) = 0 . − .5The teal curve shows upper-limits from Ref. [91], wherea less optimistic view of the excess is given (although thelimits still clearly show an anomaly around the signalregion).The red curve is for M31’s inner galaxy from Ref. [46].These limits are obtained by assuming that all of the ob-served γ -ray emission from M31’s inner galaxy arises from standard astrophysical emission, and therefore includinga 0 . ◦ disk template (which is derived directly from thebright γ -ray emission that is observed) in the DM fit.In addition, to account for the foreground/backgroundemission, the standard IEM is fit directly to the γγ