The angular power spectrum of the diffuse gamma-ray emission as measured by the Fermi Large Area Telescope and constraints on its Dark Matter interpretation
Mattia Fornasa, Alessandro Cuoco, Jesus Zavala, Jennifer M. Gaskins, Miguel A. Sanchez-Conde, German Gomez-Vargas, Eiichiro Komatsu, Tim Linden, Francisco Prada, Fabio Zandanel, Aldo Morselli
aa r X i v : . [ a s t r o - ph . H E ] D ec The angular power spectrum of the diffuse gamma-ray emission as measured by the
Fermi
Large Area Telescope and constraints on its Dark Matter interpretation
Mattia Fornasa, ∗ Alessandro Cuoco, † Jes´us Zavala,
3, 4, ‡ Jennifer M. Gaskins,
1, 5
Miguel A. S´anchez-Conde,
6, 7
German Gomez-Vargas, Eiichiro Komatsu,
9, 10
Tim Linden,
11, 12
Francisco Prada,
13, 14, 15
Fabio Zandanel, and Aldo Morselli GRAPPA, University of Amsterdam, Science Park, 1098 XH Amsterdam, Netherlands Institute for Theoretical Particle Physics and Cosmology (TTK),RWTH Aachen University, D-52056 Aachen, Germany Dark Cosmology Centre, Niels Bohr Institute, University of Copenhagen,Juliane Maries Vej 30, 2100 Copenhagen, Denmark Center for Astrophysics and Cosmology, Science Institute,University of Iceland, Dunhagi 5, 107 Reykjavik, Iceland California Institute of Technology, Pasadena, California, 91125, United States of America The Oskar Klein Centre for Cosmoparticle Physics, AlbaNova, SE-106 91 Stockholm, Sweden Department of Physics, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden Instituto de Astrofis´ıca, Pontificia Universidad Catolica de Chile,Avenida Vicuna Mackenna 4860, Santiago, Chile Max-Planck-Institut f¨ur Astrophysik, 85740 Garching bei M¨unchen, Germany Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI),Todai Institutes for Advanced Study, The University of Tokyo, Kashiwa 277-8583, Japan University of Chicago, Kavli Institute for Cosmological Physics,Chicago, Illinois, 60637, United States of America Ohio State University, Center for Cosmology and AstroParticle Physcis (CCAPP),Columbus, Ohio, 43210, United States of America Instituto de F´ısica Te´orica, (UAM/CSIC), Universidad Aut´onoma de Madrid, Cantoblanco, E-28049 Madrid, Spain Campus of International Excellence UAM+CSIC, Cantoblanco, E-28049 Madrid, Spain Instituto de Astrof´ısica de Andaluc´ıa (IAA-CSIC),Glorieta de la Astronom´ıa, E-18008, Granada, Spain Istituto Nazionale di Fisica Nucleare, Sezione di Roma ”Tor Vergata”, I-00133 Roma, Italy (Dated: December 22, 2016)The isotropic gamma-ray background arises from the contribution of unresolved sources, includingmembers of confirmed source classes and proposed gamma-ray emitters such as the radiation inducedby dark matter annihilation and decay. Clues about the properties of the contributing sources areimprinted in the anisotropy characteristics of the gamma-ray background. We use 81 months of Pass7 Reprocessed data from the
Fermi
Large Area Telescope to perform a measurement of the anisotropyangular power spectrum of the gamma-ray background. We analyze energies between 0.5 and 500GeV, extending the range considered in the previous measurement based on 22 months of data.We also compute, for the first time, the cross-correlation angular power spectrum between differentenergy bins. We find that the derived angular spectra are compatible with being Poissonian, i.e.constant in multipole. Moreover, the energy dependence of the anisotropy suggests that the signalis due to two populations of sources, contributing, respectively, below and above ∼ I. INTRODUCTION
In 2012, the
Fermi
Large Area Telescope (LAT)Collaboration measured for the first time the auto-correlation angular power spectrum (auto-APS) of thediffuse gamma-ray emission detected far from the ∗ [email protected] † [email protected] ‡ [email protected]; Marie Curie Fellow Galactic plane [1]. In that analysis, point sources inthe first
Fermi
LAT source catalog (1FGL) [2] and aband along the Galactic plane with Galactic latitude | b | < ◦ were masked in order to isolate the contributionto the auto-APS from the so-called Isotropic Gamma-RayBackground (IGRB).The IGRB is what remains of the gamma-ray sky afterthe subtraction of the emission from resolved sources andfrom the Galactic diffuse foreground induced by cosmicrays [3, 4]. It dominates the gamma-ray sky at largeGalactic latitudes and its intensity energy spectrum isfound to be compatible with a power-law with a slope of2 . ± .
02 between 100 MeV and ∼
300 GeV, and with anexponential cut-off at higher energies [4]. These valuesfor the spectral slope and for the energy cut-off are thosefound when “model A” from Ref. [4] is used to describethe Galactic diffuse foreground emission. A differentforeground model for the Galaxy would lead to a slightlydifferent energy spectrum for the IGRB. Deviations canbe as large as 20-30% depending on energy.The IGRB is interpreted as the cumulative emission ofsources (e.g., blazars, star-forming and radio galaxies)that are too faint to be detected individually (seeRef. [5] for a recent review and the references therein) .Yet, its exact composition remains unknown. It isexpected to be isotropic on large angular scales but itcan still contain anisotropies on small angular scales.Indeed, the contribution to the IGRB from unresolvedsources imprints anisotropies in the diffuse emissionwhich can be used to infer the properties of thecontributing sources (see Refs. [6–15] among others).For example, the detection of a significant angularpower in Ref. [1] determined an upper limit to thecontribution of unresolved blazars [10, 11, 16] to theIGRB. Additional tools to reconstruct the nature ofthe IGRB are the study of its cross-correlation withcatalogs of resolved galaxies [17–22], with gravitationallensing cosmic shear [23] and with lensing of the cosmicmicrowave background radiation [24]. Complementaryinformation can also be inferred by modeling its 1-pointphoton count distribution [25–27].The detection of the auto-APS presented in Ref. [1]was based on ∼
22 months of data. Since then,
Fermi
LAT has increased its statistics by approximately afactor of 4. Therefore, we expect that an updatedmeasurement of the auto-APS will significantly improveour understanding of the IGRB. In the first part ofthis work, we perform this measurement by analyzing81 months of
Fermi
LAT data from 0.5 to 500 GeV,extending the 1-50 GeV energy range considered inRef. [1]. This enables a more precise characterization ofthe energy dependence of the auto-APS. Indeed, lookingfor features in the so-called “anisotropy energy spectrum”is a powerful way to single out different components ofthe IGRB [28]. We also compute, for the first time, thecross-correlation angular power spectrum (cross-APS) ofthe diffuse gamma-ray emission between different energybins. The cross-APS additionally enhances our ability tobreak down the IGRB into its different components sinceit provides information about the degree of correlationof the emission at different energies, which is stronger ifthe emission originates from one single source population(see, e.g., Ref. [29, 30]).In the second part of this paper, we focus ouranalysis on one possible contributor to the IGRB, namely Instead, the sum of the emission from the resolved and unresolvedsources is generally referred to as the Extragalactic Gamma-rayBackground. the emission induced by dark matter (DM). If DMis a weakly-interacting massive particle (WIMP), itsannihilation or decay could generate gamma rays. Theradiation produced in extragalactic and Galactic DMstructures could contribute to the IGRB (see Ref. [5]and references therein) and, therefore, the IGRB couldbe used to indirectly search for non-gravitational DMinteractions. Indeed, both the measurement of the IGRBenergy spectrum [4] and of its auto-APS [1] have beenalready used to set constraints on the possible DM-induced gamma-ray emission [14, 31–33].In this work, we also update the predictions for theauto- and cross-APS expected from DM annihilationor decay with respect to Ref. [34]. The distributionand properties of DM structures are modeled accordingto the results of state-of-the-art N -body cosmologicalsimulations. We also employ well-motivated semi-analytical recipes to account for the emission of DMstructures below the mass resolution of the simulations.The latter is a significant part of the expected signal, atleast in the case of annihilating DM. We take specialcare to estimate the uncertainties introduced whenmodeling the clustering of DM, especially at the smallestscales. Our predicted DM signal is then comparedto the updated Fermi
LAT measurement of the auto-and cross-APS. In the most conservative scenario, thiscomparison provides an upper limit to the gamma-rayproduction rate by DM particles, i.e. an upper limit toits annihilation cross section or a lower limit to its decaylifetime, as a function of DM mass.The paper is organized as follows: in Sec. II weprovide details on the data set that will be used inSec. III, where we describe our data analysis pipeline.We validate the latter in Sec. IV A on Monte Carlo (MC)simulations of the unresolved gamma-ray sky. In Sec. V,we present our results for the auto- and cross-APS, andwe describe the validation tests performed. Sec. VIprovides a phenomenological interpretation of our resultsin terms of one or multiple populations of gamma-raysources. In Sec. VII we focus on DM-induced gamma-ray emission: we provide details on how this signal issimulated, distinguishing among different componentsand discussing the main uncertainties affecting itscalculation. In Sec. VIII the auto- and cross-APSexpected from DM are compared to the measurementsand exclusion limits are derived. Finally, Sec. IXsummarizes our conclusions.
II. DATA SELECTION AND PROCESSING
The data analysis pipeline proceeds similarly to whatis described in Ref. [1]. We use Pass 7 Reprocessed
Fermi
LAT data taken between August 4 2008 and May 252015 (MET Range: 239557417 – 454279160), and restrictourselves to photons passing the ULTRACLEAN eventselection. Thus, we use P7REP ULTRACLEAN V15as the instrument response functions (IRFs). Weplace standard selection cuts on the
Fermi
LAT data,removing events entering the detector with a zenith angleexceeding 100 ◦ , events recorded when the Fermi
LATinstrument was oriented at a rocking angle exceeding 52 ◦ and events recorded while the Fermi
LAT was passingthrough the South Atlantic anomaly, or when it wasnot in science survey mode. Since photons which pair-convert in the front of the
Fermi
LAT detector havea better angular resolution, we split our data set intofront- and back-converting events, running each data setthrough the same data analysis pipeline. The front-converting events will represent our default deta set, withthe corresponding P7REP ULTRACLEAN V15 IRF. Toproduce flux maps we bin the resulting
Fermi
LAT eventcounts and exposure maps into
HEALPix -format maps [35] with angular bins of size ∼ ◦ ( HEALPix order10,
Nside =1024), as well as into 100 logarithmically-spaced energy bins spanning the energy range between104.46 MeV and 1044.65 GeV. The conversion of theexposure maps into
HEALPix -format maps is performedwith the GaRDiAn package [36]. Flux maps are, then,built by dividing the count map by the correspondingexposure map, in each energy bin. The flux mapsobtained with the fine energy binning are later co-addedinto 13 larger bins spanning the energy range between500 MeV and 500 GeV. This is done to ensure sufficientstatistics within each energy bin. We use the smallerenergy bins to calculate the beam window functionand the photon noise within each larger energy bin, asdescribed in Sec. III.
III. ANISOTROPY ANALYSISA. Auto- and cross-correlation angular powerspectra
An intensity sky map can be decomposed into sphericalharmonics as follows: I ( ψ ) = X ℓm a ℓ,m Y ℓ,m ( ψ ) , (1)where I ( ψ ) is the intensity from the line-of-sight direction ψ and Y ℓ,m ( ψ ) are the spherical harmonic functions. Theauto-APS C ℓ of the intensity map is given by the a ℓ,m coefficients as: C ℓ = 12 ℓ + 1 ℓ X m = − ℓ | a ℓ,m | . (2)Similarly, the cross-APS between two intensity maps I i and I j is constructed from the individual a iℓ,m and a jℓ,m http://healpix.jpl.nasa.gov coefficients, obtained from the decomposition in the twoenergy bins, independently: C ijℓ = 12 ℓ + 1 ℓ X m = − ℓ a iℓ,m a j⋆ℓ,m . (3)The auto- and cross-APS are computed with specificnumerical tools as, e.g., HEALPix and
PolSpice [37].However, before applying Eqs. 2 and 3, the data setmust be prepared, accounting for possible masking,foreground subtraction and pixelization. Additionallythe calculations are complicated by the finite angularresolution of the instrument. In the followingsubsections, we summarize how these aspects are takeninto consideration.
B. Masking
We apply a mask to the all-sky data to reducecontamination from Galactic diffuse foregrounds andfrom sources already detected in the third
Fermi
LATsource catalog (3FGL) [38]. The mask applied in ourdefault analysis excludes low Galactic latitudes ( | b | < ◦ ). We also mask each point-like source in 3FGL witha disk whose radius depends on the flux detected fromthe source between 0.1 and 100 GeV: for the 500 brightestsources we consider a disk with a radius of 3.5 ◦ , for thefollowing 500 sources a disk with a radius of 2.5 ◦ , a diskwith a radius of of 1.5 ◦ for the following 1000 sources,and, finally, a radius of 1.0 ◦ for the remaining objects.Validation of the choice for the mask will be performed inSec. V C. The 3FGL catalog contains 3 extended sourcesat moderate and high latitudes: Centaurus A and theLarge and Small Magellanic Clouds. Centaurus A andthe Large Magellanic Cloud are each masked excluding a10 ◦ -region from their center in the catalog. We employa 5 ◦ -mask for the Small Magellanic Cloud. The fraction f sky of the sky outside the mask is 0.275.We also consider an alternative mask that covers thesame strip around the Galactic plane but only the sourcesin the second Fermi
LAT source catalog (2FGL) [39]. Inthis case, we mask all the sources with a 2 ◦ -radius disk.The validation for this choice is performed in Sec. V Cand, in this case, f sky = 0 . C. Foreground cleaning
Despite applying a generous cut in Galactic latitude,some Galactic diffuse emission remains visible in theunmasked area of the sky map, particularly at lowenergies (see Fig. 1). To reduce this contaminationfurther, we perform foreground cleaning by subtracting
FIG. 1. Intensity maps (in cm − s − sr − ) in Galactic coordinates for energies between 1.0 and 2.0 GeV, shown unmasked ( top )and after applying the default mask removing sources in 3FGL, as described in Sec. III B ( bottom ). Data used here follow thedefault processing (see Sec. II), but they include both front- and back-converting events. Both maps have been smoothed witha gaussian beam with σ = 0 . ◦ and their projection scheme is Mollweide.FIG. 2. Same as the bottom panel of Fig. 1 but with our model for the Galactic foreground subtracted (see Sec. III C). Theresiduals have been smoothed with a gaussian beam with σ = 1 ◦ . The projection scheme is Mollweide. Multipole0 200 400 600 800 1000 ) l p i x , ( W ) l bea m ( W FIG. 3. Pixel and beam window functions as a functionof energy and for different data selections. The grayshort-dashed line shows the pixel window function of a
HEALPix map with
Nside =1024. The pixel windowfunction is independent of energy and IRFs. The solidand long-dashed lines show the beam window functions forP7REP ULTRACLEAN V15 IRFs. The solid lines are forfront-converting events and the long-dashed ones for back-converting events. The different colors stand for 4 differentrepresentative energies. a model of the Galactic diffuse emission. We usethe recommended model for Pass 7 Reprocessed dataanalysis, i.e. gll iem v05 rev1.fit . Details of thederivation of the model are described in Ref. [40]. Thisforeground model, together with an isotropic component,is fitted to the data in the unmasked region of the skyand in each one of the 13 coarser energy bins, using GaRDiAn . The default mask is adopted when fitting thediffuse components. The resulting best-fit model is thensubtracted from the intensity maps in each energy bin toobtain residual intensity maps, on which the anisotropymeasurements are performed. Fig. 2 shows an exampleof the residual intensity map for the data in the energybin between 1 and 2 GeV.We investigate the impact of foreground cleaning onthe auto- and cross-APS measurements in Sec. V C.
D. Noise and beam window functions
We calculate the auto- and cross-APS of the intensitymaps using the
PolSpice package [37] to deconvolve theeffect of the mask on the spectra and to provide thecovariance matrix for the estimated C ℓ .Both the finite angular resolution of the instrument(given by its point-spread function, PSF) and the finite http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundMod-els.html angular resolution of the map (i.e., the pixelizationscheme) suppress the measured auto- and cross-APS atlarge multipoles (i.e. small angular scales). This effectis described using the beam window function W beam ℓ andpixel window function W pix ℓ , respectively. We note thatthey affect the signal but not the noise term C N (seeRef. [1]). We use the beam and pixel window functionsto correct the suppression at large multipoles so that ourestimation for the auto- and cross-APS is as follows : C signal ,ijℓ = C Pol ,ijℓ − δ ij C i N ( W beam ,iℓ W beam ,jℓ )( W pix ℓ ) , (4)where the i and j indexes run from 1 to 13 and labelemission in different energy bins. The case i = j corresponds to the auto-APS and the one with i = j tothe cross-APS between energy bins. Also, C Pol ,ijℓ is theAPS delivered by
PolSpice , which is already correctedfor the effect of masking. The noise term δ ij C i N is equal tozero for the cross-APS since it is due to shot noise fromthe finite statistics of the gamma-ray events, which isuncorrelated between different energy bins. We compute C i N from the shot noise C k N of the 100 finely-gridded intensity maps, where C k N = h n kγ, pix / ( A k pix ) i Ω pix , (5)where n kγ, pix and A k pix are the number of observed eventsand the exposure, respectively, in each pixel and for the k -th finely-gridded energy bin. The averaging is doneover the unmasked pixels. Ω pix is the pixel solid angle,which is the same for each pixel. See Appendix A fora derivation of Eq. 5. The noise term C i N for the auto-APS in the i -th large energy bin is given by the sum ofthe noise terms in Eq. 5 of all the finely-gridded energybins covered by the i -th bin. We note that Eq. 5 ismore accurate than the shot noise used in Ref. [1], i.e. C N = h n γ, pix i / (Ω pix h A i ).The beam window function is computed as follows: W beam ℓ ( E ) = 2 π Z − d cos θP ℓ (cos( θ ))PSF( θ ; E) , (6)where P ℓ (cos( θ )) are the Legendre polynomials andPSF( θ ; E) is the energy-dependent PSF for a given set ofIRFs, with θ denoting the angular distance in the PSF.We use the gtpsf tool in the Science Tools packageto calculate the effective PSF, as a function of energy,averaged over the actual pointing and live-time historyof the LAT. The beam window functions are calculatedseparately for the P7REP ULTRACLEAN V15 front-and back-converting events. Finally, the pixel window In the remainder of the paper, we commonly refer to thisestimator simply by C ℓ instead of C signal ℓ . function W pix ℓ is computed using the tools provided inthe HEALPix package for
Nside ∼ ∼ ℓ > ∼ Fermi
LAT variessignificantly over the energy range considered in thisanalysis, and in some cases within the individual energybins used when computing the auto- and cross-APS, it isnecessary to calculate an effective beam window functionfor each energy bin. Therefore, for the i -th energy bin,we define the average window function h W beam ,iℓ i byweighting Eq. 6 with the intensity spectrum of the eventsin that bin outside the mask: h W beam ,iℓ i = 1 I bin Z E max ,i E min ,i d E W beam ℓ ( E ) d I ( E )d E , (7)where I bin ≡ R E max ,i E min ,i d E (d I/ d E ) and E min ,i and E max ,i are the lower and upper bounds of the i -th energy bin.We approximate the energy spectrum of the data byusing the measured differential intensity d I/ d E outsidethe mask in each intensity map for the finely-binnedenergy bins. IV. MONTE CARLO VALIDATION OF THEBINNING OF THE APS AND OF ITSPOISSONIAN FITA. Auto-correlation angular power spectrum
In this section we describe in detail the procedureused to bin the auto-APS estimated in Eq. 4 into largemultipole bins. Binning is required in order to reduce thecorrelation among nearby C ℓ due to the presence of themask.In contrast with the analysis of Ref. [1], in thepresent work the binned spectra C l are taken to be theunweighted average of the individual C ℓ in the bin. Also,the error σ ℓ on C l is computed by averaging all theentries of the covariance matrix provided by PolSpice inthe block corresponding to the bin under consideration.A dedicated set of MC simulations of all-sky data areproduced to validate these choices and to additionallytest alternative binning schemes. The MC validationprocedure is described below.
1. Monte Carlo simulations
The simulations are performed for a single energy binfrom 1 to 10 GeV. We assume an underlying populationof sources with a power-law source-count distribution,i.e., dN/dS = A ( S/S ) − α . The parameters A , S and α are fixed to the values 3 . × cm s sr − ,10 − cm − s − and 2.0, respectively, in agreement withthe best-fit results of Ref. [26]. We consider sourceswith fluxes (in the energy range between 1 and 10 GeV)from 10 − cm − s − to 10 − cm − s − . The upper valueis roughly equal to the 3FGL sensitivity threshold. Inthis way, the level of anisotropy expected from thesesources is roughly equal to that observed in the datawhen masking the 3FGL sources. The lower value is notcrucial since the auto-APS is dominated by the sourcesjust below the detection threshold. From the sourcecount distribution dN/dS , we create a realization of thesource population, producing about 40,000 objects andassigning them random positions on the sky. This createsa map with a Poissonian (i.e., constant in multipole)auto-APS, C P , whose value can be computed by summingtogether the squared flux, Φ i , of all the simulated sourcesdivided by 4 π : C P = P i Φ i / π . This is equivalent tothe usual way of calculating C P by integrating S dN/dS over the range in flux mentioned above. The resultingPoissonian auto-APS C P is 3 . × − cm − s − sr − .This is the nominal auto-APS that we want to recoverby applying our analysis pipeline to the simulations.We use the exposure (averaged in the energy rangebetween 1 and 10 GeV) for 5 years of data-taking toconvert the intensity map into a counts map. Themap is also convolved with the average PSF for theP7REP ULTRACLEAN V15 IRFs for front-convertingevents (averaged in the 1-10 GeV range, assuming anenergy spectrum ∝ E − . ). The result is a HEALPix -formatted map with resolution
Nside =1024 containingthe expected emission, in counts, from the simulatedsources. Purely isotropic emission is also included byadding an isotropic template to the map, which wasalso convolved with the IRFs and normalized to give thenumber of counts expected from the IGRB measured inthe 1-10 GeV energy range, including the contaminationfrom residual cosmic rays. For simplicity we did notmodel the Galactic foregrounds. This final map is thenPoisson-sampled pixel-by-pixel 200 times to yield 200different realizations of the expected counts. The auto-APS of each map is calculated with
PolSpice , afterapplying the default mask used in the analysis of thereal data, i.e., excluding the region with | b | < ◦ andthe sources in 3FGL, even though the simulation doesnot include those sources. Finally, noise subtractionand beam correction are also applied as described inSec. III D. sr] -2 sr -2 s -4 [cm l C0 1 2 3 4 5 6 7 8 -18 × N u m be r o f M C r ea li z a t i on s Unweighted average -2l σ = l Weighted average with weights w P Nominal CGaussian distribution with estimated by PolSpice σ FIG. 4. Comparison between different methods to bin theauto-APS measured in the bin between ℓ = 243 and 317from the MC simulations described in the text. The nominal C P is represented by the vertical grey line. The solid blackhistogram shows the distribution of the measured C ℓ forthe 200 simulated realizations, where the binned auto-APSis computed by an unweighted average. The dashed bluehistogram denotes the case of a weighted average with weightsgiven in Eq. 8. The solid red curves is a Gaussian distributioncentered on the nominal C P and with a standard deviation of9 . × − cm − s − sr − , as estimated with PolSpice .
2. Binning validation
We first validate our recipe to determine the binnedauto-APS. In this case, the standard analytic error σ ℓ oneach C ℓ (assuming that C ℓ follows a χ ℓ +1 distribution[41]) is: σ ℓ = s ℓ + 1) f sky (cid:18) C ℓ + C N W ℓ (cid:19) , (8)with W ℓ = W beam ℓ W pix ℓ . We test three approaches toobtain C ℓ : i ) computing the weighted average of the C ℓ in each multipole bin, using w ℓ = σ − ℓ as weight, ii )computing the weighted average of the C ℓ in the bin,with a weight w ℓ = σ − ℓ , defining σ ℓ as in Eq. 8 but only with the noise term C N /W ℓ and iii ) computing theunweighted average of the C ℓ in the bin. Note that in thefirst approach, the weight w ℓ depends on the data via the C ℓ term in Eq. 8, while in the second and third methodsthere is no dependence on the estimated auto-APS. Thefirst method is the one employed in Ref. [1].In Fig. 4, we show a histogram of the binned C ℓ in thebin between ℓ = 243 and 317 for the 200 MC realizations.The nominal C P is denoted by the grey vertical solidline. The solid black histogram refers to the case inwhich no weights are used (method iii ), while the dashedblue histogram is for the weighted average with weightsfrom Eq. 8 (method i ). The results for method ii (i.e., weighted average but with only the noise term in Eq. 8)are not plotted but they are similar to the solid blackhistogram. It is clear that binning the data by meansof a weighted average which includes the data C ℓ itselfgives a result which underestimates the nominal C P . Onthe other hand, using the unweighted average (as wedo in the current analysis) or weighting using only thenoise term gives results compatible with the input. Theintuitive reason for this bias can be traced to the fact thatmethod i ) uses the measured auto-APS in the estimationof the error: at each multipole the measured auto-APSfluctuates up and down significantly. If we use Eq. 8 withthe measured C ℓ to weight the data at each multipole, adownward fluctuation of C ℓ is assigned a smaller error barand, thus, a larger weight. This will lead to a downwardbiased C ℓ . Finally, the histograms also show that thedistribution of the C ℓ obtained from the MC realizationsis, to a good approximation, Gaussian. Indeed, itagrees well with the solid red curve representing aGaussian distribution centered on the nominal C P andwith a standard deviation of 9 . × − cm − s − sr − (seebelow).To assign an error σ ℓ to the binned auto-APS C ℓ we also test three methods: i ) the unweighted averageof σ ℓ from Eq. 8 in the bin, ii ) the weighted averageof σ ℓ from Eq. 8 with weight w ℓ = σ − ℓ and iii ) theaverage of the covariance matrix computed by PolSpice in the bin . Differently from the estimation of C ℓ , thethree methods for the estimation of σ ℓ produce similarresults. Thus, we decide to choose method iii ) as ourstandard prescription. This has also the advantage that,by averaging different blocks of the covariance matrixprovided by PolSpice , one can build a covariance matrixfor the binned auto-APS. The average of σ ℓ from method iii ) in the multipole bin between ℓ = 243 and 317over the 200 MC realizations is 9 . × − cm − s − sr − ,i.e. the value considered in Fig. 4 for the standarddeviation of the red curve. Our validation with MCsimulations shows that our estimate of the errors isrealiable and that higher-order effects, e.g. those relatedto the bispectrum and trispectrum discussed in Ref. [42],can be neglected. It remains interesting, nonetheless, tounderstand if a small bispectrum and trispectrum can beused to independently constrain the sources contributingto the IGRB.
3. Poissonian fit validation
Having validated the binning procedure for themeasured auto-APS, we are now interested in fitting thebinned auto-APS C ℓ with a constant value. Indeed, a PolSpice returns the covariance matrix of the beam-uncorrected C ℓ , denoted here by V ℓℓ ′ . In method iii ) the error σ ℓ is definedas P ℓℓ ′ V ℓℓ ′ / ( W ℓ W ℓ ′ ∆ ℓ ), where the sum runs over the ℓ, ℓ ′ inside each multipole bin and ∆ ℓ is the width of the bin. Poissonian APS C P (i.e. an APS that is constant inmultipole) is a natural expectation for the anisotropiesinduced by unclustered unresolved point sources. Onepossibility is to infer C P by minimizing the following χ function: χ ( C P ) = X ℓ ( C ℓ − C P ) σ ℓ , (9)where C ℓ and σ ℓ are the the binned data and their errors,as described in the previous section.A second possibility is to consider a likelihood function L that, up to a normalization constant, can be writtenas follows:log L ( C P ) = − X ℓ log( σ ℓ ) − X ℓ ( C ℓ − C P ) σ ℓ . (10)This expression for the likelihood takes into account thefact that σ ℓ also depends on C P , since σ ℓ in Eq. 10 isdefined as the average of σ ℓ = 2(2 ℓ + 1) f sky (cid:18) C P + C N W ℓ (cid:19) , (11)over the specific multipole bin. In fact, for largemultipoles, the expected χ ℓ +1 distribution of a given C ℓ can be approximated by a Gaussian for which themean and the standard deviation are not independent butrelated as in Eq. 11. Thus, the main difference betweenthe χ minimization (as in Eq. 9) and the likelihoodmethod is that, in the latter, σ ℓ depends on C P . Ignoringsuch a dependence may bias the result of the fit.The two methods described above are used todetermine the best-fit C P for the 200 MC realizationsdescribed above, by considering 10 C ℓ in 10 bins inmultipole uniformly spaced in log ℓ between ℓ = 49 and706. As we discuss in Sec. V, this multipole rangeexcludes the large angular scales where the reconstructed C ℓ are most uncertain due to possible contaminationof the Galactic foreground, and the high-multipolerange where the effect of the window functions becomestoo severe. The results are summarized in Fig. 5:the vertical grey line is the nominal C P , while thesolid black histogram shows the distribution of the C P determined by maximazing the log L of Eq. 10 if thebinned C ℓ are computed with no weights. This approachproduces a distribution that is approximately Gaussianand centered on the nominal C P . On the other hand,if the binned C ℓ are computed with the weights fromEq. 8, then the maximization of log L underestimatesthe Poissonian auto-APS (long-dashed blue histogramin Fig. 5). Making use of the χ function in Eq. 9instead of the log L in Eq. 10 gives similar results, i.e.an unbiased distribution for C P if the binned C ℓ arecomputed without weights (short-dashed pink line) andan underestimation of the nominal C P when weights are sr] -2 sr -2 s -4 [cm P C1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 -18 × N u m be r o f M C r ea li z a t i on s log(L)-fit with no weightslog(L)-fit with weights-fit with no weights χ P Nominal C estimated σ Gaussian distr. with from the log(L) fit with no weights
FIG. 5. Comparison between different methods to measurethe Poissonian C P in the MC simulations, given the binned C ℓ . The nominal C P is represented by the vertical greyline. The solid black histogram shows the distribution ofthe Poissonian C P for the 200 simulated realizations obtainedby maximazing the log L in Eq. 10 over the multipole rangefrom 49 to 706. The binned C ℓ in Eq. 10 are computedwith no weights. If the weighted average is considered,the distribution of C P is shown by the long-dashed bluehistogram, which is clearly biased low. The short-dashedpink histogram shows the distribution of C P computed bythe minimization of the χ in Eq. 9 from C ℓ binned withno weights. The solid red curve is a Gaussian distributioncentered on the nominal C P and with a standard deviation of2 . × − cm − s − sr − , as estimated from the log L method. included (not shown in Fig. 5) . The error associated tothe best-fit C P corresponds to the 68% confidence-level(CL) region. We note that the log L approach yieldsslightly smaller errors and we decide to adopt this as ourstandard way to measure the Poissonian auto-APS in thefollowing. The average of the error on the best-fit C P overthe 200 MC realizations is 2 . × − cm − s − sr − , i.e.the value used as the standard deviation for the Gaussianfunction plotted as the solid red line in Fig. 5, which iscentered on the nominal C P . B. Cross-correlation angular power spectrum
Similar checks to what is described above for the auto-APS are performed for the cross-APS between two energy Note that applying the log L or the χ approach to the unbinned C ℓ provided by PolSpice also leads to an underestimation of C P . bins. In this case, the standard analytical error is: σ ℓ = 1(2 ℓ + 1) f sky " C ℓ + C ,ℓ + C , N W ,ℓ ! C ,ℓ + C , N W ,ℓ ! , (12)where C ℓ is the cross-APS and C ,ℓ and C ,ℓ are the auto-APS for the two energy bins. Similarly, W ,ℓ and W ,ℓ are the window functions for the two energies consideredand C , N and C , N are the two photon noises. Aftertesting different averaging schemes, we decide to use thesame method as for the auto-APS case, i.e. to bin thecross-APS with an unweighted average and to estimate σ ℓ by computing the block-average of the covariance matrixprovided by PolSpice .Similarly, we tested the likelihood and χ approachto derive the Poissonian best-fit C P to the cross-APSdata. For the likelihood approach, σ ℓ is now defined asthe average of Eq. 12 after having replaced C ℓ by C P .As for the auto-APS, we find compatible results betweenthe two methods, with the likelihood approach providingslightly smaller errors. Therefore, in the following,we will quote Poissonian cross-APS derived with thismethod.We end this section by noting that the proper wayto estimate C P for the cross-APS would be to use thelikelihood method but replacing the auto-APS C ,ℓ and C ,ℓ in Eq. 12 with their Poissonian estimates C , P and C , P , and to perform a joint likelihood fit to all threequantities, i.e. C P , C , P and C , P . However, thisapproach would not provide results that are significantlydifferent than the ones obtained as described above. Infact, at present, the noise terms in Eq. 12 dominateover the signal terms, reducing the effect of covariancebetween energy bins . V. MEASURED AUTO- ANDCROSS-CORRELATION ANGULAR POWERSPECTRA OF THE ISOTROPIC GAMMA-RAYBACKGROUND
Following the analysis described in the previoussection, we measure the auto- and cross-APS in 13 energybins spanning the energy range between 500 MeV and 500GeV.
A. Auto-correlation angular power spectra
The auto-APS of the IGRB is shown in Fig. 6for two representative energy bins. The auto-APS for all 13 energy bins considered is shown The noise terms in Eq. 12 are a factor of 4-5 larger than C P , C , P and C , P . Therefore, not performing the joint likelihood fitas described in the text generates an error of, at most, 10-20%on σ ℓ . The effect on the estimated best-fit Poissonian auto- andcross-APS will be even smaller. y -axis range ofFig. 6 and in Appendix B has been chosen to betterillustrate the auto-APS in the multipole range of interest,i.e. between ℓ =49 and 706, divided into 10 bins equallyspaced in log ℓ . The red circles indicate the auto-APS forour reference data set (i.e., P7REP ULTRACLEAN V15front events) and the default mask covering the regionwith | b | < ◦ and 3FGL sources, as described inSec. III B. Instead, the blue triangles refer to the samedata set but using the default mask covering only 2FGLsources (see Sec. III B). Note that the blue triangles aresystematically higher than the red circles, due to theanisotropy power associated with the sources that arepresent in 3FGL but still unresolved in 2FGL.Fig. 7 shows the auto-APS for the two same energybins but over a broader multipole range, i.e. from ℓ = 10to 2000. This illustrates the behavior of the auto-APSabove and below the signal region used in our analysis,i.e. between ℓ = 49 and 706. At large scales (i.e., lowmultipoles), there might be some residual contaminationfrom the Galactic foregrounds. This motivates our choiceof neglecting the APS below ℓ = 49. In Sec. V C theeffect of foreground contamination is discussed in moredetail. On the other hand, at high multipoles and at lowenergies (left panel), the size of the error bars increasesdramatically due to the strong signal suppression causedby the beam window functions. Our signal regionneglects any C ℓ above 706. At high energies (right panel),the effect of the beam window function is more modest,even up to ℓ = 2000 (see Fig. 3). In principle, for highenergies, we could consider a signal region in multipolethat extends to smaller scales. However, we prefer towork with a window in multipole that is independent ofthe energy bin and, therefore, we choose the value of ℓ = 706 as a reasonable compromise.Note that each individual data point in Figs. 6 and 7can be negative, since our auto-APS estimator quantifiesthe excess of power with respect to the photon noise C N . We fit the auto-APS (between ℓ = 49 and 706)in each energy bin to a constant value, in order todetermine the Poissonian C P (the possibility of a non-constant C ℓ is considered later). The fit is performedas discussed in the previous section. The best-fit C P C P as a pink band. Thesignificance of the measured Poissonian auto-APS can bequantified by computing the Test Statistics (TS) of thebest-fit C P , defined as the difference between the − L of the best fit and the − L of the null hypothesis.The latter is obtained from Eq. 10 by setting C P tozero. Assuming Wilks’ theorem, TS is distributed as0 Multipole s r ] - s r - s - [ c m l C -3-2-10123 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -3-2-10123 -18 × Energy bin [1.38-1.99] GeV Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -10-8-6-4-20246810 -21 × Energy bin [50.00-95.27] GeV
FIG. 6. Auto-APS of the IGRB for 2 representative energy bins (between 1.38 and 1.99 GeV in the left panel and between 50.0and 95.27 GeV in the right panel) and for the reference data set (P7REP ULTRACLEAN V15 front events) using the referencemask which excludes | b | < ◦ and 3FGL sources (red circles). The blue triangles show the same but masking the sources in2FGL. Data have been binned as described in Sec. IV A. The solid red line shows the best-fit C P for the red data points, withthe pink band indicating its 68% CL error. The dashed blue line corresponds to the best-fit C P for the blue data points. Notethat only the results in our signal region (i.e. between ℓ = 49 and 706) are plotted and that the scale of the y -axis varies inthe two panels. Also, the blue triangles have been slightly shifted horizontally with respect to the red circles to increase thereadibility of the plots. This will happen also in many of the following plots. Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -10-8-6-4-20246810 -18 × Energy bin [1.38-1.99] GeV Multipole s r ] - s r - s - [ c m l C -20-15-10-505101520 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -20-15-10-505101520 -21 × Energy bin [50.00-95.27] GeV
FIG. 7. Same as Fig. 6, but showing a wider range in multipole, going from ℓ = 10 to 2000. The two dashed grey vertical linesindicate the lower and upper bounds of the multipole range used for the present analysis. Note the different scale of the y -axisin each panel. χ distribution with 1 degree of freedom and, thus, itcan be used to estimate the significance associated to C P . For the default data set masking 3FGL sources, thesignificance of the measured auto-APS C P is larger than3 σ for all energy bins up to 21.8 GeV, except between5.00 and 10.45 GeV. The significance of the detectionis reported in italics in Tabs. I and II. In the case ofthe mask around 3FGL sources, the highest significancein the auto-APS is 6.3 σ and it is reached in the secondenergy bin, i.e. between 0.72 and 1.04 GeV. The way the auto- and cross-APS depend on the energy(i.e. the so-called “anisotropy energy spectrum”) isan informative observable that can provide insight intothe emission causing the anisotropic signal. In fact,in the case that the auto-APS is produced by a singlepopulation of sources, the anisotropy energy spectrumallows their energy spectrum to be reconstructed [28, 43,1 Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 Masking sources in 3FGLMasking sources in 2FGL
FIG. 8. Anisotropy energy spectra for the auto-APS using thereference data set with the default 3FGL mask (red circles)in comparison with the case in which we use the default maskaround 2FGL sources (blue triangles). . If more than one class of objects are responsible forthe signal, then, by detecting features in the anisotropyenergy spectrum, it may be possible to identify energyregimes where the different classes dominate the signal.The measured anisotropy energy spectrum for theauto-APS is shown in Fig. 8. In the figure, the datapoints are weighted by E / ∆ E where E is the log-center of the energy bin and ∆ E is the width of thebin. This weighting is introduced in order to comparethe anisotropy energy spectrum directly with the squaredintensity energy spectrum of the sources responsible forthe anisotropy signal. Fig. 8 compares the auto-APS C P for the case of the mask excluding 3FGL sources (redcircles) to that of the mask excluding 2FGL sources (bluetriangles). As already mentioned, the amplitude of theauto-APS is lower when we exclude the sources in 3FGL.In both data sets, the low-energy part of the spectrumappears generally consistent with a power law, while afeature is apparent around 7 GeV. We comment furtheron the structure of the anisotropy energy spectrum inSec. VI. B. Cross-correlation angular power spectra
Two examples of the cross-APS between energy binsare shown in Fig. 9. The left panel is for the cross-APS between bins at low energies. A clear correlationis detected in the multipole range of interest (boundedby the vertical grey lines in the figure). Note the effect The anisotropy energy spectrum traces the intensity energyspectrum of the sources responsible for the anisotropy signal onlyif the clustering of the source population is independent of energy. of the beam window function on the error bars at highmultipoles, as in Fig. 7. The right panel shows the cross-APS between two high-energy bins. This combinationdoes not correspond to a significant detection, as thebest-fit C P is compatible with zero at a 2 σ level.The best-fit C P for the cross-APS between the i -thand the j -th energy bins are shown in Appendix C,multiplied by E i E j / ∆ E i ∆ E j and for all the possiblecombinations of energy bins. Cross-APS C P is detectedin most combinations of energy bins, with the ones failingto yield a detection mainly involving the two highestenergy bins. Tabs. I and II report the detected cross-APSwith their significance The largest detection significanceis 7.8 σ for the case of the cross-APS between the energybin from 1.99 and 3.15 GeV and the energy bin between3.15 and 5.0 GeV.The tables also report in bold the χ associated withthe best-fit C P according to the definition in Eq. 9.Fig. 10 shows the distribution of the 91 χ of best-fit C P in the 91 independent combinations of the 13 energy bins.The solid black line refers to the case when all sourcesin the 3FGL are masked and the dashed blue line whenonly sources in the 2FGL are masked. Both distributionsare compatible with that of a χ distribution with 9degrees of freedom (i.e. the 10 data points inside thesignal region in multipole minus 1 fitted parameter). Thelatter is represented by a solid red line in Fig. 10. Only3 (4) combinations of energy bins have a χ larger than16.9 (that would correspond to a p -value of 0.05) whenmasking 3FGL (2FGL) sources.Together with the auto-APS in Fig. 8, the cross-APS provides an important handle to characterize theemission responsible for the anisotropy signal. Inparticular, if the latter is due to only one class ofunresolved sources, the auto-APS C i,i P allows us toreconstruct their energy spectrum and the cross-APS canbe predicted as C i,j P = q C i,i P C j,j P . Alternatively, if wedefine the so-called cross-correlation coefficients r i,j as C i,j P / q C i,i P C j,j P , any deviation from 1 when i = j canbe interpreted as an indication of multiple source classescontributing to the signal. In Fig. 11, we show thecross-correlation coefficients corresponding to the best-fit C i,j P for the data set obtained masking 2FGL sources(left panel) and masking 3FGL sources (right panel).In the former case, it is clear that the cross-correlationcoefficients of low-energy bins are systematically smallerthan 1, when correlated with high-energy bins. This isin qualitative agreement with the findings of Ref. [16], inwhich the auto-APS measured in Ref. [1] was explainedby the sum of two different populations of unresolvedblazars at low energies, while, above ∼
10 GeV, the signalwas compatible with only one source class. Figs. 33 and Note that in some cases the best-fit C P is negative. However,whenever that happens the estimated error is large and themeasurement is compatible with zero. Multipole s r ] - s r - s - [ c m l C -4-3-2-101234 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -4-3-2-101234 -18 × =[0.72-1.04] GeV j =[0.5-0.72] GeV, E i E Multipole s r ] - s r - s - [ c m l C -1-0.500.51 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -1-0.500.51 -21 × =[199.05-500] GeV j =[50.00-95.27] GeV, E i E FIG. 9. Cross-APS for two representative combinations of energy bins as indicated above the two panels. The default data setand 3FGL mask are used to compute the red circles, while the blue triangles are for the mask excluding 2FGL sources. Thevertical dashed grey lines denote the bounds of the multipole range used in the analysis. The best-fit C P is shown as a solidred line in the case of the mask around 3FGL sources, with a pink band indicating its 68% CL error. The dashed blue linecorresponds to the best-fit C P for the blue data points (i.e. masking 2FGL sources). (best-fit) χ Masking sources in 3FGLMasking sources in 2FGL χ FIG. 10. Normalized distribution of the χ (defined as inEq. 9) for the best-fit Poissonian C P for all 91 independentcombinations of energy bins. The solid black line is for thecase when 3FGL sources are masked and the dashed blue lineis for the mask covering 2FGL sources. The solid red curve isa χ distribution with 9 degrees of freedom.
34 in Appendix D show, for each energy bin i , how thecross-correlation coefficents r i,j depend on energy E j .When 3FGL sources are masked (right panel) thesituation is less clear as errors are larger (especially athigh energies) and the estimated C P more uncertain. Wefurther discuss about the nature of our auto- and cross-APS in Sec. VI. Note that in some cases the coefficients r i,j shown in Fig. 11 are larger than 1, since only the best-fit values are plotted. They are, however, compatible with 1, within their uncertainty. Also, some coefficientsare negative (and they are associated with a black pixel).Although within the error bars these negative r i,j areactually compatible with 0, we note that negative valuesare allowed, in the case of anti-correlations between twoenergy bins.We finish this section by studying whether the binnedauto- and cross-APS C ℓ are better described by an APSthat changes with multipole, instead of a constant value.We fit the binned C ℓ with a power law , i.e., C ℓ = A ( ℓ/ℓ ) − α , with ℓ = 100. We leave the normalization A free to vary independently in all the 91 combinations ofenergy bins but we consider one common slope , i.e. α .The best-fit value of α is − . ± .
08 and it correspondsto a χ per degree of freedom of 0.91. This should becompared to the global (i.e. for all the 91 combinationsof energy bins) Poissonian fit, with also has a χ perdegree of freedom of 0.91. Therefore, we cannot deduceany preference for the power-law scenario. The fit is performed with
Minuit2 v5.34.14,http://lcgapp.cern.ch/project/cls/work-packages/mathlibs/minuit/index.html If the auto- and cross-APS are interpreted as produced bya population of unresolved sources, they can be expressed interms of the 3-dimensional power spectrum of the density fieldassociated with the sources of the gamma-ray emission [5, 45, 46].The latter determines the dependence on multipole, hence theshape of the auto- and cross-APS. Normally, the 3-dimensionalpower spectrum is only mildly dependent on the gamma-rayenergy, which is encoded in the “window function”. Therefore,the APS associated with different combinations of energy binsis expected to have approximately an energy-independent shape.It is therefore reasonable to assume a constant α . TABLE I. Best-fit Poisson auto- and cross-APS C P for the default data set and for the mask covering all sources in 3FGL, in units of cm − s − sr − . The numbers initalics indicate the significance of the detection in units of standard deviation (see text), while the numbers in bold give the χ associated with the corresponding C P .The entries marked in grey correspond to the auto-APS. Energy 0.50- 0.72- 1.04- 1.38- 1.99- 3.15- 5.00- 7.23- 10.45- 21.83- 50.00- 95.27- 199.05-bin [GeV] 0.72 1.04 1.38 1.99 3.15 5.00 7.23 10.45 21.83 50.00 95.27 199.05 500.000.50- (7 . ± . × − , . ± .
62) (3 . ± . × − × − , , . ± .
34) (1 . ± .
21) (5 . ± . × − × − × − , , , . ± .
86) (1 . ± .
16) (3 . ± .
88) (4 . ± . × − × − × − × − , , , , . ± .
21) (8 . ± .
18) (3 . ± .
62) (2 . ± .
48) (2 . ± . × − × − × − × − × − , , , , , . ± .
32) (3 . ± .
75) (1 . ± .
39) (1 . ± .
29) (1 . ± .
19) (8 . ± . × − × − × − × − × − × − , , , , , , . ± .
82) (1 . ± .
46) (9 . ± .
35) (6 . ± .
73) (5 . ± .
15) (2 . ± .
64) (6 . ± . × − × − × − × − × − × − × − , , , , , , , . ± .
06) (6 . ± .
36) (5 . ± .
69) (1 . ± .
26) (3 . ± .
81) (2 . ± .
46) (1 . ± .
25) (3 . ± . × − × − × − × − × − × − × − × − , , , , , , , , . ± .
58) (8 . ± .
18) (4 . ± .
60) (4 . ± .
16) (3 . ± .
75) (2 . ± .
41) (1 . ± .
23) (5 . ± .
56) (1 . ± . × − × − × − × − × − × − × − × − × − , , , , , , , , , . ± .
35) (4 . ± .
82) (2 . ± .
91) (2 . ± .
66) (1 . ± .
42) (1 . ± .
23) (3 . ± .
31) (3 . ± .
88) (2 . ± .
79) (1 . ± . × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , . ± .
69) ( − . ± .
93) (7 . ± .
63) (1 . ± .
34) (1 . ± .
22) (0 . ± .
18) (1 . ± .
66) (1 . ± .
44) (1 . ± .
39) (7 . ± .
22) (2 . ± . × − × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , , . ± .
04) (1 . ± .
58) (1 . ± .
88) (0 . ± .
10) (0 . ± .
35) (1 . ± .
48) (2 . ± .
15) (1 . ± .
84) ( − . ± .
52) (3 . ± .
38) (3 . ± .
09) (4 . ± . × − × − × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , , , . ± .
16) ( − . ± .
43) (0 . ± .
71) ( − . ± .
26) ( − . ± .
02) (2 . ± .
44) ( − . ± .
51) (1 . ± .
71) (2 . ± .
52) (1 . ± .
33) ( − . ± .
32) ( − . ± .
74) (0 . ± . × − × − × − × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , , , , Energy 0.50- 0.72- 1.04- 1.38- 1.99- 3.15- 5.00- 7.23- 10.45- 21.83- 50.00- 95.27- 199.05-bin [GeV] 0.72 1.04 1.38 1.99 3.15 5.00 7.23 10.45 21.83 50.00 95.27 199.05 500.000.50- (1 . ± . × − , . ± .
07) (6 . ± . × − × − , , . ± .
38) (3 . ± .
22) (1 . ± . × − × − × − , , , . ± .
30) (2 . ± .
19) (1 . ± .
09) (1 . ± . × − × − × − × − , , , , . ± .
22) (2 . ± .
13) (9 . ± .
66) (8 . ± .
49) (6 . ± . × − × − × − × − × − , , , , , . ± .
14) (1 . ± .
08) (5 . ± .
40) (4 . ± .
30) (3 . ± .
20) (2 . ± . × − × − × − × − × − × − , , , , , , . ± .
87) (3 . ± .
47) (2 . ± .
25) (1 . ± .
18) (1 . ± .
12) (8 . ± .
68) (3 . ± . × − × − × − × − × − × − × − , , , , , , , . ± .
63) (2 . ± .
35) (1 . ± .
18) (1 . ± .
13) (1 . ± .
08) (6 . ± .
45) (2 . ± .
25) (1 . ± . × − × − × − × − × − × − × − × − , , , , , , , , . ± .
60) (2 . ± .
33) (1 . ± .
17) (1 . ± .
12) (9 . ± .
77) (6 . ± .
42) (2 . ± .
23) (1 . ± .
16) (2 . ± . × − × − × − × − × − × − × − × − × − , , , , , , , , , . ± .
42) (6 . ± .
89) (3 . ± .
94) (3 . ± .
69) (3 . ± .
42) (2 . ± .
23) (9 . ± .
27) (6 . ± .
85) (8 . ± .
75) (4 . ± . × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , . ± .
77) (1 . ± .
97) (1 . ± .
48) (1 . ± .
35) (1 . ± .
22) (6 . ± .
21) (4 . ± .
67) (3 . ± .
45) (3 . ± .
41) (1 . ± .
23) (9 . ± . × − × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , , . ± .
10) (1 . ± .
61) (6 . ± .
02) (8 . ± .
22) (3 . ± .
42) (3 . ± .
77) (1 . ± .
43) (1 . ± .
29) (1 . ± .
26) (7 . ± .
43) (3 . ± .
74) (1 . ± . × − × − × − × − × − × − × − × − × − × − × − × − , , , , , , , , , , , , − . ± .
74) ( − . ± .
14) (0 . ± .
57) (0 . ± .
17) ( − . ± .
44) (8 . ± .
08) (1 . ± .
27) (2 . ± .
54) (4 . ± .
38) (4 . ± .
58) (4 . ± .
95) ( − . ± .
54) (0 . ± . × − × − × − × − × − × − × − × − × − × − × − × − × − , , < , , , , , , , , , , < , Energy [GeV]1 10 E ne r g y [ G e V ]
10 00.20.40.60.811.21.4
Masking 2FGL sources
Energy [GeV]1 10 E ne r g y [ G e V ]
10 00.20.40.60.811.21.41.61.8
Masking 3FGL sources
FIG. 11. Cross-correlation coefficients between energy bins. Each pixel in the panels corresponds to a pair ( i, j ) of energy binsand it is colored according to the cross-correlation coefficent r i,j . By construction the panels are symmetric with respect to thediagonal. The panel on the left refers to the default data set with a mask that covers the sources in 2FGL, while the one on theright is for the mask covering 3FGL sources. Cross-correlation coefficents below 1 indicate that the signal is due to multiplepopulations of sources. C. Validation studies
We note that the uncertainties reported in thelast section only include statistical errors. It istherefore important to estimate any systematic errorsas, e.g., those related to the analysis (such as theforeground cleaning and the use of the mask) or to thecharacterization of the instrument, which may affect theeffective area and beam window functions. We discusspossible sources of systematic errors in the followingsections.
1. Foreground cleaning
The Galactic diffuse emission is bright, especially atlow energies, and generally displays an approximatesymmetry around the disk of the Galaxy, leading toexcess power at low multipoles, corresponding to largeangular scales. The measured auto-APS calculatedboth with and without performing foreground cleaningis shown in the left panel of Fig. 12 (red dots andblue triangles, respectively) for a selected energy binat low energy. The default mask and data set areused. The effect of foreground cleaning is dramatic atlow multipoles, significantly reducing the measured C ℓ below ℓ ∼
50. On the other hand, our analysis onlyconsiders multipoles larger than ℓ = 49, where the effectof foreground cleaning is smaller, although still important enough to be non-negligible. Above ℓ ∼ ℓ between 49 and 706(red circles, the same as in Fig.8), is compared to thebest-fit C P for the case without foreground cleaning butperforming the fit only between ℓ = 143 and 706 (bluetriangles), i.e. neglecting the first four bins in multipoleinside our signal region. Errors at low energies are largerfor the uncleaned case than for the cleaned one. Thisis due to the fact that, at low energies, only the few C ℓ with ℓ < ∼
300 play a role in the determination of thebest-fit C P , since at larger ℓ the beam suppression is toostrong. Therefore, cutting the signal region at ℓ = 143means that the best-fit C P is determined only by veryfew data. Nonetheless, the two cases are found to be ingood agreement within their uncertainties at all energies.From this we can conclude that the foreground cleaningis effective even down to ℓ ∼
50, therefore validating ourchoice for the signal region in multipole.Also, from the left panel of Fig. 12, it is clearthat the binned APS C ℓ without foreground cleaningis characterized by much larger error bars than withcleaning, at least for ℓ < ∼ C ℓ aresimply the square root of the diagonal elements of thecovariance matrix, while the full covariance matrix is6 Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × with foreground cleaningwithout foreground cleaningPoissonian fit (with cleaning) -10-8-6-4-20246810 -18 × Energy bin [1.04-1.38] GeV Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -16 -19 -18 -17 -16 (49,706) ∈ With foreground cleaning, l (143,706) ∈ Without foreground cleaning, l
FIG. 12.
Left:
Auto-APS in the energy bin between 1.04 and 1.38 GeV, comparing the data with (red circles) and without(blue triangles) foreground cleaning. The solid red line indicates the best-fit C P for the case with foreground cleaning, and thepink band its 68% CL error. The two dashed vertical lines mark our signal region in multipole. Right:
Poissonian auto-APSas a function of the energy for the case with foreground cleaning and a signal region between ℓ = 49 and 706 (red circles) andfor the case without foreground cleaning and a signal region between ℓ = 143 and 706 (blue triangles). Default data selectionand 3FGL mask are used. Multipole M u l t i po l e -3 -2 -1 With foreground cleaning
Multipole M u l t i po l e -3 -2 -1 Without foreground cleaning
FIG. 13. Normalized covariance matrix ( σ i,j / √ σ i,i σ j,j ) of the binned C ℓ shown in the left panel of Fig. 12, i.e., for the energybin between 1.04 and 1.38 GeV, default data selection and default mask covering 3FGL sources. The left panel shows the casewith foreground cleaning, while the right panel is for the uncleaned case. The comparison between the two panels indicatesthat large covariances are present in the case without foreground cleaning up to multipoles ℓ ∼ plotted in Fig. 13, for the data between 1.04 and 1.38GeV, with (left panel) and without foreground cleaning(right panel). Each pixel in the panels correspondsto a pair of bins in multipole and its color providesthe covariance between those two bins. We do notdirectly plot the covariance matrix but, instead, eachelement σ i,j is divided by the square root of the product of the corresponding diagonal elements √ σ i,i σ j,j .The main difference between the two panels is at lowmultipoles, where the case without foreground cleaning ischaracterized by a large covariance among different bins.This large covariance causes the diagonal terms at ℓ < ∼ ℓ < ∼
30 are characterized by larger C ℓ ℓ ∼ ℓ < ∼
50 and itconsiderably removes the coupling between multipoles,leading to smaller and weakly correlated estimatederrors. It also justifies the use of Eqs. 9 and 10 forthe determination of the Poissonian APS, since theyare valid only under the hypothesis that covariances arenegligible .
2. Data selection
Next we consider the impact of our choice ofdata set. As described in Sec. II, our analysisis based on P7REP ULTRACLEAN V15 front events.For comparison, we now show the results for twodifferent event selections using Pass 8 data, i.e., themost recent revision of the event-level
Fermi
LATreconstruction analysis [47]. In particular, we will usethe P8R2 ULTRACLEANVETO V6 class, designed toreduce the cosmic-ray contamination significantly. Weconsider separately two event selections, i.e. only Pass8 front-converting events and the so-called PSF3 events.PSF3 is a new selection available with Pass 8 data andit is characterized by an improved angular resolution.The effective area for the PSF3 events is roughly afactor of two smaller than that for the front events,since PSF3 represents the quartile of events with thebest angular resolution, while the front events constituteapproximately half the total events gathered by the LAT.The same analysis pipeline applied to the Pass 7 datais employed to the Pass 8 data, including foregroundcleaning with the same Galactic diffuse model (refittedto the Pass 8 events outside the mask).In Fig. 14 we compare the measured auto-APS in oneenergy bin between the default data set (i.e., Pass 7,denoted by red circles) and the two Pass 8 selections:Pass 8 front-converting events in the left panel (orangesquares) and Pass 8 PSF3 in the right panel (bluetriangles). The auto-APS is shown over the multipolerange between ℓ = 10 and 2000. The two Pass 8 datasets are overall in good agreement with the default Pass 7data set in the multipole range used for analysis, markedby the two grey vertical dashed lines in the figure.In Fig. 15 we show the anisotropy energy spectrafor the three data sets discussed above. TheirPoissonian auto-APS agree well within the measurementuncertainties in the various energy bins. The sharp dropin C P around ∼ Note that in Ref. [1] the covariance between multipole bins wasnot discussed. significant in the Pass 8 PSF3 data and absent in thePass 8 front data, suggesting that the feature in the Pass7 data may be the result of a statistical fluctuation. Also,with Pass 8, the auto-APS around 70 GeV has a largervalue than with Pass 7, although the difference is onlyat the 2 σ level and, thus, not very significant. We stressthat this is only a qualitative comparison and a morethorough analysis of the Pass 8 data should be performed.With Pass 8, the measurement of the auto-APS andcross-APS is expected to improve in several ways, e.g.taking advantange of the new PSF classes (from PSF0 toPSF3), especially at low energies where the measurementuncertainties in the Pass 7 data are dominated by thesuppression induced by the beam window functions and(potentially) by the leaking from bright sources outsidethe mask (see Sec.V C 3). Also, new data selections areavailable with Pass 8, characterized by different balancesbetween effective area and cosmic-ray contamination.In fact, the improvement expected from using Pass 8PSF3 or Pass 8 front data can already been seen in thereduction of the error bars for the blue triangles andorange squared in Fig. 15, with respect to the red circles,especially at around 100 GeV. A detailed study with Pass8 is beyond the scope of the present analysis and is leftfor future work.We further investigate the impact of event selectionby comparing the results obtained from the Pass 7 datausing front data only (i.e., our default choice) to theresults obtained using both the front and back data.Including back-converting events in the analysis has theadvantage of increasing the statistics by enlarging theeffective area by a factor of ∼
2. However, the averagePSF for the front+back data set is poorer than for thefront events alone, leading to a larger suppression dueto the beam window function and to a stronger leakageoutside the mask from bright point-like sources. In thiscomparison it should be kept in mind that the data setsare not independent since, by definition, the front+backdata set contains all the front-converting events. Also,it is important to note that due to the poorer PSF ofthe front+back data set, our source-masking scheme maynot be sufficiently effective for that data set, particularlyat low energies where the PSF is broadest (see also thediscussion in Sec. V C 3).The left panel of Fig. 16 shows the auto-APS C ℓ ina specific energy bin. Red circles refer to the Pass7 front data set and the blue triangles to the Pass 7front+back one. The right panel indicates the Poissonianauto-APS as a function of energy, with the same colorcode. The measured C ℓ is in good agreement betweenthe two data sets at all multipoles in our signal region.The same is true for the Poissonian C P , except in thelowest energy bin, where the front+back data yields asignificantly higher C P . This discrepancy is consistentwith the possibility that, for the front+back data set, themask employed here (covering all sources in the 3FGL) isnot big enough to get rid of the emission of the sources atlow energies. Note that, also in this case, the significance8 Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × Pass 7 (front)Pass 8 (ULTRACLEAN front)Poissonian fit (Pass 7 front)Poissonian fit (Pass 8, ULTRACLEAN front) -10-8-6-4-20246810 -18 × Energy bin [1.04-1.38] GeV Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × Pass 7 (front)Pass 8 (ULTRACLEAN PSF3)Poissonian fit (Pass 7, front)Poissonian fit (Pass 8, ULTRACLEAN PSF3) -10-8-6-4-20246810 -18 × Energy bin [1.04-1.38] GeV
FIG. 14.
Left:
Comparison of the auto-APS measurement in the energy bin between 1.04 and 1.38 GeV between the defaultPass 7 data set (red circles) and the Pass 8 front event selection (orange squares).
Right: same as the left panel but thecomparison is between the Pass 7 data (red circles) and the Pass 8 PSF3 event selection (blue triangles). The solid red linemarks the best-fit C P for the default Pass 7 data set, with the pink band indicating its 68% CL error. The dashed orange(blue) line gives the best-fit C P for Pass 8 front (PSF3) in the left (right) panel. The vertical grey dashed lines mark the signalregion between ℓ = 49 and 706. The default mask covering 3FGL sources is applied. Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 Pass 7 (front)Pass 8 (ULTRACLEAN PSF3)Pass 8 (ULTRACLEAN front)
FIG. 15. Poissonian auto-APS as a function of energy for thePass 7 (red circles), Pass 8 front (orange squares) and Pass 8PSF3 (blue triangles) data sets. The default mask covering3FGL sources is applied. of the dip at ∼
3. Mask around resolved sources
We now investigate the effect of any possible leakage ofemission outside the mask around the resolved sources.We recall that our default mask excludes (in addition toa latitude cut of | b | < ◦ ) a disk with a radius of 3.5 ◦ around the 500 brightest sources in the 3FGL catalog, a disk with a radius of 2.5 ◦ around the following 500sources, one with a radius of 1.5 ◦ for the following 1000sources and, finally, a disk with a region with a 1.0 ◦ -radius around all the remaining ones. This is what werefer to as our default mask when covering the sourcesin 3FGL. However, in order to validate our choice, weconsider four additional masks. They are defined asfollows: • ◦ -mask excludes a disk with a radius of 4 ◦ aroundthe 500 brightest sources in the 3FGL catalog,a disk with a radius of 3 ◦ for the following 500sources, one with a radius of 1.5 ◦ for the next1000 sources and one with a radius of 1 ◦ for theremaining ones; • ◦ -mask excludes a disk with a radius of 3.5 ◦ around the 500 brightest sources in the 3FGLcatalog, a disk with a radius of 2.5 ◦ for the following500 sources, one with a radius of 2.0 ◦ for the next1000 sources and a disk with a radius of 1.5 ◦ forthe remaining ones; • ◦ -mask covers a disk with a radius of 2.0 ◦ aroundthe 500 brightest sources and a disk with a radiusof 1.0 ◦ around the remaining sources; • ◦ -mask excludes a disk with a radius of 1.0 ◦ around each source.Our default mask is located between the 2 ◦ -mask andthe 3.5 ◦ -mask, in terms of masked area. The specificdetails of the masks considered are not the result of an9 Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × frontfront and backPoissonian fit (front)Poissonian fit (front and back) -10-8-6-4-20246810 -18 × Energy bin [1.04-1.38] GeV Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 frontfront and back FIG. 16.
Left:
Comparison of the auto-APS in the energy bin between 1.04 and 1.38 GeV between the default data set whichuses only front events (red circles) and the front+back data set (blue triangles). The solid red line marks the best-fit C P forthe default front-converting data, with the pink band indicating its 68% CL error. The dashed blue line gives the best-fit C P for the front+back data set. The vertical grey dashed lines mark the signal region between ℓ = 49 and 706. Right:
Poissonianauto-APS as a function of energy for the front events (red circles) and the front+back ones (blue triangles). The default maskcovering 3FGL sources is applied for both panels.
Multipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -18 × default mask-mask ° ° ° Poissonian fit (2 -mask) ° Poissonian fit (1 -10-8-6-4-20246810 -18 × Energy bin [1.04-1.38] GeV Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 default mask-mask ° ° FIG. 17.
Left:
Comparison of the auto-APS in the energy bin from 1.04 to 1.38 GeV among the case with the default maskcovering the sources in 3FGL (red circles), the case with the 2 ◦ -mask (blue triangles) and the one with the 1 ◦ -mask (orangesquares). The solid red line marks the best-fit C P for the default mask, with the pink band indicating its 68% CL error. Thelong-dashed blue line gives the best-fit C P when the 2 ◦ -mask is employed and the short-dashed orange one for the case withthe 1 ◦ -mask. The vertical grey dashed lines mark the signal region between ℓ = 49 and 706. Right:
Poissonian auto-APS asa function of energy for the case with the default mask (red circles), the 2 ◦ -mask (blue triangles) and the 3.5 ◦ -mask (greensquares). a-priori analysis and, thus, they are somewhat subjective.However, our goal is to identify a reasonable mask that isas small as possible without suffering from leakage frompoint-like sources. As proved in the following, our defaultmask provides a suitable choice.Above few GeV, where the PSF is narrower, we expectthe 1 ◦ -mask to be sufficient to exclude the emission of thesources detected in 3FGL. However, at low energies some leakage may appear. Results are summarized in Fig. 17.The left panel shows the measured auto-APS in the energy bin between 1.04 and 1.38 Gev, for the 1.0 ◦ -mask(orange squares), for the 2.0 ◦ -mask (blue triangles) andfor the default one (red circles). It is clear that there isa significant contamination due to power leakage outsidethe 1.0 ◦ -mask, especially at ℓ <
50, but up to ℓ ∼
80. Theother two more aggresive masks give consistent results in0this energy bin. In the right panel, we plot the anisotropyenergy spectrum for the 2 ◦ -masks (blue triangle), forthe default one (red circle) and for the 3.5 ◦ -mask (greensquares). While, at high energies, the three cases yieldconsistent results, the 2 ◦ -mask shows still an excess ofpower in the first energy bin. On the other hand, resultsfor the 3.5 ◦ -mask are consistent with our default mask.The anisotropy energy spectrum for the 4 ◦ -mask (notshown in Fig. 17 for clarity) is also consistent with thedefault case. This validates our choice of the latter asour fiducial mask when dealing with 3FGL sources .A similar validation is performed on the mask coveringthe sources in 2FGL. We find that cutting a 1 ◦ -diskaround all 2FGL sources leads to some power excess atlow energies. However, extending the mask by coveringa disk with a radius of 2 ◦ for all sources is enough toget rid of the leakage and there is no need of moreaggressive masks as for the case of 3FGL sources. Thisis probably due to the fact that, when masking 2FGLsources, the measured power spectra are intrinsicallylarger than when masking sources in 3FGL (see Sec. V).Thus, the contimation from leakage has relatively a minorimpact. D. Effect of the gamma-ray emission from the Sun
Steady gamma-ray emission from the Sun was detectedin the
Fermi
LAT data in Ref. [48] from 0.1 to 10GeV. Later, Ref. [49] extended the detection up to 100GeV, also establishing that the flux varies with timeand anti-correlates with Solar activity. Gamma rays areproduced from the interaction of cosmic rays with theSolar atmosphere [50], as well as from Inverse-Compton(IC) scattering of cosmic-ray electrons and positrons withSolar photons [51–53].The emission is quite difficult to see by eye because,even if quite significant, it is spread over the path followedby the Sun in the sky, i.e. the ecliptic. However, itmay still induce some feature in the auto- and cross-APS. We test this possibility by masking the region of1 . ◦ above and below the ecliptic. The auto- and cross-APS obtained after having introduced this additionalmask are compatible with our default case within theiruncertainty. Thus, we conclude that the effect of the Sunon the measured anisotropies is negligible . We also test an additional mask that covers exactly the sameregion of sky as our default mask for 3FGL sources but it alsomasks the region around Loop-I and the Galactic Lobes. Thebest-fit C P with this more aggressive cut are compatible with thedefault Poissonian C P in Fig. 8, within their statistical errors. The high-energy emission of the Moon peaked at about 200 MeV,with a similar intensity than the Sun [54]. However, at higherenergies, it has an energy spectrum that is steeper than that ofthe Sun. Therefore, above 500 MeV, the effect of the Moon onthe APS measurement is expected to be negligible.
Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -16 -19 -18 -17 -16 This work (new energy bins)This work (old energy bins)Ackermann et al. (2012)
FIG. 18. Poissonian auto-APS as a function of energy forthe default data set in this analysis, with 13 energy bins(blue triangles) and with the 4 energy bins used in Ref. [1](red circles). Note that data are obtained using the maskthat covers a 2 ◦ -circle around each source in 2FGL. The greysquares denote the measurement from Ref.[1] using the samemask. E. Comparison with previous measurement
We conclude this section by comparing our newmeasurement to the previous (indeed, the first)anisotropy measurement from Ref. [1]. Our currentanalysis includes many improvements with respect to theoriginal one, both from the perspective of the data set(as we now use Pass 7 Reprocessed events and IRFs,compared to the Pass 6 events used in Ref. [1]) andin terms of analysis method, including an improvedcalculation of the noise term C N , the deconvolution ofthe mask (performed now with PolSpice ) and a MC-validated procedure to bin the auto-APS in multipolesand to estimate its error. The improved data setalso allows us to measure the auto-APS with betterprecision over a larger multipole range covering thewindow between ℓ = 49 and 706, while the analysisin Ref. [1] was restricted to ℓ = 155 − C P in ourcurrent measurement compared to the original one, we1find good consistency with Ref. [1]. The only exceptionis the highest energy bin of the original analysis, which islower than the current measurement and inconsistent atabout 3 σ . Many factors may lead to the small systematicincrease of the new C P in the first 3 bins and to thelarger difference in the last energy bin. However, weattribute this trend primarily to the way the data arebinned in multipole and to the way the Poissonian fit C P is determined. As discussed in Sec. IV A, in thisanalysis we follow a different procedure with respect tothe original analysis in Ref. [1], after having verified thatthe latter can lead to a downward bias of both the C ℓ and the best-fit C P .We end by noting that a concern about the auto-APSin Ref. [1] being somewhat underestimated was raisedalready in Ref. [55]. However, in that case it was claimedthat the correct anisotropy should have been a factor 5-6larger than the measured one for each energy. In the lightof the present analysis, this is true only for the highestenergy bin, while for the others the difference is only of20-30%, and not significant within error bars. VI. INTERPRETATION IN TERMS OFSOURCE POPULATIONS
In this section we provide a phenomenologicalinterpretation of our measurement in terms of differentpopulations of unresolved sources. The main observablesthat we consider are the results of the Poissonian fitsto the auto- and cross-APS, i.e., the anisotropy energyspectrum.We assume the sources responsible for the signalto be point-like and unclustered [7, 10, 16, 19], andthat they give rise only to a Poissonian auto- andcross-APS. We also assume each population to becharacterized by a common intensity energy spectrum F a ( E ). The index a runs over the number of sourcepopulations contributing to the signal. The contributionof the a -th source population to the auto-APS will beproportional to F a ( E i ) , while it will be proportionalto F a ( E i ) F a ( E j ) in the case of the cross-APS. Ourchoice of interpreting the auto- and cross-APS data ina phenomenological way is motivated by the desire to bemodel-independent. Alternative interpretations in termsof physically-motivated models of astrophysical gamma-ray emitters are on-going. We start by considering onesingle source population with a power-law spectrum, i.e. F ( E ) = A (cid:18) EE (cid:19) − α , (13)with E = 1 GeV. We fit both the best-fit Poissonianauto- and cross-APS taken from Tab. I, i.e. for the maskaround 3FGL sources. This scan and the following onesare performed with MultiNest − in order to provide agood sampling of the likelihood. The prior probability is chosen to be flat for all free parameters, between -15.0 and -5.0 for log ( A ) (and all the normalizations,measured in cm − s − sr − ), between 0.0 and 5.0 for theslopes and between 5.0 GeV and 500.0 GeV for the energybreaks (see later). The best-fit solution is reported inTab. III and is represented by a solid magenta line in theleft panel of Fig. 19. The best fit has a χ per degreeof freedom which is 1.52, corresponding to a p -value of0.001.Alternatively, we also consider a broken power lawparametrized as follows: F ( E ) = (cid:26) A ( E/E ) − α if E ≥ E b A ( E /E b ) + α − β ( E/E ) − β otherwise . (14)In this case, the best-fit is reported in Tab. III and shownas a solid blue line in the left panel of Fig. 19. Its χ perdegree of freedom is 1.36 with a p -value of 0.01.Then, we allow for the possibility of two independentpopulations, starting with the case of two power laws.The best-fit values are reported in Tab. III and the modelis represented by the solid yellow line in the right panel ofFig. 19: one population explains the data points below afew GeV and another one reproduces the data at higherenergies. The best fit has a χ per degree of freedom of1.47, corresponding to a p -value of 0.003.We also consider the possibility of two broken powerlaws. With a χ per degree of freedom of 1.10, itrepresents the best description to the data. The model isshown as a solid black line in both panels of Fig. 19: onebroken power law reproduces the data at low energies(short-dashed black line) and the other one at higherenergies (long-dashed black line). The best-fit solutionfor the cross-APS is shown in Figs. 31 and 32 inAppendix C.Finally, we also test the hypothesis of one populationemitting as a power law and one as a broken power law.This interpretation is characterized by a χ per degreeof freedom of 1.16 (with a p -value of 0.14) and it isrepresented in the right panel of Fig. 19 by a solid greenline. The fit is slightly worse than the case with twobroken power laws, especially around 3-4 GeV.The difference between the χ of the best-fit solutionfor a model with one population and the same quantityfor the model with two populations can be used asa TS to determine whether we can exclude the one-source-population scenario. From the values of the χ in Tab. III, the exclusion is at 95% CL in all cases .We can test how the different interpretations performalso in a Bayesian framework. Indeed, we can define the In the comparison between one and two populations of sources,if the number of additional degrees of freedom is 2 (as for thecase in which the sources in the second population emit as powerlaws), then the 95% CL exclusion corresponds to a ∆ χ of 5.99.On the other hand, if the second population emits as a brokenpower law, the number of additional degrees of freedom is 4 andthe 95% CL limit is obtained for a ∆ χ of 9.49. Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 Masking sources in 3FGLOne power-lawOne broken power-lawTwo broken power-laws
Energy [GeV]1 10 s r ] - s r - s - c m [ G e V - E ) ∆ ( P C E -19 -18 -17 -19 -18 -17 Masking sources in 3FGLTwo power-lawsOne power-law and one broken power-lawTwo broken power-laws
FIG. 19. Anisotropy energy spectrum for the default data set masking the sources in 3FGL (red circles). The different linescorrespond to the best-fit models to the measured auto- and cross-APS with one or two populations of unresolved sources. Thesolid magenta line and the solid blue one (left panel) are for one population emitting as a power law or as a broken power law,respectively. The solid yellow line (right panel) is for two populations with power-law energy spectra. The solid green line (rightpanel) shows the best-fit in the case of one population emitting as a power law and another as a broken power law. Finally,the thicker solid black line (present in both panels) represents the case of two populations emitting as broken power laws. Thisis the scenario that best fits the data. In this case, the contribution of the two components are shown as short-dashed andlong-dashed black lines.TABLE III. Best-fit values for the parameters defining the populations assumed to describe the measured auto- and cross-APS.See the text for the definition of the parameters. The normalizations ( A , A and A ) are measured in cm − s − sr − and theenergy breaks ( E b , E b, and E b, ) are measured in GeV. Errors are given at 68% CL. The table also indicates the number ofdegrees of freedom N dof (i.e., the number of fitted data points minus the number of free parameters), the χ of the best-fitsolution, the χ of the best-fit point per degree of freedom and the corresponding p -value. N dof χ χ /N dof p -value One power law log ( A ) α − . +0 . − . . +0 . − .
89 135.31 1.52 0.001
One broken power law log ( A ) α β E b − . +0 . − . . +0 . − . > .
74 92 . +16 . − .
87 118.57 1.36 0.010at 68% CL
Two power laws log ( A ) α log ( A ) α − . +0 . − . . +0 . − . − . +0 . − . . +0 . − .
87 127.60 1.47 0.003
Two broken power laws log ( A ) α β E b, log ( A ) α β E b, − . +0 . − . . +0 . − . > .
49 3 . +1 . − . − . +0 . − . . +0 . − . > .
86 84 . +10 . − .
83 91.58 1.10 0.240at 68% CL at 68% CL
One power law and one broken power law log ( A ) α log ( A ) α β E b, − . +0 . − . . . . − . +0 . − . . +0 . − . > . . .
85 98.86 1.16 0.140at 68% CL
Bayes Factor B as the ratio of the so-called “evidence”for two competing models (given the data) and it canbe used to discriminate between them. In particular,with a ln B = 0 .
5, there is not a preference between theinterpretation with one or two power laws (according tothe Jeffrey’s scale [59]), while, with a ln B = 3 .
1, there is a weak preference for the two-broken-power-laws solutionover the one-broken-power-law one.3
VII. SIMULATING THE GAMMA-RAYEMISSION INDUCED BY DARK MATTER
From this section onwards we focus our attentionon the DM-induced gamma-ray emission: we firstsummarize how we simulate this component, and thenwe analyze our mock gamma-ray sky maps by computingtheir auto- and cross-APS. This will constitute ourprediction for the APS associated with DM that willbe compared to the measured auto- and cross-APSpresented in the previous sections.The simulated DM signal needs to account for all DMstructures (halos and subhalos) around us, including theemission generated in the halo of our own Milky Way(MW). We divide the DM auto- and cross-APS intodifferent components that are discussed separately in thefollowing subsections (from Sec. VII A to Sec. VII E). Wefollow closely the semi-analytical procedure developedin Ref. [34], i.e., we directly employ catalogs of DM(sub)halos from N -body simulations and complementthem with well-motivated recipes to account for theemission of DM halos and subhalos below the massresolution of the simulations. As in Ref. [34], we makeuse of the Millennium-II and Aquarius simulations, fromthe Virgo Consortium [60–62] to simulate the Galacticand extragalactic components, respectively.We take particular care in estimating the systematicuncertainties associated with the DM auto- and cross-APS. In particular, each time we introduce a quantitythat is not well determined, we consider a reasonablerange of variability for it and determine its impact onthe final DM signal.We separately consider gamma-ray emission producedby annihilations or decays of DM particles, organizingour predictions in the form of HEALPix maps with
Nside =512. This corresponds to 3145728 pixels and anangular size of approximately 0 . ◦ . The order is lowerthan the one used in the data analysis (see Sec. III).However, note that we will only compare our predictionsfor the DM signal to the measured spectra below ℓ = 706,i.e. for angular scales larger than 0.25 ◦ .The gamma-ray flux (in units of cm − s − ) produced byDM annihilations in the i -th energy bin and coming fromthe pixel centered towards direction n j can be written asfollows:Φ( E i , n j ) = h σ ann v i πm χ Z ∆Ω j d Ω n Z . . dz c (1 + z ) H ( z ) ρ χ ( z, n j ) Z E i +1 E i dE γ dN ann γ ( E γ (1 + z )) dE exp[ − τ EBL ( E γ (1 + z ))] , (15)where the integration d Ω n extends over the pixel centeredon n j . For redshifts higher than ∼
2, the evolution of theDM density field, combined with the larger comovingvolume probed attenuates the signal to a negligiblelevel. The interaction with the Extragalactic BackgroundLight (EBL) additionally reduces the emisson from large redshifts . The EBL attenuation is modeled in Eq. 15by the factor exp( − τ EBL ( E γ (1 + z ))), which is takenfrom Ref. [64]. The thermal average of the cross sectiontimes the relative velocity and the mass of the DMparticles are expressed by h σ ann v i and m χ , while c and H ( z ) are the speed of light and the Hubble parameter.The function ρ χ ( z, n j ) denotes the DM density atredshift z towards the direction n j . The photon yield dN ann γ /dE determines the number of photons producedper annihilation. Different mechanisms of gamma-rayproduction contribute to dN ann γ /dE . We specify whichcontribution is included when we discuss the differentcomponents of the total DM signal.In the case of decaying DM, the expected gamma-rayemission is written as follows:Φ( E i , n j ) = 14 π m χ τ Z ∆Ω j d Ω n Z . . dz cH ( z ) ρ χ ( z, n j ) Z E i +1 E i dE γ dN decay γ ( E γ (1 + z )) dE exp[ − τ EBL ( E γ (1 + z ))] . (16)Contrary to Eq. 15, Eq. 16 depends linearly on the DMdensity and it features the DM decay lifetime τ , insteadof h σ ann v i . A. Extragalactic resolved main halos and subhalos(EG-MSII)
We label halos and subhalos as “resolved” if theyare present in the Millennium-II catalog [60] with amass larger than 6 . × M ⊙ /h . We employ thesame procedure used in Ref. [34] to fill the regionbelow z = 2 .
15 with copies of the original Millennium-II simulation box (see Refs. [29] and [34] for furtherdetails). This provides a possible realization of thedistribution of resolved extragalactic DM halos andsubhalos along the past light-cone. The sky map oftheir emission (i.e., what we call “EG-MSII” in thefollowing) is obtained by determining, for each pixel inthe map, which DM structures fall inside the angulararea of the pixel (completely or partially, according totheir size) and by summing together their gamma-rayflux. In the case of annihilating DM, the annihilationrate of a DM halo or subhalo is computed from V max and r max (i.e., the maximal circular velocity and the distancefrom the center of the halo where this occurs) and byassuming that all DM structures are characterized bya Navarro-Frenk-White (NFW) density profile [65]. Adifferent choice of density profile would affect the overallintensity of the DM-induced emission (by a factor aslarge as 10, between extreme cases such as the Moore Refs. [29, 34, 63] show that more than 90% of the emission isproduced below z = 2. N -body simulationswere performed assuming cosmological parametersfavored by WMAP 1. Adopting the most recentvalues in agreement with the Planck mission [68] couldmodify the clustering and abundance of DM structuresin the simulations. However, it was shown that theincreased matter density Ω m and the decreased linearfluctuation amplitude σ (with respect to WMAP 1) havecompensating effects [69] and, therefore, we neglect thedependence of our results on the cosmological parameters(see also Ref. [70]).As in Refs. [29, 34], the way the copies of theMillennium-II simulation box are positioned around theobserver is a random process. Changing their orientationmodifies the distribution of resolved DM halos andsubhalos, affecting the shape of the auto- and cross-APSfor EG-MSII. Ref. [34] showed that this is just a 10%effect that can be neglected in comparison with othersources of uncertainty that will be mentioned later.For EG-MSII, the photon yield dN ann γ /dE includes theprimary gamma-ray emission (taken from Ref. [71]), i.e.hadronization of particles produced in the annihilation,final state radiation and internal bremsstrahlung. Wealso consider secondary emission, namely the photonsup-scattered by the IC of DM-induced electrons ontothe cosmic microwave background (see Ref. [34] fordetails). In the case of decaying DM, the photon yield isdetermined as dN decay i ( E ) /dE = dN ann i (2 E ) /dE , where i stands for either photons or electrons. B. Extragalactic unresolved main halos(EG-UNRESMain)
The emission of unresolved main halos (i.e. with amass smaller than 6 . × M ⊙ /h ) all the way downto the mass of the smallest self-bound halos M min , isreferred to as “EG-UNRESMain”. M min depends onthe nature of the DM particle and on its interactionswith normal matter but, at least within the contextof supersymmetric WIMPs, values between 10 − M ⊙ /h and 1 M ⊙ /h are reasonable, while M min = 10 − M ⊙ /h has become a popular benchmark [72, 73].In order to estimate EG-UNRESMain, we assumethat unresolved main halos share the same clusteringproperties of main halos with a mass between 1 . × M ⊙ /h and 6 . × M ⊙ /h . These main halos are justbelow our threshold of resolved DM structures. They arebarely resolved in the Millennium-II simulations (witha number of particles between 20 and 100) and theypopulate a regime in mass where the linear halo bias reaches a plateau [60]. Assuming that this remains trueeven below 1 . × M ⊙ /h , it is reasonable to think thatunresolved main halos will have a similar linear halo biasas the main halos with a mass between 1 . × M ⊙ /h and 6 . × M ⊙ /h . Therefore, their emission canbe accounted for simply by artificially enhancing theemission of the main halos between 1 . × M ⊙ /h and6 . × M ⊙ /h . Such an enhancement is implementedas follows: Z . × M ⊙ /hM min dM dn h ( M ) dM L ih ( M ) , (17)where dn h /dM is the main-halo mass function and L ih ( M ) is the gamma-ray flux produced by a singlemain halo with mass M . The index i stands for“ann” or “decay”, accordingly. The mass function isassumed to follow a power law in mass, down to M min .Its normalization and slope are fixed by fitting theabundance of main halos above the mass resolution ofMillennium-II, separately in the different snapshots of thesimulation, in order to reproduce the redshift dependenceof dn h /dM .Accounting for the contribution of main halos belowthe mass resolution of Millennium-II by enhancing theemission of the main halos with a mass between 1 . × M ⊙ /h and 6 . × M ⊙ /h is equivalent to assumingthat the two populations of DM halos share the samespatial distribution. Following the formalism introducedin Ref. [45] this is also equivalent to assuming that theyare characterized by the same 2-halo term.In the case of annihilating DM, the computation of L ann h in Eq. 17 is very sensivite to the concentration ofthe halo c ( M, z ). Contrary to Ref. [34], we consideronly the concentration model described in Ref. [74]:this model allows for c ( M, z ) to flatten as M decreases.Consequently, it agrees with the results of the recent N -body simulations in Refs. [75, 76] and is a more accuratemodel than the ones from Ref. [77, 78]. At z = 0and for M min = 10 − M ⊙ /h , Eq. 17 is 24 times largerthan the emission of all main DM halos with a massbetween 1 . × M ⊙ /h and 6 . × M ⊙ /h . For M min = 10 − M ⊙ /h ( M min = 1 M ⊙ /h ), the number is28 (15).For decaying DM, L decay h ( M, z ) = M and theenhancement (with respect to the emission of DM haloswith a mass between 1 . × M ⊙ /h and 6 . × M ⊙ /h ) is 6.7, 6.5 and 5.8 for M min = 10 − , 10 − and 1 M ⊙ /h , respectively. C. Extragalactic unresolved subhalos
In order to account for the emission of the subhalos ofunresolved halos, we modify Eq. 17 as follows: Z . × M ⊙ /hM min dM dn h ( M ) dM L ih ( M ) B i ( M, z ) . (18)5The additional term B i ( M, z ) is the so-called boostfactor, describing how much the emission of main halosincreases when the contribution of their subhalos isincluded. B decay is equal to 1 for all DM halosand redshifts, while the value of B ann ( M, z ) is quiteuncertain. We consider two scenarios that we believebracket the current uncertainty on B ann ( M, z ): • LOW scenario: this prescription is the same as inRef. [34] and it is motivated by the parametrization(performed in Ref. [79] and extended in Ref. [80]) ofthe probability P ( ρ, r ) of finding a value of the DMdensity between ρ and ρ + dρ in the data of the ViaLactea II N -body simulation. For this scenario andat z = 0, the overall enhancement in Eq. 18 (withrespect to the emission of DM main halos with amass between 1 . × M ⊙ /h and 6 . × M ⊙ /h )is 160, 88 and 23 for M min = 10 − , 10 − and 1 M ⊙ /h , respectively. • HIGH scenario: this is the same as the LOW recipebelow 10 M ⊙ /h while it predicts a boost factor5 times larger above that mass. Indeed, recentresults favor boost factors that are larger than theones of the LOW framework. Ref. [81] developeda semi-analytic model that accounts for the massaccretion rate of subhalos in larger host halos. Themodel also describes the effect of tidal strippingand dynamical friction experienced by subhalos.Including these effects increases the boost factorby a factor of 2-5, relative to the LOW scenario. Asimilar increase is expected when one accounts forthe fact that the concentration of subhalos changesaccording to the distance of the subhalo from thecenter of the host halo [82]. Finally, Refs. [83, 84]developed a new statistical method to describe thebehaviour of DM particles in collapsed structures,based on the modelling of the so-called ParticlePhase-Space Average Density (P SAD). Ref. [85]demonstrated that, when computed in the case ofDM subhalos, the P SAD is universal over subhalosof halos with a mass that goes from that of dwarfgalaxies to that of galaxy clusters. Employing areasonable parametrisation of the P SAD, Ref. [85]found boost factors that are as much as a factor of5 larger than the LOW case, at least for massiveDM halos.The boost factor predicted in Ref. [85] for DM haloswith a mass below 10 M ⊙ /h is, however, moderate(see their Fig. 5 ). Thus, when we compute theemission of unresolved main halos from 10 − M ⊙ /h to 10 M ⊙ /h as in Eq. 18, and we include the boost Note that the boost factor in Ref. [85] is defined in a differentway than in Ref. [34]. Whenever we take some information fromRef. [85], we translate it into the same definition used in Ref. [34]. factor of Ref. [85], we find a very similar result tothat of the LOW scenario. However, predictionsare different for massive DM halos: for objectswith a mass larger than 10 M ⊙ /h a boost factor5 times larger than for the LOW case is viable.We assume that such an increment is the same forall masses above 10 M ⊙ /h and we do not concernourselves with which mechanism (or combinationsof mechanisms) is responsible for it among theones mentioned above, since those studies agree onan increase of this magnitude. For the case with M min = 10 − M ⊙ /h ( M min = 1 M ⊙ /h ), the HIGHboost factor is defined to be a factor 8.4 (2.1) largerthan the LOW one (see Fig. 5 of Ref. [84] ).As in the case of resolved structures (EG-MSII), theemission of unresolved halos and subhalos is computedincluding the primary gamma-ray emission and thatresulting from the IC scattering off the cosmic microwavebackground. The total extragalactic signal is defined asthe sum of EG-MSII and EG-UNRESMain, boosted forthe emission of unresolved subhalos. We refer to the totalemission as “EG-LOW” and “EG-HIGH”, depending onthe subhalo boost factor scheme employed. D. The smooth halo of the Milky Way(GAL-MWsmooth)
As in Ref. [34], the emission of the smooth halo of theMW (called “GAL-MWsmooth”) is modeled by assumingthat the MW halo follows an Einasto profile [86]:log (cid:18) ρρ s (cid:19) = − α (cid:20)(cid:18) rr s (cid:19) α − (cid:21) . (19)The parameters in the above equation that provide thebest fit to the data of the highest resolution halo, Aq-A-1, in the Aquarius simulation [61, 62] are: ρ s =7 . × h M ⊙ / Mpc , r s = 11 .
05 kpc /h and α =0 . . × M ⊙ /h , defined as the amount of DM containedin a sphere with an average density of 200 times thecritical density of the Universe. Observationally, themass of the MW DM halo remains uncertain: Fig. 1of Ref. [87] shows how different methods (including, e.g.,MW mass modeling, dynamics of different tracers andthe study of the orbits of Andromeda and the MW)suggest values that go from 5 . × M ⊙ to 2 . × M ⊙ .Halo Aq-A-1 described above is on the higher end ofthis range. In order to account for the uncertaintyon the MW mass, we assume that ρ s can vary fromits nominal value of 7 . × h M ⊙ / Mpc down to Fig. 5 of Ref. [84] refers to the boost factor of MW-like DM halo.Assuming that similar results apply for all DM halos with a masslarger than 10 M ⊙ /h is therefore an approximation. . × h M ⊙ / Mpc . The latter corresponds to a MWmass that is 1/4 of the value of Aq-A-1 . Note thatthe intensity of the DM-induced gamma-ray emission isproportional to ρ s and to ρ s for an annihilating anddecaying DM candidate, respectively. On the otherhand, the auto- and cross-APS of GAL-MWsmooth areproportional to ρ s for annihilation-induced gamma raysand to ρ s for decaying DM. Thus, the uncertainty on theMW halo mass is a major systematic for the predictedsignal we are interested in.Ref. [61] also showed that a NFW profile provides areasonable fit to the Aq-A-1 data. We do not considerthis alternative here, since the difference compared to theEinasto profile described above would be evident onlybelow ∼ ◦ from the center of the MW, i.e. in a regionlocated inside the mask considered in the data analysis(see Sec. III B and Refs. [34] and [88]).In the case of GAL-MWsmooth, the photon yieldis computed including primary emission and secondaryemission from IC. The latter is computed using a fullmodel for the interstellar radiation field of the MW, asdescribed in Ref. [34], and not just the IC scatteringoff the cosmic microwave background. We also includethe hadronic emission, produced in the interactions ofDM-induced protons with the interstellar medium (seeRef. [34]). E. The subhalos of the Milky Way (GAL-AQ)
The last component to be considered is called GAL-AQ and it accounts for the emission of the subhalosof the DM halo of the MW. We derive this componentfrom the subhalo catalog produced in the Aquarius N -body simulation. Structures with a mass larger than1 . × M ⊙ are treated as “resolved”. We placethe observer at a distance of R = 8 . . × M ⊙ ), as they do not contributeto the auto- and cross-APS as argued in Refs.[34] and[89] .As in Ref. [34], only the primary gamma-ray emissionis considered when computing GAL-AQ. Since the observer is located approximately at a distance of R =8.5 kpc from the center of the MW, the best-fit Einastoprofile to Aq-A-1 corresponds to a local DM density of 0.45GeV/cm . An uncertainty of a factor 4 on ρ s would generatea variability of the same size on the local DM density. However, unresolved Galactic subhalos are expected tocontribute to the intensity of the DM-induced emission. Thisshould be kept in mind when, in Fig. 23, we compute ourpredictions for the DM-induced gamma-ray intensity.
Depending on the exact position of the observer onthe sphere with radius R and centered on the GalacticCenter, the distribution of resolved subhalos changes,and so does the intensity of GAL-AQ and its auto- andcross-APS. We estimate this variability by producing100 realizations of GAL-AQ, changing, each time, theposition of the observer on the sphere. We compute theauto- and cross-APS for each realization and note that,for annihilating DM, the 10% quantile of the distributionof the auto-APS (at ℓ = 400) among the 100 realizationsis a factor ∼ ∼ ℓ = 49 ( ℓ = 706). We suspect that thedistribution gets more peaked (i.e. less variable) at largemultipoles because it becomes more sensitive to the innerstructure of DM subhalos (which is constant among therealizations), instead of their distribution in the sky. Inthe case of decaying DM, the 10% and 90% quantiles arealways less than a factor 1.5 away from the median (seethe bands in the right panel of Fig. 20). The variabilityinduced by changing the position of the observer is animportant component in the total uncertainty of our DMpredictions and it will be considered in the followingsections.When discussing GAL-MWsmooth (Sec. VII D), weconsidered the effect of allowing the MW mass todecrease by a factor 4 with respect to the nominal valueof Aq-A-1. This has an impact also on GAL-AQ asthe number of subhalos in a DM structure is found tobe proportional to the mass of the host halo [90, 91].If we define k as the fraction by which we decreasethe MW mass, we consider 16 values of k , from 0.0to 0.25. For each k , we randomly remove a fraction k of the subhalos in the Aquarius catalog, to simulate alighter MW DM halo. For each value of k , we produce100 realizations of GAL-AQ for different positions ofthe observer. We compute the auto- and cross-APS foreach of the realizations. In Fig. 20, the lines show themedian of the distribution of the auto-APS of GAL-AQ(for a fixed multipole) as a function of k (and, thus,as a function of the MW mass). The left panel is forannihilating DM and the right one for decaying DM. Thesolid black line is for the auto-APS at ℓ = 400, while thelong-dashed red (short-dashed blue) one is for ℓ = 49( ℓ = 706). The coloured band (when present) denotesthe scatter between the 10% and the 90% quantiles in thedistribution among the 100 realizations. For annihilatingDM (left panel), the band becomes larger as the mapis populated by less and less subhalos and, therefore,it depends more and more on their distribution. Thevariability induced by our partial knowledge of the MW Here and for all the sky maps simulating the DM-inducedemission, the auto- and cross-APS are computed on the maskedgamma-ray sky with the anafast routine of
HEALPix afterhaving subtracted the monopole and dipole contributions. ] [M M400 600 800 1000 1200 1400 1600 1800 × s r ] - s r - s - APS o f G A L - A Q a t a f i x ed m u l t i po l e [ c m -28 -27 -26 -28 -27 -26 l=400 (median)10% and 90% quantiles (l=400)l=49 (median)l=706 (median) Annihilation ] [M
M400 600 800 1000 1200 1400 1600 1800 × s r ] - s r - s - APS o f G A L - A Q a t a f i x ed m u l t i po l e [ c m -27 -26 -25 -24 -27 -26 -25 -24 l=400 (median)10% and 90% quantiles (l=400)l=49 (median)10% and 90% quantiles (l=49)l=706 (median)10% and 90% quantiles (l=706) Decay
FIG. 20. Dependence of the APS of the GAL-AQ component on the MW mass.
Left:
The lines show the auto-APS at a fixedmultipole ( ℓ = 49 for the long-dashed red line, ℓ = 400 for the solid black one and ℓ = 706 for the short-dashed blue one)as a function of the mass of the DM halo of the MW, in our simulation of GAL-AQ described in the text. The auto-APS iscomputed between 0.5 and 0.72 GeV, for m χ = 2 .
203 TeV with a thermal annihilation cross section h σ ann v i = 3 × − cm s − and for annihilations into b ¯ b . For each value of the MW mass, 100 realizations of GAL-AQ are computed for different positionsof the observer. The lines refer to the median of the distribution of the corresponding auto-APS, while the grey band denotesthe variability between the 10% and the 90% quantiles of the distribution. The band is present only for the case with ℓ = 400for clarity. Right : The same as in the left panel but for decaying DM. The auto-APS is computed in the same energy bin, forthe same m χ and a decay lifetime of 2 × s. mass is another important source of uncertainty that willbe considered in the following sections.For some values of the DM mass, annihilation crosssection and decay lifetime, the gamma-ray flux of someDM subhalos in GAL-AQ may exceed the Fermi
LATsource sensitivity threshold. These DM subhalos wouldappear as resolved sources in the sky and they would beincluded in the 3FGL catalog. Since the auto- and cross-APS are measured masking the sources in 3FGL, DMsubhalos that are bright enough to be detected should beneglected when simulating GAL-AQ. Being very bright,they may be responsible for a significant fraction of theauto- and cross-APS of GAL-AQ. Thus, neglecting themmay affect significantly our predictions for GAL-AQ, asnoticed in Ref. [33] . In order to test this, we define theso-called particle physics factors Φ annPP and Φ decayPP , whichgathers all the terms in Eqs. 15 and 16 that do not dependon the DM distribution. More precisely:Φ annPP = h σ ann v i m χ Z m χ E thr E dN ann γ dE dE (20) One should also check that none of the DM halos or subhalosin EG-MSII are bright enough to be detected individually. Wedo not perform such a test because, even if some DM structureswere to be removed, this would hardly affect the prediction forthe auto- and cross-APS of EG-LOW and EG-HIGH. and Φ decayPP = 1 m χ τ Z m χ / E thr E dN decay γ dE dE, (21)where we choose a reference energy E thr of 0.1 GeV.We consider a reasonable range for the particle physicsfactors that goes from 10 − to 10 − cm s − GeV − for Φ annPP and from 10 − and 10 − s − for Φ decayPP 23 .This is divided in 50 logarithmic bins and, for eachbin, we build 100 realizations of GAL-AQ, varying theposition of the observer. For each particle physics factorand for each realization, we identify the subhalos (ifany) with an energy flux above 0.1 GeV that is largerthan the sensitivity flux in 3FGL at | b | > ◦ , i.e.3 × − erg cm − s − [92]. We consider the energyflux and not the number flux, since the Fermi
LATsensitivity, expressed in terms of the energy flux, ismore independent of the shape of the gamma-ray energyspectrum than when it is expressed by the number flux.Also, we use the sensitivity obtained for point sources,noting that the majority of the emission in a DM clumpcomes from a region within its scale radius r s . In a typicalrealization of GAL-AQ, almost all resolved DM subhaloshave a r s that corresponds to an angular size smaller For a DM mass of 200 GeV and annihilations (decays) into b ¯ b ,the range mentioned above corresponds to a variation between3 . × − cm s − and 3 . × − cm s − for h σ ann v i (between1 . × s and 1 . × s for τ ). Fermi
LAT at the energies of interest here.Thus, DM subhalos in GAL-AQ are rarely extended andthe use of the point-source sensitivity is well motivated.In Fig. 21, the solid lines show the median (over the 100realizations) of the number of subhalos that have beenexcluded because they are too bright, as a function ofthe annihilation (red line, bottom axis) and decay (blueline, top axis) particle physics factor. At the upper endof the range considered for Φ PP , this correction affectsbetween 500 and 2000 DM subhalos for annihilating anddecaying DM, respectively. These numbers correspondto approximately 1-2% of the total number of subhalosconsidered in the Aquarius catalog . The colored bandsindicate the variability associated with the 10% and 90%quantiles among the 100 realizations. The dashed verticallines are included as a reference and they correspondto the particle physics factor for an annihilating DMcandidate with a mass of 200 GeV and h σ ann v i =10 − cm s − (dashed red line) and for a decaying DMcandidate with the same mass and τ = 2 × s (dashedblue line). In both cases annihilations/decays into b ¯ b are considered and the values chosen for h σ ann v i and τ correspond approximately to their exclusion limits (for m χ = 200 GeV and for the REF benchmark scenario, seelater) as they will be computed in the following sections.This tells us that, for a given DM mass, the allowedregion in the parameter space of decaying DM wouldhave almost no DM subhalos that are too bright. Onthe other hand, the impact of bright subhalos may beimportant in the case of annihilating DM and this effectwill be considered when deriving the exclusion limits on h σ ann v i .In Fig. 22 we see the effect on the auto-APS ofneglecting the DM subhalos in Aquarius that are toobright. The left (right) panel shows the auto-APS ata specific multipole, for an annihilating (decaying) DMcandidate as a function of Φ annPP (Φ decayPP ). The auto-APSis multiplied by (Φ annPP ) − and by (Φ decayPP ) − , respectively,so that deviations with respect to a horizontal lineindicate how much the auto-APS is suppressed due tothe excluded subhalos. Note that the default maskcovering the sources in 3FGL is used when computingthe auto-APS. The solid black line is for ℓ = 400, whilethe red and blue ones are for ℓ = 49 and ℓ = 706.They indicate the median over the 100 realizations, whilethe grey band (sometimes difficult to see because it istoo narrow) represents the variablity between the 10%and 90% quantiles. As we anticipated in Fig. 21, the Note that only 1010 unidentified sources are present in 3FGL [38].Therefore, a particle physics factor that yields more than 1010DM subhalos with a flux larger than the 3FGL source sensitivityshould be excluded. As we will see in the following sections, theregion in the parameter space of DM that is not excluded bythe measured auto- and cross-APS does not correspond to thoseextreme values of particle physics factor. ] -1 GeV -1 s [cm annPP Φ -30 -29 -28 -27 -26 -25 N u m be r o f e xc l uded s ubha l o s ] -1 [s decayPP Φ -30 -29 -28 -27 -26 -25 -24 Annihilation10% and 90% quantilesDecay10% and 90% quantiles -1 s cm -24 v)=10 ann σ =200 GeV and ( χ for m annPP Φ s × =2 τ =200 GeV and χ for m decayPP Φ FIG. 21. The solid lines show the number of the DM subhalosin GAL-AQ with an energy flux above 0.1 GeV that is largerthan 3 × − erg cm − s − , i.e., the point-source sensitivityof Fermi
LAT in the 3FGL catalog [92]. The solid lines denotethe median over 100 independent realizations differing by theposition of the observer and they are plotted as a function ofthe particle physics factor, in the case of an annihilating DMcandidate (red line and bottom axis) and for a decaying one(blue line and top axis). In both cases, annihilations/decaysinto b ¯ b are considered. The colored bands indicate the 10%and 90% quantiles among the 100 realizations. For reference,the Φ annPP for m χ = 200 GeV, h σ ann v i = 10 − cm s − andannihilation into b ¯ b is marked by the dashed red line. Finally,the dashed blue line corresponds to the Φ decayPP for m χ = 200GeV, τ = 2 × s and decaying into b ¯ b . effect of neglecting bright subhalos starts to be importantaround 10 − cm s − GeV − in the left panel and around10 − s − in the right panel. The same values of theParticle Physics factor marked by the vertical lines inFig. 21 are plotted in Fig. 22 by the solid grey lines.The effect of DM subhalos that are too bright isaccounted for by defining the following quantity: κ (Φ PP , ℓ ) = C ℓ (Φ PP ) C ℓ (Φ minPP ) , (22)where Φ minPP = 10 − cm s − GeV − for annihilating DMand Φ minPP = 10 − s − for decaying DM. κ is computedby using the median over the 100 realizations. Wewill employ it as a correction factor to account for thebright DM subhalos that should be masked and it willbe multiplied by the APS of GAL-AQ with all the DMsubhalos. F. Results
In this section we define some benchmark cases thatwe will use in the following to discuss our main results:9 ] -1 GeV -1 s [cm annPP Φ -30 -29 -28 -27 -26 -25 ] - s r G e V - [ c m - ) ann PP Φ ( l C l=400 (median)10% and 90% quantiles (l=400)l=49 (median)l=706 (median) -1 s cm -24 v)=10 ann σ =200 GeV and ( χ for m annPP Φ Annihilation ] -1 [s decayPP Φ -30 -29 -28 -27 -26 -25 -24 s r ] - s r - [ c m - ) de c a y PP Φ ( l C l=400 (median)10% and 90% quantiles (l=400)l=49 (median)l=706 (median) s × =2 τ =200 GeV and χ for m decayPP Φ Decay
FIG. 22. Dependence of the APS of the GAL-AQ component on the Particle Physics factor defined in the text.
Left:
The solidlines show the DM-induced APS (computed above 0.1 GeV, for a fixed multipole and multiplied by (Φ annPP ) − ) as a functionof Φ annPP , neglecting the DM subhalos that would be detected individually according to the Fermi
LAT sensitivity threshold inthe 3FGL catalog [92]. The black, red and blue lines are for ℓ = 400, ℓ = 49 and ℓ = 706. They indicate the median overthe 100 realizations with different positions for the observer, while the grey band (only for the case with ℓ = 400) shows thevariabilty between the 10% and 90% quantiles. For reference, the value of Φ annPP for m χ = 200 GeV, h σ ann v i = 10 − cm s − andannihilation into b ¯ b is marked by the grey vertical line. Right:
The same as in the left panel, but for a decaying DM candidate.The vertical grey line is the particle physics factor of a DM candidate with m χ = 200 GeV, τ = 2 × s and decaying into b ¯ b .Note that the default mask covering 3FGL sources is employed when computing the auto-APS. • REF: this is our reference case and it is constructedby summing EG-MSII and EG-LOW, with M min =10 − M ⊙ . We also include GAL-MWsmooth (forthe nominal value of the MW DM halo, taken fromAq-A-1, of 1 . × M ⊙ /h ) and the median ofGAL-AQ over the 100 realizations produced for thenominal MW DM halo mass; • MAX: we build this case by maximizing allthe uncertainties considered (and discussed inthe previous sections). Thus, we take it as agood estimate of the largest signal that can beassociated with DM (for a given value of theparticle physics factors and of M min ). The MAXbenchmark is defined by summing EG-MSII, EG-HIGH (for M min = 10 − M ⊙ ), GAL-MWsmooth(for a nominal mass of the MW DM halo) and the90% quantile among the 100 realizations of GAL-AQ relative to a 1 . × M ⊙ /h MW; • MIN: contrary to MAX, this benchmark is obtainedby tuning all the uncertainties considered aboveto their minimal configuration. In particular, wesum EG-MSII, EG-LOW (for M min = 10 − M ⊙ ),GAL-MWsmooth (for a MW mass that is 1/4 ofthe nominal value of Aq-A-1) and the 10% quantileof the 100 realizations of GAL-AQ for a MW massthat is 1/4 of the value of Aq-A-1. In order to discuss the effect of changing M min , we alsocompute the MIN and MAX benchmarks for M min =10 − M ⊙ and M min = 1 M ⊙ .Fig. 23 shows our predictions for the intensity of theDM-induced emission, averaged over the whole sky. Theleft panel is for annihilating DM with a mass of 212GeV and h σ ann v i = 3 × − cm s − , while the rightpanel is for a decaying candidate with the same mass and τ = 2 × s. Annihilations and decays produce gammarays via b ¯ b . The solid black line is for the REF scenario,while the red and blue ones are for the MIN and MAXbenchmark (for M min = 10 − M ⊙ ). Thus, the grey bandbetween the red and blue lines indicates how much ourpredictions change when accounting for the uncertaintiesmentioned above.For annihilating DM, in the case of the REFbenchmark, the emission is contributed, almost equally,by EG-LOW and GAL-MWsmooth. Thus, the differencebetween REF and MAX comes entirely from thedifferent subhalo boost employed to describe unresolvedextragalactic DM structures (see Sec. VII C). The boostfactor is larger at higher redshifts and, therefore, atenergies close to m χ , where the emission is dominatedby nearby sources, the red line approaches the blackone. On the other hand, in the LOW scenario thecontribution of GAL-MWsmooth is suppressed and,therefore, the intensity of the LOW benchmark is almosthalved compared to REF.Subhalo boosts do not affect the predictions fordecaying DM and, therefore, the REF and the MAX0 Energy [GeV]1 10 G e V ] - s r - s - ( E ) / d E [ c m Φ d E -10 -9 -8 -7 -10 -9 -8 -7 REF scenarioMIN scenarioMAX scenario (MIN) min
Uncertainty on M (MAX) min
Uncertainty on M
Annihilation
Energy [GeV]1 10 G e V ] - s r - s - ( E ) / d E [ c m Φ d E -10 -9 -8 -7 -10 -9 -8 -7 REF scenarioMIN scenarioMAX scenario
Decay
FIG. 23.
Left:
The predicted energy spectrum of the gamma-ray emission induced by the annihilation of a DM particle witha mass of 212 GeV, h σ ann v i = 3 × − cm s − and annihilation into b ¯ b . The black line stands for the REF scenario, whilethe red and blue ones are for the MAX and MIN cases. Thus, the grey band determines the variability between the MIN andMAX scenarios. In all cases M min = 10 − M ⊙ (see text for details). The red and blue shaded areas indicate how the MAX andMIN benchmarks change if we let M min vary between 10 − and 1 M ⊙ . Right : The same as in the left panel but for a decayingDM particle with a mass of 212 GeV and a decaying lifetime of 2 × s. The red and black lines overlap. benchmarks in the right panel overlap . The lowerintensity of the LOW case is, as before, due to thesuppression of GAL-MWsmooth.In the left panel, the blue and red shaded areas indicatehow our predictions for MIN and MAX change whenallowing M min to vary in the range mentioned above.These uncertainty bands get larger for smaller energies,as the signal becomes sensitive to the emission at higherredshifts. Changing the minimal DM halo mass has avery minor effect on decaying DM and, thus, the shadedbands are not present.The predictions of Fig. 23 can be compared withFig. 12 of Ref. [34]: the main difference is the factthat our brightest configuration (i.e. the MAX scenario)predicts almost one order of magnitude less gamma-rayflux than in Ref. [34], given the new definition of theHIGH subhalo boost factor. On the other hand, ourpredictions are compatible with the results of Ref. [31].In Fig. 24 we show the predicted auto-APS in theenergy bin between 1.38 and 1.99 GeV, for the sameDM candidates considered in Fig. 23. Note that thesecorrespond to particle physics factors that are quite lowand therefore none of the DM subhalos in GAL-AQ wouldbe resolved by Fermi
LAT. We can then neglect the κ correction discussed in Sec. VII E. The left panel isfor annihilating DM and the right one for decaying DM.The predicted intensity APS is shown separately for thedifferent components discussed above. We also includethe Poissonian auto-APS measured by Fermi
LAT in The contribution of GAL-AQ is subdominant. the same energy bin (solid black line) and its estimatederror (grey band), for the mask around 3FGL sources.To compare the measurement with the predicted DMsignal, the latter needs to be corrected for the presenceof the mask. In the analysis of the
Fermi
LAT data,this is done automatically by
PolSpice (see Sec. III B),while we correct our predictions by dividing the auto-and cross-APS obtained from the masked sky by f sky ,i.e., the fraction of the sky outside the default maskdefined in Sec. III B. Such a recipe is based on theassumption that the masked region is characterised bythe same clustering properties of the unmasked one. Wetest this hypothesis by computing the auto-APS of oursimulated sky maps with and without the mask. For theextragalactic signal we find that, indeed, dividing themasked auto-APS by f sky , we reproduce the unmaskedone. On the other hand, for GAL-AQ, the maskedauto-APS is 0.43 times smaller than the unmasked one(approximately at all multipoles). This factor is largerthan f sky , which suggests that the distribution of DMsubhalos outside the mask is slightly more isotropicthan the distribution inside the mask. This is to beexpected since DM subhalos are more clustered towardsthe center of the host halo. Finally, the auto-APS ofGAL-MWsmooth is significantly different if we includethe mask: the intensity of the auto-APS decreases as weare covering the region which produces the majority ofthe emission. Also, the morphology of the mask inducessome spurious fluctuations on the auto-APS. A moresophisticated algorithm should be employed to correctfor these features. Alternatively, the wiggles would bereduced by smoothing the edges of the mask. However,1 Multipoles100 200 300 400 500 600 700 s r ] - s r - s - [ c m sky /f l C -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 Fermi-LAT measurement (1.38-1.99 GeV)GAL-MW smoothUncertainty on MW massGAL-AQ (median)GAL-AQ (10% and 90% quantiles)Adding uncertainty on MW massEG-LOWEG-HIGHUncertainty on boost factor
Annihilation M -6 =10 min , Mb, b -1 s cm -26 × v>=3 ann σ =212 GeV, < χ m Multipoles100 200 300 400 500 600 700 s r ] - s r - s - [ c m sky /f l C -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 Fermi-LAT measurement (1.38-1.99 GeV)GAL-MW smoothUncertainty on MW massGAL-AQ (median)GAL-AQ (10% and 90% quantiles)Adding uncertainty on MW massEG-LOWEG-HIGH
Decay M -6 =10 min , Mb, b -1 s × =2 τ =212 GeV, χ m FIG. 24.
Left:
Auto-APS in the energy bin between 1.38 and 1.99 GeV, for a DM candidate with a mass of 212 GeV, h σ ann v i = 3 × − cm s − and annihilation into b ¯ b . The auto-APS is divided by f sky to correct for the presence of the maskdescribed in Sec. III B. The black solid line and the grey band indicate the Poissonian auto-APS measured in this energy binand for the mask around 3FGL sources (see Sec. V). The solid blue line is the median of the auto-APS for GAL-AQ over the100 realizations with different positions for the observer and the blue band shows the variability between the 10% and 90%quantiles. The uncertainty band on GAL-AQ extends downwards (shaded blue area) if we account for an uncertainty of a factor4 in the value of the mass of MW DM halo. The red and purple lines show the auto-APS for EG-LOW and EG-HIGH, for M min = 10 − M ⊙ . The green line stands for the GAL-MWsmooth component and the green band accounts for an uncertaintyof a factor 4 in the mass of MW DM halo. The wiggles in this component are due to the mask applied to cover the Galacticplane (see text for details). Right : The same as in the left panel but for a decaying DM particle with a mass of 212 GeV and adecaying lifetime of 2 × s. The red and and purple lines overlap. we note that, in the signal region defined in Sec. III,the GAL-MWsmooth component is not responsible forthe majority of the signal and, therefore, using thereconstructed auto-APS of GAL-MWsmooth would notconsiderably change our results. Therefore, for both theGAL-AQ and the GAL-MWsmooth, we simply apply the1 /f sky correction.In the left panel of Fig. 24 we note that the dominantcontribution is GAL-AQ: the solid blue line is the medianover the 100 realizations with different observers (in thecase of a MW DM halo with the same mass as Aq-A-1), while the filled blue band shows the variabiltybetween the 10% and 90% quantiles. If we also includethe possibility that the MW DM halo may be up to 4times lighter (see Sec. VII E and Fig. 20), the uncertaintyband extends downwards to include the shaded blueband. Over the signal region, GAL-AQ is not constantand it decreases by approximately a factor of 10. Theextragalactic signal is plotted in red and purple, fora LOW and HIGH subhalo boost, respectively. Thisuncertainty gives rise to the pink band that coversapproximately one order of magnitude. The extragalacticcomponent becomes nearly constant for ℓ > ∼
300 but, overthe whole signal region, it decreases by a factor of 10.Finally, the GAL-MWsmooth is plotted in green and thegreen band indicates how much the signal decreases whenthe mass of MW DM halo is allowed to decrease by upto a factor of 4 with respect to the value of Aq-A-1. If we had considered M min = 10 − M ⊙ /h instead,the intensity of EG-LOW and EG-HIGH would havebeen approximately 4 times larger, while it would havedecreased by a factor 50 if we had considered M min =1 M ⊙ /h . However, since the EG-LOW and EG-HIGHare not dominant components, the effect of changing M min on the total DM signal will not be as large.In the right panel we follow the same color coding: themain difference with respect to the case of annihilatingDM is the fact that the extragalactic contributiondominates the signal for most of the measured signalregion. There is no uncertainty associated with the boostfactor and, therefore, the red and purple lines coincide.As in the left panel, the auto-APS is nearly constantfor ℓ > ∼
300 and it decreases by a factor ∼
50 overall.Another important difference, with respect to the case ofDM annihilation, is the fact that the auto-APS of GAL-AQ is much steeper, decreasing by a factor ∼
600 from ℓ = 49 to ℓ = 716.Independently of how the different components aresummed together , producing the different REF, MINand MAX scenarios described above, the total signal In principle, one should include the cross-correlation betweenthe different components considered. We do not expect anycorrelation between extragalactic and Galactic emission. Thecross-correlation between GAL-MWsmooth and GAL-AQ was not
Poissonian but decreasesat smaller angular scales. This will be crucial whencomparing our predictions to the
Fermi
LAT data.
VIII. USING THE AUTO- AND CROSS-APS TOCONSTRAIN DARK MATTER ANNIHILATIONAND DECAY
In this section we compare the predictions for theDM-induced auto- and cross-APS obtained in Sec. VIIwith the updated
Fermi
LAT measurement of theIGRB auto- and cross-APS presented in Sec. V. Sucha comparison will allow us to determine whether thedata are consistent with a DM interpretation or how wecan use them to constrain the nature of DM. We followtwo complementary approaches that will be describedseparately in the following subsections. Neither methodfinds a significant detection of DM in the auto- andcross-APS data and, therefore, the measurement is usedto derive exclusion limits on the intensity of the DM-induced gamma-ray emission, as a function of m χ . A. Conservative exclusion limits
This first strategy is motivated by the desire to beconservative. In particular, the DM-induced APS, forany energy bin or combination of bins, must not exceedthe measured data. For a certain benchmark case (amongREF, MIN and MAX) and for a certain value of M min ,imposing that constraint will translate into upper limitson the intensity of the DM-induced signal or, by fixing m χ and the annihilation/decay channel, into upper limitson h σ ann v i for annihilating DM and into lower limits on τ for decaying DM. We refer to these exclusion limits as“conservative”.We consider 60 values of m χ between 5 GeV and 5TeV. For each value of m χ , we compute the DM-inducedauto- and cross-APS in the 13 energy bins defined inSec. II, for the 3 benchmarks described in Sec. VII F,for 3 annihilation/decay channels (i.e., b ¯ b , τ + τ − and µ + µ − ) and for 3 values of M min (i.e., 10 − , 10 − and1 M ⊙ ). The APS associated with DM for energy bins i and j is averaged over the signal region in multipoleand we require it to be smaller than the Poissonian APSmeasured for that pair of energy bins plus 1.64 timesits error: h C i,jℓ, DM i < C i,j P + 1 . σ C i,j P . Assuming thatthe measured C P has a Gaussian probability distributionwith a central value of C P and a standard deviationof σ C P , values further away than 1.64 times σ C P from computed in Ref. [93] and it is at least one order of magnitudebelow the auto-correlation of GAL-AQ. We neglect it here. See Ref. [34] on how to compute the emission for multipleannihilation/decay channels without having to recompute, foreach case, the mock sky maps from the N -body simulations. the central value correspond to a cumulative probabilitydistribution larger than 0.95. Thus, excluding themprovides a 95% CL exclusion bound. For each m χ , wetake the most stringent limit among all the combinationsof energy bins.Fig. 25 shows the upper limits on the h σ ann v i (leftpanel) and the lower limits on τ (right panel), as afunction of m χ , for annihilations/decays into b ¯ b and forthe different benchmark scenarios considered above. Theblack line is for REF while the blue and red ones arefor MIN and MAX. Thus, the grey band represents ourtotal systematic astrophysical uncertainty (for M min =10 − M ⊙ ) and it is as large as approximately a factor of5 or 2, for annihilating and decaying DM, respectively.The limits have some wiggles because, depending on theDM mass, the emission peaks at different energies, anddifferent combinations of energy bins are responsible forthe exclusion limit. Solid lines are obtained consideringall the possible combinations of energy bins, while forthe black, blue and red dashed ones only the auto-APSis employed. The figure shows that, at large DM massand both for annihilating and decaying DM, the exclusionlimits are driven by the cross-APS and not by the auto-APS, approximately for m χ >
200 GeV for annihilatingDM and for m χ >
700 GeV for decaying DM.In the left panel, the red and blue shaded areas accountfor the effect of changing M min between 10 − M ⊙ and1 M ⊙ and they are computed only for the MIN andMAX scenarios. The effect is more important for theMIN case since the emission from extragalactic DMstructures contributes more to the total signal in thiscase (see Fig. 24). If we include the variability on M min in our budget for the systematic uncertainty, thesystematic error grows to a factor of 40. Compared tothe conservative upper limits on h σ ann v i derived by theintensity of the IGRB in Ref. [31], our uncertainty isapproximately a factor of 2 larger. The long-dashed greyline is the thermal annihilation cross section computedin Ref. [94]. The line marks the beginning of the regionwhere, for WIMP DM, one can find annihilation crosssections that correspond to a relic DM abundance inagreement with the Planck data [68]. Unfortunately, ourconservative upper limits do not probe this region, asthey are, at least, a factor 3 away from it. Also, the REFupper limit is approximately two orders of magnitudehigher than the upper limit derived from the observationof 15 dwarf spheroidal galaxies performed by Fermi
LAT[95] and included here as a grey dash-dotted line. Finally,it is a factor 2 higher than the conservative limits thatcan be derived from the intensity of the IGRB (short-dashed grey line) [31], at least below 100 GeV. Abovethis value, the IGRB intensity leads to an even morestringent exclusion.In the right panel, it can be seen that the lower limitson τ derived here from the auto- and cross-APS are afactor of 5 less stringent than the lower limits obtainedin Fig. 6 of Ref. [96] from the IGRB intenstity, in theirconservative scenario. The REF scenario is also at least3 [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)Fermi LAT IGRB intensity (Ackermann et al., conservative, 2015)REF scenarioMAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 bAnnihilation, b [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, conservative, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioMAX scenarioMIN scenario bDecay, b FIG. 25. Conservative exclusion limits on annihilating and decaying DM from the new APS measurement.
Left:
The solidlines show the upper limits on h σ ann v i derived from the auto- and cross-APS measured in Sec. III, as a function of m χ , for M min = 10 − M ⊙ and annihilations into b ¯ b . The limits follow the conservative approach described in the text. The black lineis for the REF scenario, while the red and blue ones are for MAX and MIN, respectively. The grey band between the MINand MAX scenario represents our estimated total astrophysical uncertainty for M min = 10 − M ⊙ , accounting for all the sourcesof uncertainty mentioned in Sec. VII. The red and blue shaded bands describe the effect of changing M min between 10 − M ⊙ and 1 M ⊙ , for the MAX and MIN scenario. In the case of the black, red and blue dashed lines, the upper limits are derivedonly considering the measured auto-APS and neglecting the cross-APS. For comparison, the long-dashed grey line marks theannihilation cross section for thermal relics from Ref. [94] and the dash-dotted grey line the upper limit obtained in Ref. [95]from the combined analysis of 15 dwarf spheroidal galaxies. Finally, the short-dashed grey line shows the conservative upperlimit derived in Ref. [31] from the intensity of the IGRB. Right : The same as in the left panel but for the lower limits on τ fordecaying DM. The short-dashed grey line represents the lower limit obtained in Fig. 6 of Ref. [96] from the IGRB intensity,while the dash-dotted grey one is obtained from the combined analysis of 15 dwarf spheroidal galaxies in Ref. [97]. a factor of 5 from the dash-dotted grey line showing thelower limits obtained from the combined analysis of 15dwarf spheroidal galaxies in Ref. [97].Figs. 35 and 36 in Appendix E show the exclusionlimits on h σ ann v i and τ in the case of the τ - and µ -channels. B. Fit to the data and realistic exclusion limits
In this section we describe our analysis of the auto-and cross-APS using a 2-component model that includesa Poissonian term and a DM-induced one which, as wenoticed in Fig. 24, deviates from a Poissonian behaviour.The Poissonian component is interpreted as the APS ofunresolved astrophysical sources, even if we do not try topredict its amplitude in terms of a specific model. This2-component model will be used to fit the
Fermi
LATAPS as a function of multipole.The fit minimizes the χ defined as: χ = X i,j,ℓ [ C ℓi,j − C i,jℓ, DM − C i,j P , astro ] [ σ ℓi,j ] , (23)where the i, j indexes in the sum extend over all the 91independent combinations of energy bins and the ℓ index runs over the 10 bins in multipoles contained in the signalregion. C ℓi,j indicates the APS measured in the ( i, j )combination of energy bins and in the ℓ multipole bin,while C i,jℓ, DM and C i,j P , astro are the DM and Poissoniancomponents of our model in the same combination ofenergy bins and in the same multipole bin. Finally, σ ℓi,j is the experimental error associated to C ℓi,j and providedby PolSpice . The DM APS C i,jℓ, DM are computed for thesame 60 values of m χ as in the previous section, the same3 annihilation/decay channels, 3 benchmark scenariosand 3 values of M min . The only remaining parameterneeded to calculate C i,jℓ, DM is either h σ ann v i or τ : theywill be fixed to a specific value every time we compute χ . On the other hand, the 91 independent values of C i,j P , astro in Eq. 23 are left free in the fit. Putting theDM term to zero in Eq. 23 defines our null hypothesis.In that case, the fit to the Fermi
LAT data leads usto the C P estimators discussed in Sec. III, whose auto-APS is plotted in Fig. 8. Including the DM component,for a fixed m χ , annihilation/decay channel, benchmarkscenario and M min , we repeat the minimization of χ inEq. 23, for different values of h σ ann v i and τ .We show an example in Fig. 26, for the case of theREF scenario for a DM candidate with a mass of 768.1GeV, h σ ann v i = 6 . × − cm s − annihilating into b ¯ b .The value of the annihilation cross section corresponds4 Multipole s r ] - s r - s - [ c m l C -50-40-30-20-1001020304050 -21 × Measured auto-APSBest-fit solution for the 2-components modelPoissonian fit (null hypothesis) best-fit for the 2-components model
P,astro C -50-40-30-20-1001020304050 -21 × Energy bin [10.45-21.83] GeV M -6 =10 min REF scenario, M b, b -1 s cm -24 × v)=6.1 ann σ =768.1 GeV, ( χ m FIG. 26. Example of a fit to the binned APS C ℓ in ourparticular energy bins, in terms of the 2-component modeldescribed in the text. The red circles show the measuredauto-APS in the energy bin between 10.4 and 21.8 GeV, asa function of multipole. The solid red line is the Poissonianbest-fit APS in the null hypothesis and the pink band denotesits estimated 68% CL error. The dashed blue line denotes thebest-fit of the Poissonian component when DM is includedin the fit, for a DM mass of 768.1 GeV and a h σ ann v i of6 . × − cm s − , annihilation into b ¯ b and a REF scenariowith M min = 10 − M ⊙ . The best-fit signal (Poissonian plusDM component) is plotted by means of the blue triangles. approximately to the exclusion upper limit for that valueof DM mass, as will be computed later. The red circlesshow the measured auto-APS as a function of ℓ inthe signal region for one reference energy bin, i.e., theone between 10.4 and 21.8 GeV (when masking 3FGLsources). The solid red line with the pink band denotesthe best-fit C P in that energy bin for the null hypothesis(i.e., without DM), while the dashed blue line is thebest-fit Poissonian component C P , astro when the fit isdone with the 2-component model (i.e. including DM).The dashed line is lower than the solid one, since atthese energies part of the signal is explained by DM and,therefore, there is less need of a Poissonian component.Energy bins not localized near the peak of the DMemission are only slightly affected by the inclusion of theDM term in the fit. The best-fit configuration for thetwo-components model is plotted by blue triangles: theinclusion of the DM term makes it multipole-dependentso that it decreases by a factor of ∼ χ ofthe best-fit point with respect to the null hypothesis, atleast for DM masses above few hundreds of GeV. This isprobably due to the fact that the measured auto-APSis slightly multipole-dependent. We can quantify theimprovement in the fit provided by the DM component by building a 2-dimensional grid in ( m χ , h σ ann v i ) forannihilating DM and in ( m χ , τ ) for decaying DM andplotting the TS ∆ χ , i.e. the difference between the χ ofthe best fit for the null hypothesis and the χ of the bestfit in the case of the 2-component model. This is shownin Fig. 27, where the left panel refers to annihilating DMand the right one to decaying DM (for the b channeland a REF scenario with M min = 10 − M ⊙ ). In bothpanels, the closed area indicate the region where the 2-component model is preferred over the null hypothesisat a 68% CL. The 90% and 95% CL regions are openand bounded by the corresponding white lines . Thistells us that including the DM component in the modelprovides a better fit to the auto- and cross-APS measuredin Sec. III, with a significance between 1 and 1.6 σ . Thisis too small to consider as significant. Thus, we concludethat the data do not significantly prefer the addition ofa DM component and we use the measured auto- andcross-APS to derive constraints on the DM signal.The contour plots for the τ - and the µ -channel can beseen in Appendix E. In both cases, the 68% CL region isthe only closed one.For each value of m χ , M min , annihilation/decaychannel and benchmark scenario, the exclusion limits on h σ ann v i and τ are derived by scanning on h σ ann v i and τ until we find the values that correspond to a best fitwith a TS ∆ χ of 3.84 with respect to the null hypothesis.Such a value is derived assuming that ∆ χ follows a χ probability distribution with one degree of freedom (i.e., h σ ann v i or τ ) and noting that values larger than 3.84 falloutside 5% of the cumulative distribution probability.This recipe provides the 95% CL exclusion limits on h σ ann v i and τ that are summarised in the left and rightpanels of Fig. 28, respectively.In the left panel, as in Fig. 25, the black, blue andred solid lines correspond to the REF, MIN and MAXscenarios. The difference between MIN and MAX coversslightly more than a factor of 5. The blue and redshaded regions around the solid lines of the same colorindicate how the upper limits change when we leave M min free to vary. This extends the range of the totalsystematic uncertainty to approximately a factor of 20.For comparison, the black dashed line is the REF upperlimit in its conservative version (from Fig. 25). Fittingthe data with the 2-component model generates exclusionlimits that are approximately a factor of 10 stronger,at least at low DM masses. As the mass increases,the method employed in this section starts to performprogressively worse and the solid black line gets closerto the dashed one. This is due to the fact that, for m χ >
150 GeV, the data slightly prefer the interpretationwith DM as opposed to the null hypothesis. The figure The 68% CL ares is obtained by identifying the region where thebest-fit solution for the 2-component model has a TS ∆ χ of 2.30larger than the null hypothesis. The values are 4.61 and 5.99 forthe 90% and 95% CL regions. mass [GeV]10 ] - s v > [ c m ann σ < -27 -26 -25 -24 -23 -22
10 -2-10123456
68% CL90% CL95% CL bAnnihilation, b mass [GeV]10 [ s ] τ
10 -2-10123456
68% CL90% CL95% CL bDecay, b
FIG. 27. ∆ χ between the best-fit solution for the 2-component scenario and the best fit in the null hypothesis. Resultspresented here refer to the REF scenario for annihilation/decay into b ¯ b and M min = 10 − M ⊙ . Left:
Each point in the( m χ , h σ ann v i ) parameter space is colored according to its ∆ χ , i.e. the difference between the χ of the best fit to the auto-and cross-APS in terms of the 2-component model and the χ of the best fit of the null hypothesis (i.e. no DM). The closedwhite contour marks the 68% CL region. The 90% and 95% CL regions are below the white open curves labelled “90% CL”and “95% CL” respectively. Right : The same as the left panel but for decaying DM. also includes the thermal cross section from Ref. [94] asa long-dashed grey line: our upper limit for the REFcase is slightly above it, below 10 GeV. It is also morethan a factor 10 weaker than the upper limit derivedfrom the observation of 15 dwarf spheroidal galaxies inRef. [95]. Finally, the short-dashed grey line indicatesthe exclusion limits obtained in Ref. [98] by studying theintensity of the IGRB with a 2-component model that,similarly to what is done here, includes both a genericmodel-independent astrophysical contribution and a DMone. Our REF limit is slightly stronger than the short-dashed grey line for m χ <
30 GeV, suggesting that thestudy of the IGRB anisotropies could in principle be amore effective way of constraining DM than the IGRBintensity. However, for larger DM masses, our limit getsworse due, again, to the fact that the data slightly preferan interpretation that includes DM.The same color coding is used in the right panel fordecaying DM. With no dependence on M min , the bandof the systematic uncertainty covers a factor 2, and theREF upper limit is even one order of magnitude above theconservative one, at least at 60-70 GeV. For larger massesour limit worsens for the same reason as in the left panel.As in Fig. 25, the short-dashed grey line is the lower limitobtained from the analysis of the IGRB intensity fromRef. [96]. The line refers to the case in which the IGRBis modeled in terms of a component with a power-lawenergy spectrum and a DM contribution. Above 20 GeV, where both lines are available, the analysis of the IGRBintensity is always more powerful than the anisotropystudy performed here. Finally, the dot-dashed grey lineis the lower limit obtained from the analysis of 15 dwarfspheroidal galaxies in Ref. [97]. Our REF scenario isalways below this line, at least by a factor of 2.In Appendix E we include the exclusion limits for the τ - and µ -channels.When fitting with the 2-component model, the 91 C i,j P , astro can vary independently and they react to thepresence of the DM component by reproducing themeasured APS in those combinations of energy binswhere the DM component is subdominant. Therefore,it may difficult to interpret a best-fit set of C i,j P , astro interms of one or more populations of actual astrophysicalsources, e.g. unresolved blazars or star-forming galaxies.A more physical approach can be obtained by consideringthe phenomenological description presented in Sec. VI.In this case, the astrophysical component in the 2-component fit is described by means of the two-broken-power-law scenario (see Tab. III). The latter depends on8 free parameters instead of 91. We employ this revisedversion of the 2-component model to fit the binned auto-and cross-APS C ℓ in all the combinations of energy bins.The exclusion limits on h σ ann v i and τ are obtained byfinding the configuration that yields a χ that is largerby 3.84 than the best-fit χ of the null hypothesis (i.e.no DM). The resulting upper limits are compatible with6 [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 -20 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)Fermi LAT IGRB intensity (Ajello et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 -20 bAnnihilation, b [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, power-law, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario bDecay, b FIG. 28. Exclusion limits on annihilating and decaying DM from the fit to the binned C ℓ in terms of the 2-component model. Left:
The solid lines show the upper limits that can be derived on h σ ann v i as a function of m χ (for annihilations into b ¯ b quarksand M min = 10 − M ⊙ ) by fitting the Fermi
LAT data with a 2-component model that includes astrophysical sources and DM(see text for details). The black, blue and red lines correspond to the REF, MIN and MAX scenario. The blue and red shadedareas indicate how the MIN and MAX upper limits change when leaving M min free to vary between 10 − M ⊙ and 1 M ⊙ . Theblack dashed line is the REF upper limit in the conservative case, from Fig. 25, while the long-dashed grey line is the thermalannihilation cross section from Ref. [94]. The dot-dashed grey line is the upper limits derived in Ref. [95] from the combinedanalysis of 15 dwarf spheroidals, while the short-dashed grey line comes from the analysis of the IGRB intensity performed inRef. [98]. Right : The same as in the left panel but for the lower limits on τ , in the case of decaying DM. The short-dashedgrey line represents the lower limit obtained in Ref. [96] from the IGRB intensity. The line is taken from Fig. 5 of Ref. [96],where the IGRB is interpreted in terms of a component with a power-law emission spectrum and a DM contribution. Finally,the dot-dashed grey line is the upper limit from the analysis of 15 dwarf spheroidal galaxies performed in Ref. [97]. the ones showed in Fig. 28 and, therefore, we decided notto show them. IX. DISCUSSION AND CONCLUSION
In this paper we measure the auto-correlation andcross-correlation angular power spectrum (APS) of thediffuse gamma-ray emission detected by
Fermi
LAT athigh Galactic latitudes in 81 months of observation. Themeasurement builds on a similar analysis based on 22months of data and published in Ref. [1]. With respectto the latter, this work takes advantage of the largerstatistics, as well as of the improved event reconstructionachieved for Pass 7 Reprocessed events and instrumentresponse functions. Other improvements, with respect toRef. [1], consist of a revised method for binning the datain multipole and to compute the Poissonian auto- andcross-APS. We also correct the estimate of the photonnoise and we employ a different method to account forthe effect of the mask. Finally, we consider a more recentmodel of the diffuse Galactic foreground associated withthe Milky Way (MW) disk.The second part of the paper focuses on the auto-and cross-APS expected from annihilation or decayof DM. We employ a hybrid approach to model thedistribution of DM, making use of catalogs of DM halos and subhalos from state-of-the-art N -body simulationsand combining them with analytical recipes to accountfor DM structures below the mass resolution of thesimulations. The methodology follows what was done inRef. [34]. Compared to the latter, this work discards thepossibility of very large subhalo boost factors inducedby na¨ıve power-law extrapolations of the concentrationparameter to low halo masses. We also account for theuncertainty associated with the mass of the MW, and wecorrect for the possibility of having very bright GalacticDM subhalos that would be individually resolved asgamma-ray sources.The main results of this papers are summarized in thefollowing list. • Detection of auto- and cross-APS: because of theinstrumental improvements and of the refinementsin the analysis mentioned before, the measurementpresented here probes a larger energy range(compared to the original analysis in Ref. [1]),between 0.5 and 500 GeV, divided in 13 energybins. We also compute, for the first time, thecross-APS between different energy bins. We detectsignificant auto-APS in almost all the energy binsbelow 21.83 GeV. Significant cross-APS is alsomeasured in most combinations of energy bins (seeTabs. I and II).7 • Independence on angular multipole: our resultscover a larger range in multipoles than the originalanalysis, i.e., from ℓ = 49 and 706. In thismultipole range, the detected auto- and cross-APS are consistent with being Poissonian, i.e.constant in multipole. An alternative ℓ -dependentmodel is also employed to fit the data but thereis no significant preference for the ℓ -dependentmodel over the Poissonian interpretation. If futuredata sets were able to detect a non-Poissonianbehaviour, it would represent the first detection ofscale dependence in gamma-ray anisotropies. Sucha result would provide valuable insight into thenature of the Isotropic Gamma-Ray Background(IGRB), e.g. an upper limit on the contributionof sources like blazars or misaligned active galacticnuclei, which are associated with a Poissonian APS.It would also probe other possible sources like star-forming galaxies or Dark Matter (DM) structures,from which we expect a ℓ -dependent auto- andcross-APS [8]. • Detection of multiple source classes: the anisotropyenergy spectrum (i.e. the dependence of the auto-and cross-APS on the energy) is not featurelessand it is best fitted by two populations ofsources with broken-power-law energy spectra.The interpretation in terms of only one sourcepopulation (whether emitting as a power lawor broken power law) is excluded at 95% CL.This suggests that the auto- and cross-APS resultfrom a class of objects emitting mainly at lowenergies with a soft energy spectrum ∝ E − α with α ∼ α ∼ • Presence of an high-energy cut-off: our best-fit interpretation shows a cut-off at around 85GeV. This may be related to the absorption ofthe extragalactic background light (EBL), since asimilar feature is detected in the intensity energyspectrum of the IGRB in Ref. [4] . If we were The cut-off in the IGRB energy spectrum detected in Ref. [4] is atslightly higher energies, namely around 200 GeV, depending onthe model of the diffuse Galactic foreground employed. Notice,however, that the measurement in Ref. [4] is performed maskingthe sources in the 2FGL, while the value of E b =85 GeV quotedabove refers to the measurement after masking the sources in the3FGL. Thus, it is possible for the two cut-offs to be located atdifferent energies. able to confirm that the cut-off is associated withthe EBL, this would be the first time that theabsorption by the EBL is detected via anisotropies.One way to achieve such a confirmation wouldbe to detect a significant cross-correlation betweenthe same data set employed here and a catalog oftracers of the Large-Scale Structure of the Universe.The possibility of binning the catalog in redshiftwould allow us to perform a tomographic analysisand to select the emission coming from differentcomoving distances [21, 22]. Alternatively, thecut-off may be an intrinsic feature of the energyspectrum of the sources responsible for the auto-and cross-APS at high energies. • Systematic uncertaintiesin the anisotropies induced by annihilating DM: in the case of an annihilating DM candidate, anuncertainty of a factor 4 in the mass of the MWinduces a variation of a factor ∼
30 in the auto-and cross-APS associated with Galactic subhalos.For a MW mass of the order of 10 M ⊙ , Galacticsubhalos dominate the expected anisotropic signalfrom DM. If the MW is less massive, i.e., few times10 M ⊙ , the extragalactic component starts to beimportant, at least for large subhalo boost models.For DM annihilations occuring in extragalacticDM halos and subhalos, the uncertainties on thesubhalo boost factor (for a fixed M min ) inducean uncertainty of a factor ∼
20 on the expectedauto- and cross-APS from this component. Thegamma-ray emission produced by DM annihilationsin the smooth halo of the MW generates a negligibleanisotropic signal outside the adopted mask. Theoverall uncertainty on the predicted DM-inducedAPS (for a fixed M min ) is of a factor of 20, similarto the one estimated in Ref. [31] for the intensity ofall-sky gamma-ray emission. Changing the value of M min from 10 − to 1 M ⊙ approximately doublesthe systematic uncertainty. • Systematic uncertaintiesin the anisotropies induced by decaying DM:
Inthe case of decaying DM, the extragalactic signaldominates the expected auto- and cross-APS andthe prediction is independent of the value of thesubhalo boost factor. Decays in the smooth halo ofthe MW or in its subhalos are subdominant. Theoverall uncertainty (for a fixed value of M min ) isless than a factor of 2. Varying M min over the rangementioned before has a negligible effect in the caseof decaying DM. • Conservative exclusion limits on DM: requiringthat the DM-induced auto- and cross-APS doesnot exceed the measurement in any energy bin orcombination of energy bins yields an upper limit on h σ ann v i that is at least a factor of 2 less stringentthat the one obtained in Ref. [31] from the analysis8of the IGRB energy spectrum (for the REF scenarioand the b channel). In the case of annihilations into b ¯ b , the constraint on the annihilation cross sectionreaches a value as low as 10 − cm s − for a DMmass of 5 GeV and, therefore, it is approximatelytwo orders of magnitude less constraining thanthe one inferred from the observation of dwarfspheroidals galaxies. For decaying DM, the lowerlimit on τ is a factor of 5 weaker than the onefrom the IGRB intensity [96] and, at least, a factorof 5 weaker than the one from the analysis ofdwarf spheroidal galaxies [97] (for the REF caseand decays into b ¯ b ). • Exclusion limits from the 2-component fit: fittingthe data with a 2-component model that includesDM provides more constraining exclusion limits.The resulting upper limit for annihilating DM (inthe REF scenario for a M min of 10 − M ⊙ /h andannihilations into b ¯ b ) is still a factor of 10 lessconstraining than the combined analysis of dwarfspheroidals from Ref. [95]. However, below a DMmass of 30 GeV, it is slightly better than whatwas derived in Ref. [98] from the analysis of theIGRB intensity energy spectrum in terms of a 2-component model. For decaying DM, the lowerlimits on τ are, at least, a factor of 2 less stringentthan those obtained from the IGRB intensityenergy spectrum [96] or from the combined analysisof dwarf spheroidals [97].The exclusion limits on DM, although they do notexclude new regions of the DM parameter space, arecomplementary to those computed from the intensity ofthe IGRB or from the observation of dwarf spheroidalsand, therefore, they provide independent information.Also, they are expected to become more stringentas the measurement of the auto- and cross-APS willimprove during the next years. Beside making use ofthe data collected after May 2015, future analyses willrely on Pass 8 data, benefiting from the new eventclasses and data selections available (see Sec. V C).Also, future catalogs of sources, deeper than 3FGL, willexplore faint sources that are now unresolved and willimprove our modelling of those source classes. It willcertainly be interesting to complement the measurementof gamma-ray anisotropies performed here with a similarobservation at higher energies (which will be possible inthe near future with the Cherenkov Telescope Array [102,103]) or in the sub-GeV regime (with future satelliteslike ASTROGAM and ComPair [104]). Finally, in amulti-messenger perspective, the study of gamma-rayanisotropies can be interfaced with similar analyses onthe high-energy neutrinos recently discovered by IceCube[105–108]. Since the same sources that contribute to the http://astrogam.iaps.inaf.it/ IGRB (e.g. blazars, star-forming or radio galaxies) areexpected to emit also neutrinos, the auto- and cross-APSmeasured in this work represents a useful indication forthe minimal level of anisotropies that can be found inthe distribution of neutrinos. A quantitative estimate ofIceCube prospects to detect anisotropies can be found inRef. [109].
ACKNOWLEDGMENTS
The
Fermi
LAT Collaboration acknowledges generousongoing support from a number of agencies and institutesthat have supported both the development and theoperation of the LAT as well as scientific data analysis.These include the National Aeronautics and SpaceAdministration and the Department of Energy in theUnited States, the Commissariat `a l’Energie Atomiqueand the Centre National de la Recherche Scientifique /Institut National de Physique Nucl´eaire et de Physiquedes Particules in France, the Agenzia Spaziale Italianaand the Istituto Nazionale di Fisica Nucleare in Italy,the Ministry of Education, Culture, Sports, Science andTechnology (MEXT), High Energy Accelerator ResearchOrganization (KEK) and Japan Aerospace ExplorationAgency (JAXA) in Japan, and the K. A. WallenbergFoundation, the Swedish Research Council and theSwedish National Space Board in Sweden.Additional support for science analysis during theoperations phase is gratefully acknowledged from theIstituto Nazionale di Astrofisica in Italy and the CentreNational d’´Etudes Spatiales in France.We thank Dr. Shin’ichiro Ando for useful discussionand for providing us with the lower limits on the decaylifetime from Ref. [96]. We also thanks Prof. John F.Beacom for the valuable discussion during the last stagesof the work.MF gratefully acknowledges support from theNetherlands Organization for Scientific Research (NWO)through a Vidi grant (P.I.: Dr. Shin’ichiro Ando), fromthe the Leverhulme Trust and the project MultiDarkCSD2009-00064. The Dark Cosmology Centre is fundedby the DNRF. JZ is supported by the EU under a MarieCurie International Incoming Fellowship, contract PIIF-GA-2013-62772. JG acknowledges support from NASAthrough Einstein Postdoctoral Fellowship grant PF1-120089 awarded by the Chandra X-ray Center, which isoperated by the Smithsonian Astrophysical Observatoryfor NASA under contract NAS8-03060, and from a MarieCurie International Incoming Fellowship, contract PIIF-GA-2013-628997. GAGV is supported by CONICYTFONDECYT/POSTODOCTORADO/3160153, theSpanish MINECO’s Consolider Ingenio 2010 Programmeunder grant MultiDark CSD2009-00064 and also partlyby MINECO under grant FPA2012-34694. MASC is aWenner-Gren Fellow and acknowldeges the support of theWenner-Gren Foundantions to develop his research. EKthanks J.U. Lange for useful discussion on the bias in C P Appendix A: Derivation of the photon noise C N Let n i be the number of photons in a pixel i , ¯ n i beits expectation value, and δn i ≡ n i − ¯ n i be a fluctuationaround the mean. The photon flux is given by n i /A i ,where A i is the exposure in pixel i . Then, the photonflux per unit solid angle is given by n i / ( A i Ω pix ) whereΩ pix is the solid angle of pixel i .The spherical harmonics coefficients a ℓ,m are a ℓ,m = Z d Ω i δn i A i Ω pix Y ∗ ℓ,m (Ω i ) , (A1)where Ω i denotes the direction of pixel i . The expectationvalue of the product between two coefficients is h a ℓ,m a ∗ ℓ ′ ,m ′ i = Z d Ω i Z d Ω j h δn i δn j i A i A j Ω Y ∗ ℓm (Ω i ) Y ℓ ′ m ′ (Ω j ) . (A2)If δn i is purely Poisson noise, h δn i δn j i / Ω =(¯ n i / Ω pix ) δ (Ω i − Ω j ) where δ is the Dirac delta function.Thus: h a ℓ,m a ∗ ℓ ′ ,m ′ i = Z d Ω i ¯ n i A i Ω pix Y ℓm (Ω i ) Y ∗ ℓ ′ m ′ (Ω i ) . (A3)Now, we calculate the diagonal element, i.e., C ℓ ≡ P m h a ℓ,m a ∗ ℓ,m i / (2 ℓ + 1), obtaining: C ℓ = Z d Ω i π ¯ n i A i Ω pix . (A4)The latter is independent on multipole ℓ and equivalentto the definition of C k N in Eq. 5. Appendix B: Auto-correlation angular spectra for allthe energy bins
Figs. 29 and 30 show the binned auto-APS C ℓ obtainedas described in Sec. III, for all 13 energy bins considered.The auto-APS is shown only within the signal region,i.e. between a multipole of 49 and 706. Red circles referto the data set obtained with our default mask coveringthe sources in 3FGL and the solid red line marks thecorresponding best-fit C P . The pink band denotes the68% CL error on C P Appendix C: Anisotropy energy spectrum for thecross-correlation angular power spectrum
Figs. 31 and 32 show the best-fit C P Appendix D: The dependence on energy of thecross-correlation coefficients
Figs. 33 and 34 show the cross-correlation coefficients r i,j defined in Sec. V B in terms of the best-fit auto-and cross-APS C P . Each panel shows r i,j at a specificenergy E i , as a function of E j . Red circles refer tothe mask covering 3FGL sources and blue triangles tothe mask around 2FGL sources. The solid black lineshows the cross-correlation coefficents corresponding tothe best-fit solution discussed in Sec. VI in the case of twopopulations of unresolved sources with broken-power-lawenergy spectra and masking 3FGL sources. The fact thatthe blue triangles decrease with energy in the first panels,while they increase towards 1 in the last panels indicatesthe lack of correlation between low and high energies.The same trend is noted for the red circles, but with alower significance. Appendix E: Exclusion limit on Dark Matter for the τ and µ channel Sec. VIII shows exclusion limits on the DM h σ ann v i and τ in the case of annihilations/decay into b ¯ b . Herewe calculate the same exclusion limits for two additionalchannels. Fig. 35 shows the upper limits on h σ ann v i (left panel) and on τ (right panel) as a function of theDM mass m χ , in the case of annihilations/decays into τ + τ − . The exclusion limits are obtained following theconservative approach described in Sec. VIII A. The solidblack line refers to the REF scenario, while the solidred and solid blue ones stand for the MAX and MINbenchmark. The solid red and solid black lines almostexactly overlap in the right panel. The dashed black, blueand red lines are obtained considering only the auto-APSmeasurement. The red and blue shaded band indicatethe variability of the exclusion limits in the MAX andMIN scenario when M min is left free to vary between 1 M ⊙ and 10 − M ⊙ . The long-dashed grey line in the leftpanel shows the thermal annihilation cross section, ascomputed in Ref. [94], while the dot-dashed grey line isthe upper limit obtained in Ref. [95] from the analysis of0 Multipole s r ] - s r - s - [ c m l C -80-60-40-20020406080 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -80-60-40-20020406080 -18 × Energy bin [0.50-0.72] GeV Multipole s r ] - s r - s - [ c m l C -30-20-100102030 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -30-20-100102030 -18 × Energy bin [0.72-1.04] GeVMultipole s r ] - s r - s - [ c m l C -5-4-3-2-1012345 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -5-4-3-2-1012345 -18 × Energy bin [1.04-1.38] GeV Multipole s r ] - s r - s - [ c m l C -3-2-10123 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -3-2-10123 -18 × Energy bin [1.38-1.99] GeVMultipole s r ] - s r - s - [ c m l C -1-0.8-0.6-0.4-0.200.20.40.60.81 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -1-0.8-0.6-0.4-0.200.20.40.60.81 -18 × Energy bin [1.99-3.15] GeV Multipole s r ] - s r - s - [ c m l C -0.5-0.4-0.3-0.2-0.100.10.20.30.40.5 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -0.5-0.4-0.3-0.2-0.100.10.20.30.40.5 -18 × Energy bin [3.15-5.00] GeV
FIG. 29. Auto-APS of the IGRB for the first 6 energy bins and for the reference data set (P7REP ULTRACLEAN V15 frontevents) using the reference mask which excludes | b | < ◦ and 3FGL sources (red circles). The blue triangles show the samebut masking the sources in 2FGL, instead. Data have been binned as described in Sec. IV A. The solid red line shows thebest-fit C P for the red data points, with the pink band indicating its 68% CL error. The dashed blue line corresponds to thebest-fit C P for the blue data points. The energy range is indicated on the top of each panel. Note that only the results in oursignal region (i.e. between ℓ = 49 and 706) are plotted and that the scale of the y -axis can vary from panel to panel. Multipole s r ] - s r - s - [ c m l C -0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 -18 × Energy bin [5.00-7.23] GeV Multipole s r ] - s r - s - [ c m l C -60-40-200204060 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -60-40-200204060 -21 × Energy bin [7.23-10.45] GeVMultipole s r ] - s r - s - [ c m l C -0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 -18 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -0.1-0.08-0.06-0.04-0.0200.020.040.060.080.1 -18 × Energy bin [10.45-21.83] GeV Multipole s r ] - s r - s - [ c m l C -15-10-5051015 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -15-10-5051015 -21 × Energy bin [21.83-50.00] GeVMultipole s r ] - s r - s - [ c m l C -10-8-6-4-20246810 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -10-8-6-4-20246810 -21 × Energy bin [50.00-95.27] GeV Multipole s r ] - s r - s - [ c m l C -2-1.5-1-0.500.511.52 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -2-1.5-1-0.500.511.52 -21 × Energy bin [95.27-199.05] GeV s r ] - s r - s - [ c m l C -0.5-0.4-0.3-0.2-0.100.10.20.30.40.5 -21 × Masking sources in 3FGLMasking sources in 2FGLPoissonian fit (masking sources in 3FGL)Poissonian fit (masking sources in 2FGL) -0.5-0.4-0.3-0.2-0.100.10.20.30.40.5 -21 × Energy bin [199.05-500.00] GeV [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -16 -19 -18 -17 -16 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[0.50-0.72] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -16 -19 -18 -17 -16 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[0.72-1.04] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -16 -19 -18 -17 -16 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[1.04-1.38] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[1.38-1.99] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[1.99-3.15] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[3.15-5.00] GeV i E FIG. 31. Depedence of the cross-APS on the energy. Each panel shows the best-fit Poissonian C P for the cross-APS betweenthe i -th and the j -th energy bins, as a function of E j . Red circles are for the reference data set (P7REP ULTRACLEAN V15front events) using the default mask masking 3FGL sources, while the blue triangles show the result for the same data set andfor the default mask excluding 2FGL sources. The first 6 energy bins are shown in this figure and E i is indicated in the top ofeach panel. The solid black line is the best-fit solution when data are fitted assuming two independent populations of sourceswith broken-power-law energy spectra. The short-dashed and long-dashed black lines show the two populations independently. [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[5.00-7.23] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[7.23-10.45] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[10.45-21.83] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[21.83-50.00] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[50.00-95.27] GeV i E [GeV] j Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -19 -18 -17 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[95.27-199.05] GeV i E [GeV]Energy E1 10 s r ] - s r - s - c m [ G e V - ) j E ∆ ( - ) i E ∆ ( P C j E i E -20 -19 -18 -17 -20 -19 -18 -17 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs=[199.05-500.00] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[0.50-0.72] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[0.72-1.04] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[1.04-1.38] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[1.38-1.99] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[1.99-3.15] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[3.15-5.00] GeV i E FIG. 33. Dependence of the cross-correlation coefficients on the energy. Each panel shows the cross-correlation coefficients r i,j defined in Sec. V B between the i -th and the j -th energy bins, as a function of E j . Red circles are for the reference data set(P7REP ULTRACLEAN V15 front events) using the default mask masking 3FGL sources, while the blue triangles show theresult for the same data set and for the default mask excluding 2FGL sources. The first 6 energy bins are shown in this figureand E i is indicated in the top of each panel. The solid black line shows the r i,j corresponding to the best-fit solution when dataare fitted masking 3FGL sources and assuming two independent populations of sources with broken-power-law energy spectra(see Sec. VI). [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[5.00-7.23] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[7.23-10.45] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[10.45-21.83] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[21.83-50.0] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[50.0-95.27] GeV i E [GeV] j Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[95.27-199.05] GeV i E [GeV]Energy E1 10 ) j ( E i , j r -0.500.511.52 -0.500.511.52 Masking sources in 2FGLMasking sources in 3FGLBest-fit with 2 broken PLs =[199.05-500.00] GeV i E
615 dwarf spheroidal galaxies. Finally, the short-dashedgrey line derives from the analysis of the IGRB intensityperformed in Ref. [31]. On the other hand, the short-dashed grey line in the right panel of Fig. 35 is obtainedfrom the study of the IGRB intensity in Ref. [96] andthe dot-dashed grey line comes from the observation of15 dwarf spheroidal galaxies performed in Ref. [97].Fig. 36 shows the same exclusion limits as in Fig 35but for annihilations/decays into the µ + µ − . Betweenapproximately 20 and 200 GeV, the DM-induced signal isdominated by the IC emission associated with the smoothhalo of the MW. That is the reason why the solid blackand red lines overlap, since the REF and MAX scenariosonly differ in the computation of the boost factor for theextragalactic component. For the same reason the blueand red shaded bands are reduced in width. Above 200GeV, the IC emission for the extragalactic componentstarts to contribute more and the solid black and redlines deviate one from the other again.The 2-component model developed in Sec.VIII B isused to fit the measured auto-APS and cross-APS, fordifferent values of DM mass, annihilation cross sectionor decay lifetime. Fig. 37 shows the TS defined as thedifference between the χ of the best fit for the nullhypothesis (i.e. with no DM) and the χ of the best fitin the case with the DM component. The top panelsare for annihilation/decay into τ + τ − and the bottomones for the µ -channel. The ones on the right are for an annihilating DM candidate and the ones on the leftfor decaying DM. They all refer to the REF scenario with M min = 10 − M ⊙ /h . As indicated in the labels, the whitelines determine the 68%, 90% and 95% CL regions.Assuming that the measured auto- and cross-APS arewell described simply by a Poissonian component, the2-component model is used to derive exclusion limitson DM as done in Sec. VIII B but in the case ofannihilations/decays into τ + τ − (Fig. 38) and into µ + µ − (Fig. 39). In both figures the left panel is for annihilatingDM and the right one for decaying DM. The solidblack, red and blue line show the REF, MAX and MINscenario for M min = 10 − M ⊙ /h , respectively and theblue and shaded areas around the corresponding solidlines indicate how the limits change when M min is leftfree to vary. The dashed black is the exclusion limitin the conservative scenario, from Fig. 35 and 36. Inthe left panels, the long-dashed grey line is the thermalannihilation cross section from Ref. [94] and the dot-dashed line is the upper limit derived in Ref. [95] fromthe combined analysis of 15 dwarf spheroidals. Also, theshort-dashed grey line comes from the analysis of theIGRB intensity performed in Ref. [98]. On the otherhand, in the right panels, the short-dashed grey linerepresents the lower limit from Fig. 5 of Ref. [96] fromthe IGRB intensity. The dot-dashed grey line is the lowerlimit from the analysis of 15 dwarf spheroidal galaxiesperformed in Ref. [97]. [1] M. Ackermann et al. (Fermi-LAT), Phys. Rev. D85 , 083007 (2012),arXiv:1202.2856 [astro-ph.HE].[2] A. A. Abdo et al. (Fermi-LAT),Astrophys. J. Suppl. , 405 (2010),arXiv:1002.2280 [astro-ph.HE].[3] A. A. Abdo et al. (Fermi-LAT),Phys. Rev. Lett. , 101101 (2010),arXiv:1002.3603 [astro-ph.HE].[4] M. Ackermann et al. (Fermi-LAT),Astrophys. J. , 86 (2015),arXiv:1410.3696 [astro-ph.HE].[5] M. Fornasaand M. A. Snchez-Conde, Phys. Rept. , 1 (2015),arXiv:1502.02866 [astro-ph.CO].[6] P.-J.Zhang and J. F. Beacom, Astrophys. J. , 37 (2004),arXiv:astro-ph/0401351 [astro-ph].[7] S. Ando, E. Komatsu, T. Narumoto, and T. Totani,Mon.Not.Roy.Astron.Soc. , 1635 (2007),arXiv:astro-ph/0610155 [astro-ph].[8] S. Ando and V. Pavlidou,Mon. Not. Roy. Astron. Soc. , 2122 (2009),arXiv:0908.3890 [astro-ph.HE].[9] J. M.Siegal-Gaskins, R. Reesman, V. Pavlidou, S. Profumo,and T. P. Walker, Mon.Not.Roy.Astron.Soc. , 1074S(2011), arXiv:1011.5501 [astro-ph.HE].[10] A. Cuoco, E. Komatsu, and J. M. Siegal-Gaskins, Phys. Rev.
D86 , 063004 (2012),arXiv:1202.5309 [astro-ph.CO].[11] J. P.Harding and K. N. Abazajian, JCAP , 026 (2012),arXiv:1206.4734 [astro-ph.HE].[12] M. Di Mauro, F. Calore, F. Donato, M. Ajello, andL. Latronico, Astrophys. J. , 161 (2014),arXiv:1304.0908 [astro-ph.HE].[13] F. Zandanel, I. Tamborra, S. Gabici, and S. Ando,Astron. Astrophys. , A32 (2015),arXiv:1410.8697 [astro-ph.HE].[14] S. Ando and E. Komatsu,Phys. Rev.
D87 , 123539 (2013),arXiv:1301.5901 [astro-ph.CO].[15] D. Hooper,T. Linden, and A. Lopez, JCAP , 019 (2016),arXiv:1604.08505 [astro-ph.HE].[16] M. Di Mauro, A. Cuoco, F. Donato, and J. M. Siegal-Gaskins, JCAP , 021 (2014),arXiv:1407.3275 [astro-ph.HE].[17] J.-Q. Xia, A. Cuoco, E. Branchini, and M. Viel,Astrophys. J. Suppl. , 15 (2015),arXiv:1503.05918 [astro-ph.CO].[18] M. Shirasaki, S. Horiuchi, and N. Yoshida,Phys. Rev.
D92 , 123540 (2015),arXiv:1511.07092 [astro-ph.CO].[19] A. Cuoco,J.-Q. Xia, M. Regis, E. Branchini, N. Fornengo, andM. Viel, Astrophys. J. Suppl. , 29 (2015), [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)Fermi LAT IGRB intensity (Ackermann et al., conservative, 2015)REF scenarioMAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 - τ + τ Annihilation, [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, conservative, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioMAX scenarioMIN scenario - τ + τ Decay,
FIG. 35. Conservative exclusion limits on annihilating and decaying DM from the new APS measurement, for the τ channel. Left:
The solid lines show the upper limits on h σ ann v i derived from the auto- and cross-APS measured in Sec. III, as a functionof m χ , for M min = 10 − M ⊙ and annihilations into τ + τ − . The limits follow the conservative approach described in Sec. VIII A.The black line is for the REF scenario, while the red and blue ones are for MAX and MIN. The grey band between the MINand MAX scenario represents our estimated total astrophysical uncertainty for M min = 10 − M ⊙ , accounting for all the sourcesof uncertainty mentioned in Sec. VII. The red and blue shaded bands describe the effect of changing M min between 10 − M min and 1 M min , for the MAX and MIN scenario, respectively. In the case of the black, red and blue dashed lines, the upper limitsare derived only considering the measured auto-APS and neglecting the cross-APS. For comparison, the long-dashed grey linemarks the annihilation cross section for thermal relics from Ref. [94] and the dash-dotted grey line the upper limit obtained inRef. [95] from the combined analysis of 15 dwarf spheroidal galaxies. Finally, the short-dashed grey line shows the conservativeupper limit derived in Ref. [31] from the intensity of the IGRB. Right : The same as in the left panel but for the lower limitson τ for decaying DM. The short-dashed grey line represents the lower limit obtained in Fig. 6 of Ref. [96] from the IGRBintensity, while the dash-dotted grey one is obtained from the combined analysis of 15 dwarf spheroidal galaxies in Ref. [97]. [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 -20 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)Fermi LAT IGRB intensity (Ackermann et al., conservative, 2015)REF scenarioMAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 -20 - µ + µ Annihilation, [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, conservative, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioMAX scenarioMIN scenario - µ + µ Decay,
FIG. 36. Same as Fig. 35 but for annihilations/decays into µ + µ − .arXiv:1506.01030 [astro-ph.HE].[20] M. Regis,J.-Q. Xia, A. Cuoco, E. Branchini, N. Fornengo,and M. Viel, Phys. Rev. Lett. , 241301 (2015),arXiv:1503.05922 [astro-ph.CO].[21] S. Ando, A. Benoit- Lvy, and E. Komatsu, Phys. Rev. D90 , 023514 (2014),arXiv:1312.4403 [astro-ph.CO].[22] S. Ando, JCAP , 061 (2014),arXiv:1407.8502 [astro-ph.CO].[23] M. Shirasaki, S. Horiuchi, and N. Yoshida,Phys. Rev.
D90 , 063502 (2014), mass [GeV]10 ] - s v > [ c m ann σ < -27 -26 -25 -24 -23 -22
10 -4-20246
68% CL
90% CL95% CL - τ + τ Annihilation, mass [GeV]10 [ s ] τ
10 -3-2-10123456
68% CL90% CL95% CL - τ + τ Decay, mass [GeV]10 ] - s v > [ c m ann σ < -25 -24 -23 -22 -21
10 -3-2-10123456
68% CL 90% CL
95% CL - µ + µ Annihilation, mass [GeV]10 [ s ] τ
10 -3-2-10123456
68% CL90% CL95% CL - µ + µ Decay,
FIG. 37. ∆ χ between the best-fit solution for the 2-component scenario and the best fit of the null hypothesis. Resultspresented here refer to the REF scenario with M min = 10 − M ⊙ /h and annihilation/decay into τ + τ − (top panels) or µ + µ − (bottom panels). The panels on the left are for annihilating DM and the ones on the right for decaying DM. Each point in thebi-dimensional parameter space is colored according to its ∆ χ , i.e. the difference between the χ of the best fit to the auto-and cross-APS in terms of the 2-component model and the χ of the best fit of the null hypothesis (i.e. no DM). The closedwhite contour marks the 68% CL region. The 90% and 95% CL ones in the left (right) panels contain all the region below(above) the white open curves labelled “90% CL” and “95% CL”.arXiv:1404.5503 [astro-ph.CO].[24] N. Fornengo, L. Perotto, M. Regis, and S. Camera,Astrophys. J. , L1 (2015),arXiv:1410.4997 [astro-ph.CO].[25] D. Malyshevand D. W. Hogg, Astrophys. J. , 181 (2011), arXiv:1104.0010 [astro-ph.CO].[26] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo,and A. Vittino, Astrophys. J. Suppl. , 18 (2016),arXiv:1512.07190 [astro-ph.HE].[27] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo, andM. Regis, Astrophys. J. , L31 (2016), [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 -20 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)Fermi LAT IGRB intensity (Ajello et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 -20 - τ + τ Annihilation, [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, power-law, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario - τ + τ Decay,
FIG. 38. Exclusion limits on annihilating and decaying DM (for the τ channel) from the fit to the binned C ℓ in terms ofthe 2-component model. Left:
The solid lines show the upper limits that can be derived on h σ ann v i as a function of m χ (forannihilation into τ + τ − quarks and M min = 10 − M ⊙ ) by fitting the Fermi
LAT data with a 2-component model that includesastrophysical sources and DM (see text for details). The black, blue and red lines correspond to the REF, MIN and MAXscenario. The blue and red shaded areas indicate how the MIN and MAX upper limits change when leaving M min free to varybetween 10 − M ⊙ and 1 M ⊙ . The black dashed line is the REF upper limit in the conservative case, from Fig. 35, while thelong-dashed grey line is the thermal annihilation cross section from Ref. [94]. The dot-dashed line is the upper limits derivedin Ref. [95] from the combined analysis of 15 dwarf spheroidals, while the short-dashed grey line comes from the analysis ofthe IGRB intensity performed in Ref. [98]. Right : The same as in the left panel but for the lower limits on τ , in the case ofdecaying DM. The short-dashed grey line represents the lower limit obtained in Ref. [96] from the IGRB intensity. The lineis taken from Fig. 5 of Ref. [96], where the IGRB is interpreted in terms of a component with a power-law emission spectrumand a DM contribution. Finally, the dot-dashed grey line is the upper limit from the analysis of 15 dwarf spheroidal galaxiesperformed in Ref. [97]. [GeV] χ m10 ] - s v > [ c m ann σ < -26 -25 -24 -23 -22 -21 -20 v) (Steigman et al., 2012) ann σ Thermal (Fermi LAT dwarf Spheroidals (Ackermann et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario (MAX) min
Uncertainty on M (MIN) min
Uncertainty on M -26 -25 -24 -23 -22 -21 -20 - µ + µ Annihilation, [GeV] χ m10 [ s ] τ Fermi LAT IGRB intensity (Ando & Ishiwata, power-law, 2015)Fermi LAT dwarf Spheroidals (Baring et al., 2015)REF scenarioREF scenario (conservative)MAX scenarioMIN scenario - µ + µ Decay,
FIG. 39. Same as in Fig. 38 but for annihilations/decays into µ + µ − .arXiv:1605.04256 [astro-ph.HE].[28] J. M. Siegal-Gaskinsand V. Pavlidou, Phys. Rev. Lett. , 241301 (2009),arXiv:0901.3776 [astro-ph.HE].[29] J. Zavala, V. Springel, and M. Boylan-Kolchin,Mon. Not. Roy. Astron. Soc. , 593 (2010), arXiv:0908.2428 [astro-ph.CO].[30] A. Cuoco, A. Sellerholm, J. Conrad, and S. Hannestad,Mon. Not. Roy. Astron. Soc. , 2040 (2011),arXiv:1005.0843 [astro-ph.HE].[31] M. Ackermann et al. (Fermi-LAT),JCAP , 008 (2015), arXiv:1501.05464 [astro-ph.CO].[32] G. A. Gomez-Vargas, A. Cuoco, T. Linden, M. A.Sanchez-Conde, J. M. Siegal-Gaskins, T. Delahaye,M. Fornasa, E. Komatsu, F. Prada, and J. Zavala(Fermi-LAT), Proceedings, 4th Roma InternationalConference on Astro-Particle Physics (RICAP 13) ,Nucl. Instrum. Meth.
A742 , 149 (2014).[33] J. U. Lange and M. C. Chu,Mon. Not. Roy. Astron. Soc. , 939 (2015),arXiv:1412.5749 [astro-ph.CO].[34] M. Fornasa, J. Zavala, M. A. Sanchez-Conde, J. M.Siegal-Gaskins, T. Delahaye,F. Prada, M. Vogelsberger, F. Zandanel, and C. S.Frenk, Mon. Not. Roy. Astron. Soc. , 1529 (2013),arXiv:1207.0502 [astro-ph.HE].[35] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt,F. K. Hansen, M. Reinecke, and M. Bartelman,Astrophys. J. , 759 (2005),arXiv:astro-ph/0409513 [astro-ph].[36] M. Ackermann,G. Johannesson, S. Digel, I. V. Moskalenko, T. Porter,O. Reimer, and A. Strong (Fermi-LAT), ,AIP Conf. Proc. , 763 (2009).[37] A. Challinor, G. Chon, S. Colombi, H. Eric, S. Prunet,and I. Szaput,
PolSpice: Spatially InhomogeneousCorrelation Estimator for Temperature and Polarisation (Astrophysics Source Code Library, 2011).[38] F. Acero, Astrophys. J. Suppl. , 23 (2015).[39] P. L. Nolan et al. (Fermi-LAT),Astrophys. J. Suppl. , 31 (2012),arXiv:1108.1435 [astro-ph.HE].[40] F. Acero et al. (Fermi-LAT),Astrophys. J. Suppl. , 26 (2016),arXiv:1602.07246 [astro-ph.HE].[41] L. Knox, Phys. Rev.
D52 , 4307 (1995),arXiv:astro-ph/9504054 [astro-ph].[42] S. S. Campbell,Mon. Not. Roy. Astron. Soc. , 2854 (2015),arXiv:1411.4031 [astro-ph.CO].[43] B. S. Hensley, J. M. Siegal-Gaskins, and V. Pavlidou,Astrophys. J. , 277 (2010),arXiv:0912.1854 [astro-ph.CO].[44] B. S. Hensley, V. Pavlidou, and J. M. Siegal-Gaskins, Mon. Not. Roy. Astron. Soc. , 591 (2013),arXiv:1210.7239 [astro-ph.CO].[45] S. Ando and E. Komatsu,Phys. Rev.
D73 , 023521 (2006),arXiv:astro-ph/0512217 [astro-ph].[46] N. Fornengo and M. Regis, Front. Physics , 6 (2014),arXiv:1312.4835 [astro-ph.CO].[47] W. Atwood et al. (Fermi-LAT) (2013)arXiv:1303.3514 [astro-ph.IM].[48] A. A. Abdo et al. (Fermi-LAT), Astrophys. J. , 116 (2011),arXiv:1104.2093 [astro-ph.HE].[49] K. C. Y. Ng, J. F. Beacom, A. H. G.Peter, and C. Rott, Phys. Rev. D94 , 023004 (2016),arXiv:1508.06276 [astro-ph.HE].[50] D. Seckel, T. Stanev, and T. K. Gaisser,Astrophys. J. , 652 (1991).[51] E. Orlando and A. Strong,
TheMulti-Messenger Approach to High-Energy Gamma- Ray Sources , Astrophys. Space Sci. , 359 (2007),arXiv:astro-ph/0607563 [astro-ph].[52] I. V. Moskalenko, T. A. Porter, and S. W. Digel,Astrophys. J. , L65 (2006), [Erratum: Astrophys.J.664,L143(2007)], arXiv:astro-ph/0607521 [astro-ph].[53] E. Orlando and A. Strong, (2013),arXiv:1307.6798 [astro-ph.HE].[54] M. Ackermann et al. (Fermi-LAT), Phys. Rev.
D93 , 082001 (2016),arXiv:1604.03349 [astro-ph.HE].[55] A. E. Broderick, C. Pfrommer,E. Puchwein, K. M. Smith, and P. Chang (HeidelbergInstitute for Theoretical Studies, Perimeter Institutefor Theoretical Physics), Astrophys. J. , 12 (2014),arXiv:1308.0015 [astro-ph.CO].[56] F. Feroz and M. P.Hobson, Mon. Not. Roy. Astron. Soc. , 449 (2008),arXiv:0704.3704 [astro-ph].[57] F. Feroz, M. P. Hobson, and M. Bridges,Mon. Not. Roy. Astron. Soc. , 1601 (2009),arXiv:0809.3437 [astro-ph].[58] F. Feroz, M. P. Hobson, E. Cameron, and A. N. Pettitt,(2013), arXiv:1306.2144 [astro-ph.IM].[59] R. Trotta, Contemp. Phys. , 71 (2008),arXiv:0803.4089 [astro-ph].[60] M. Boylan-Kolchin, V. Springel, S. D. M. White,A. Jenkins, and G. Lemson,Mon. Not. Roy. Astron. Soc. , 1150 (2009),arXiv:0903.3041 [astro-ph.CO].[61] J. F. Navarro, A. Ludlow, V. Springel, J. Wang,M. Vogelsberger, S. D. M. White, A. Jenkins, C. S.Frenk, andA. Helmi, Mon. Not. Roy. Astron. Soc. , 21 (2010),arXiv:0810.1522 [astro-ph].[62] V. Springel, J. Wang, M. Vogelsberger, A. Ludlow,A. Jenkins,A. Helmi, J. F. Navarro, C. S. Frenk, and S. D. M.White, Mon. Not. Roy. Astron. Soc. , 1685 (2008),arXiv:0809.0898 [astro-ph].[63] S. Profumo and T. E. Jeltema, JCAP , 020 (2009),arXiv:0906.0001 [astro-ph.CO].[64] A. Dominguez et al. , Mon. Not. Roy. Astron. Soc. , 2556 (2011),arXiv:1007.1459 [astro-ph.CO].[65] J. F. Navarro, C. S. Frenk, and S. D. M. White,Astrophys. J. , 493 (1997),arXiv:astro-ph/9611107 [astro-ph].[66] B. Moore, T. R. Quinn, F. Governato, J. Stadel, andG. Lake,Mon. Not. Roy. Astron. Soc. , 1147 (1999),arXiv:astro-ph/9903164 [astro-ph].[67] A. Burkert, IAU Symposium 171: New Light onGalaxy Evolution Heidelberg, Germany, June 26-30, 1995 , IAU Symp. , 175 (1996), [Astrophys.J.447,L25(1995)], arXiv:astro-ph/9504041 [astro-ph].[68] P. A. R. Ade et al. (Planck), (2015),arXiv:1502.01589 [astro-ph.CO].[69] Q. Guo, S. White, R. E. Angulo, B. Henriques,G. Lemson, M. Boylan-Kolchin, P. Thomas, andC. Short,Mon. Not. Roy. Astron. Soc. , 1351 (2013),arXiv:1206.0052 [astro-ph.CO].[70] A. D. Ludlow, J. F. Navarro, R. E. Angulo, M. Boylan-Kolchin, V. Springel, C. Frenk, and S. D. M. White, Mon. Not. Roy. Astron. Soc. , 378 (2014),arXiv:1312.0945 [astro-ph.CO].[71] M. Cirelli, G. Corcella, A. Hektor, G. Hutsi,M. Kadastik, P. Panci, M. Raidal, F. Sala,and A. Strumia, JCAP , 051 (2011), [Erratum:JCAP1210,E01(2012)], arXiv:1012.4515 [hep-ph].[72] S. Profumo, K. Sigurdson, andM. Kamionkowski, Phys. Rev. Lett. , 031301 (2006),arXiv:astro-ph/0603373 [astro-ph].[73] T. Bringmann, New J. Phys. , 105027 (2009),arXiv:0903.0189 [astro-ph.CO].[74] M. A. Snchez-Conde and F. Prada,Mon. Not. Roy. Astron. Soc. , 2271 (2014),arXiv:1312.1729 [astro-ph.CO].[75] D. Anderhaldenand J. Diemand, JCAP , 009 (2013), [Erratum:JCAP1308,E02(2013)], arXiv:1302.0003 [astro-ph.CO].[76] T. Ishiyama, Astrophys. J. , 27 (2014),arXiv:1404.1650 [astro-ph.CO].[77] V. Springel, S. D. M. White, C. S. Frenk, J. F. Navarro,A. Jenkins, M. Vogelsberger, J. Wang, A. Ludlow, andA. Helmi, Nature , 73 (2008).[78] L. Gao,C. S. Frenk, A. Jenkins, V. Springel, and S. D. M.White, Mon. Not. Roy. Astron. Soc. , 1721 (2012),arXiv:1107.1916 [astro-ph.CO].[79] M. Kamionkowski, S. M. Koushiappas, and M. Kuhlen,Phys. Rev. D81 , 043532 (2010),arXiv:1001.3144 [astro-ph.GA].[80] M. A. Sanchez-Conde, M. Cannoni, F. Zandanel,M. E. Gomez, and F. Prada, JCAP , 011 (2011),arXiv:1104.3530 [astro-ph.HE].[81] R. Bartels and S. Ando,Phys. Rev.
D92 , 123508 (2015),arXiv:1507.08656 [astro-ph.CO].[82] . Molin, M. A. Snchez-Conde, S. Palomares-Ruiz, andF. Prada, (2016), arXiv:1603.04057 [astro-ph.CO].[83] J. Zavala and N. Afshordi,Mon. Not. Roy. Astron. Soc. , 1317 (2014),arXiv:1308.1098 [astro-ph.CO].[84] J. Zavala and N. Afshordi,Mon. Not. Roy. Astron. Soc. , 1329 (2014),arXiv:1311.3296 [astro-ph.CO].[85] J. Zavala and N. Afshordi,Mon. Not. Roy. Astron. Soc. , 986 (2016),arXiv:1508.02713 [astro-ph.CO].[86] J. Einasto, Trudy Inst. Astroz. Alma-Ata , 87 (1965).[87] W. Wang, J. Han, A. P. Cooper, S. Cole, C. Frenk,and B. Lowing,Mon. Not. Roy. Astron. Soc. , 377 (2015),arXiv:1502.03477 [astro-ph.GA].[88] G. Bertone, M. Cirelli, A. Strumia, and M. Taoso,JCAP , 009 (2009), arXiv:0811.3744 [astro-ph].[89] S. Ando, Phys. Rev. D80 , 023520 (2009),arXiv:0903.4685 [astro-ph.CO]. [90] M. Boylan-Kolchin, V. Springel, S. D. M. White, andA. Jenkins, Mon. Not. Roy. Astron. Soc. , 896(2010), arXiv:0911.4484 [astro-ph.CO].[91] A. Klypin, S. Trujillo-Gomez, and J. Primack,Astrophys. J. , 102 (2011),arXiv:1002.3660 [astro-ph.CO].[92] M. Ackermann et al. (Fermi-LAT),Astrophys. J. , 14 (2015),arXiv:1501.06054 [astro-ph.HE].[93] M. Fornasa, L. Pieri, G. Bertone, and E. Branchini,Phys. Rev.
D80 , 023518 (2009),arXiv:0901.2921 [astro-ph.CO].[94] G. Steigman, B. Dasgupta,and J. F. Beacom, Phys. Rev.
D86 , 023506 (2012),arXiv:1204.3622 [hep-ph].[95] M. Ackermann et al. (Fermi-LAT),Phys. Rev. Lett. , 231301 (2015),arXiv:1503.02641 [astro-ph.HE].[96] S. Ando and K. Ishiwata, JCAP , 024 (2015),arXiv:1502.02007 [astro-ph.CO].[97] M. G. Baring, T. Ghosh, F. S. Queiroz, and K. Sinha,Phys. Rev.
D93 , 103009 (2016),arXiv:1510.00389 [hep-ph].[98] M. Ajello et al. , Astrophys. J. , L27 (2015),arXiv:1501.05301 [astro-ph.HE].[99] M. Ajello et al. , Astrophys. J. , 108 (2012),arXiv:1110.3787 [astro-ph.CO].[100] I. Tamborra,S. Ando, and K. Murase, JCAP , 043 (2014),arXiv:1404.1189 [astro-ph.HE].[101] A. M. Brown,C. Boehm, J. Graham, T. Lacroix, P. M. Chadwick,and J. Silk, (2016), arXiv:1603.05469 [astro-ph.HE].[102] J. Ripken, A. Cuoco, H.-S. Zechlin, J. Conrad, andD. Horns, JCAP , 049 (2014),arXiv:1211.6922 [astro-ph.HE].[103] B. S. Acharya et al. , Astropart. Phys. , 3 (2013).[104] A. A. Moiseev et al. , (2015),arXiv:1508.07349 [astro-ph.IM].[105] M. G. Aartsen et al. (IceCube), Phys. Rev. Lett. , 021103 (2013),arXiv:1304.5356 [astro-ph.HE].[106] M. G.Aartsen et al. (IceCube), Science , 1242856 (2013),arXiv:1311.5238 [astro-ph.HE].[107] M. G. Aartsen et al. (IceCube), Phys. Rev. Lett. , 101101 (2014),arXiv:1405.5303 [astro-ph.HE].[108] M. G. Aartsen et al. (IceCube), (2016),arXiv:1607.08006 [astro-ph.HE].[109] M. G. Aartsen et al. (IceCube),Astropart. Phys.66