Characterizing the local gamma-ray Universe via angular cross-correlations
Simone Ammazzalorso, Nicolao Fornengo, Shunsaku Horiuchi, Marco Regis
CCharacterizing the local gamma-ray Universe via angular cross-correlations
Simone Ammazzalorso , , Nicolao Fornengo , , Shunsaku Horiuchi , and Marco Regis , Dipartimento di Fisica, Universit`a di Torino, via P. Giuria 1, I–10125 Torino, Italy Istituto Nazionale di Fisica Nucleare, Sezione di Torino, via P. Giuria 1, I–10125 Torino, Italy and Center for Neutrino Physics, Department of Physics, Virginia Tech, Blacksburg, VA 24061, USA ∗ With a decade of γ -ray data from the Fermi -LAT telescope, we can now hope to answer how wellwe understand the local Universe at γ -ray frequencies. On the other hand, with γ -ray data alone it isnot possible to directly access the distance of the emission and to point out the origin of unresolvedsources. This obstacle can be overcome by cross-correlating the γ -ray data with catalogs of objectswith well-determined redshifts and positions. In this work, we cross-correlate Fermi -LAT skymapswith the 2MPZ catalog to study the local z < . γ -ray Universe, where about 10% of the totalunresolved γ -ray background is produced. We find the signal to be dominated by AGN emissions,while star forming galaxies provide a subdominant contribution. Possible hints for a particle DMsignal are discussed. Keywords: gamma rays: unresolved diffuse background; dark matter
I. INTRODUCTION
The extragalactic γ -ray background (EGB) is definedas the the γ -ray emission remaining after the subtrac-tion of all Galactic sources from the γ -ray sky. It shouldbe sourced by various classes of extragalactic γ -ray emit-ters, including the common star-forming galaxies, Ac-tive Galactic Nuclei such as blazars, and cascades ofhigh-energy particle propagation (for a recent review, seeRef. [1]). Exotic sources, such as dark matter annihila-tion or decay, can also contribute to this signal. In theera of the Fermi -LAT satellite, much has been revealedabout the origins of the EGB. Some ∼ γ -ray sources, dominantly blazars, have been resolved [2],which explain up to half of the EGB [3], and the numberwill almost double with the upcoming FL8Y point sourcecatalog. Removing these extragalactic point sources fromthe EGB leaves a residual, the so-called unresolved (orisotropic) γ -ray background (UGRB) [4], whose originsremain debated and is the focus of this analysis.The large numbers of EGB point sources detected haveenabled increasingly sophisticated predictions for theircontributions to the UGRB [5–9]. Often, these utilizeextrapolations of multi-wavelength observations to pre-dict the source behaviors in the faint unresolved end. Inparallel, a number of new and complementary techniqueshave been developed to study the UGRB in a more di-rect way. These uniquely exploit the sub-threshold infor-mation in the spatial distribution of γ -ray photons, andinclude the techniques of anisotropy [10–18], pixel statis-tics [19–24], and spatial cross-correlation with tracers oflarge-scale structure [25–42].Galaxies provide abundant opportunities that allowpowerful probes of the local large-scale structure of theUniverse. In recent works, Fermi -LAT data were cross-correlated with a variety of galaxy catalogs, including ∗ [email protected], [email protected], [email protected],[email protected] the SDSS-DR6 quasars, SDSS-DR8 main galaxies, SDSS-DR8 luminous red galaxies, SDSS-DR12 photo- z galax-ies, NVSS radiogalaxies, WI × SC galaxies, the 2MASSgalaxies, and the 2MPZ subsample of 2MASS galaxies.Positive correlations (at the level of 3–5 σ ) were detectedon angular scales of (cid:46) ◦ with all but the luminous redgalaxy catalogs [32–35], providing valuable informationon the sources behind the UGRB and constraints on darkmatter contributions. Tomographic analyses, wherebydepth (redshift) information is also utilized in unravel-ing the sources of the correlation signals [31], have beensuccessful in increasing the significance of the measuredcorrelations with some galaxy catalogs to ∼ σ [40].In this work, we perform new analyses of the cross-correlation that focus on disentangling the astrophysicaland exotic contributions to the UGRB at low redshift.Galaxy observables, e.g., in the B - and K -bands, pro-vide proxies for the amount of astrophysical activity anddark matter, respectively. Thus they can be used topredict astrophysical background and dark matter sig-nal strengths. In order to capture this information, weexploit the plethora of multi-wavelength data availableon galaxies and perform new position cross-correlationanalyses using galaxies divided into multiple quadrantsof astrophysical and dark matter signal expectations.We work with the 2MASS Photometric Redshift cata-log (2MPZ), which consists of cross-matching 2MASSXSC, WISE and SuperCOSMOS all-sky samples, whichprovide multi-wavelength data in 8 wavelengths (B, R,I, J, H, Ks, W1, W2) for over a million galaxies withdistribution peaked at z = 0 .
07. Simply put, one ex-pects dark matter to correlate most cleanly with massiveyet astrophysically inactive targets, and also in nearbygalaxies since competing astrophysical processes peak athigher redshifts. The fact that dark matter peaks at low- z stems from three competing effects: stronger clustering(namely, higher concentration for dark matter halos) as z decreases, higher average dark matter density as z in-creases (scaling as (1 + z ) ), and dilution of the observedradiation as z (i.e., distance) increases. The first and a r X i v : . [ a s t r o - ph . C O ] N ov the latter effects win over the second (see, e.g., [31]),and the different redshift distribution of the dark mat-ter signal compared to astrophysical backgrounds is oneof the most important features making cross-correlationanalyses relevant for constraining the particle dark mat-ter nature. Generically speaking, the method can probeweakly interacting massive particle (WIMP) dark matterwith annihilation cross section around the thermal value(depending on the mass and type of analysis [32]), asconfirmed also by the work presented here.This paper is organized as followed. In Section II, wedescribe the γ -ray data and galaxy catalogs used. In Sec-tion III, we describe our analysis procedure. In SectionIV we present our results and provide discussions for theorigins of the UGRB. Section V concludes. We providedetails of our various validation checks in Appendix A,and treatment of source modeling in Appendices B andC. Throughout, we adopt the Planck cosmology with pa-rameters from Ref. [43]. II. DATA
The datasets that we employ in our cross-correlationanalyses are (i) the first 9-years data release of γ -raysfrom Fermi -LAT, for which we consider a broad energyrange running from 630 MeV to 1 TeV, and (ii) the 2MPZgalaxy catalog. The data sets and data selection are de-scribed in the next subsections.
A. Fermi-LAT
Fermi -LAT is a γ -ray pair-conversion telescopelaunched in June 2008. It offers excellent capabilities toinvestigate the nature of the extra-galactic γ -ray back-ground, covering an energy range between 20 MeV and1 TeV with remarkable angular resolution ( ∼ . ◦ above10 GeV) and rejection of charged particles background.In this work, we use 108 months of data from Au-gust 4th 2008 to July 13th 2017 ( Fermi
Mission ElapsedTime: 239557417 s – 521597050 s). The photon countsand exposure maps are produced with the LAT Sci-ence Tools version v10r0p5[44]. We select the Pass8[45] ULTRACLEANVETO event class (correspondingto the
P8R2_ULTRACLEANVETO_V6 instrument responsefunctions (IRFs)), which is recommended for diffuseemission analysis since it has the lowest cosmic-ray con-tamination.We use both back- and front-converting events. Forphoton energies below 1.2 GeV, where photon statisticsis significantly larger than at higher energies but directionreconstruction is worse, we use photons belonging to theevent class PSF3, which refers to the best-quality quartilein the reconstructed photon direction (technically, thiscorresponds to event type 32). At higher energies, forwhich the direction reconstruction is inherently betterbut photon statistics declines, we extend the selection to the three best quality quartiles PSF1 + PSF2 + PSF3(event type 56). This choice allows us to have at the sametime a very good angular resolution and a good photonstatistics in the whole energy range of our analysis.The analyses are performed on photon intensity maps,obtained by dividing the count maps by the exposuremaps and the pixel area Ω pix = 4 π/N pix . We adopt aHEALPix pixelation format with resolution parameter N side = 1024, which corresponds to N pix = 12 , , ∼ . ◦ , similar to the best an-gular resolution of the gamma-ray data. Intensity mapsare produced in 100 energy bins, evenly spaced in loga-rithmic scale between 100 MeV and 1 TeV. The resultingintensity maps are then re-binned in larger energy binsfor the cross-correlation studies. After various tests, wedecided to limit the lowest photon energy to 630 MeV,largely determined by the poorer angular resolution be-low this scale. We perform all cross-correlation analysesin 11 energy bins, evenly spaced in logarithmic scale be-tween 630 MeV and 1 TeV, and projected in HEALPixmaps with N side = 1024.An example of γ -ray map obtained by including allphotons above 1 GeV is shown in the Fig. 1 (left).
1. Masking intensity maps
Although we do not expect a correlation between theGalactic γ -ray diffuse emission and the extra-galacticmatter distribution traced by the galaxy catalogs, wenevertheless need to exclude the very bright Galacticemission, especially along the Galactic plane (in orderto reduce the noise). We therefore perform a Galacticplane cut by masking galactic latitudes | b | < ◦ and,in addition, we further subtract the Galactic foregroundemission from the data maps. Resolved point sources arealso masked, in order to leave in the intensity maps onlythe components contributing to the UGRB.The point source masks are based on Fermi -LAT cat-alogs of resolved sources. We select the sources from theFL8Y catalog: this is a preliminary source list releasedby the
Fermi -LAT Collaboration which contains almostall the pre-identified sources of the 3FGL catalog aug-mented by new ones. It includes 5523 sources in the 100MeV – 1 TeV energy range [47]. For energies above 10GeV, we additionally mask sources from the 3FHL [46]catalog, which contains 1556 objects characterized in the10 GeV – 2 TeV energy range, in order to account forhard-spectrum sources that might be not contained inthe FL8Y catalog. The mask is built by taking into ac-count both the angular resolution of the detector in eachspecific energy bin and the brightness of the source to bemasked. Specifically, for each source we mask the pix-els inside a circle of radius R around its position definedthrough the following condition: F γ ∆ E exp (cid:18) − R θ E (cid:19) > F γ ∆ E, faintest FIG. 1. Left: All-sky
Fermi -LAT photon intensity map for photon energies above 1 GeV shown in Mollweide projection andsmoothed with a Gaussian beam of size σ = 0 . ◦ for illustration purposes. Right: Galaxy counts map of the full 2MPZ catalog,in Mollweide projection (the map has been downsized to N side = 128 for illustration purposes). where F γ ∆ E is integral flux of the source in a given en-ergy bin ∆ E , F γ ∆ E, faintest is the flux of the faintest sourcein the same energy bin, and θ ∆ E is the 68% contain-ment angle in that energy bin, as provided by the Fermi -LAT point-spread-function (PSF) analysis. The thresh-old condition based on 1/5 of the flux of the faintestsource corresponds approximately to the rms in the spe-cific energy bin (sources are detected with TS ≥ Fermi -LAT PSF with increasing energy allows to set energy-dependent masks which improve (i.e., become progres-sively less constraining) as energy grows, which is impor-tant since it coincides with where photon statistics be-come reduced. Sources that are marked as “extended” inthe FL8Y/3FHL catalog, are masked with the “extensionradius” provided in the
Fermi -LAT catalog. An exampleof the mask is shown in Fig. 2 (left) for the energy bin(1 . , .
3) GeV.Foreground removal is done by using the Galactic emis-sion model gll iem v06.fits of the
Fermi -LAT Collab-oration [48]. Foreground template maps are producedin the same 100 energy bins, evenly spaced in logarith-mic scale between 100 MeV and 1 TeV, and projectedin HEALPix maps with N side = 1024 as introduced forthe intensity maps. Each template map is assigned afree normalization (and added to a free constant, repre-senting the UGRB and cosmic-ray contamination) and aPoissonian likelihood fit is performed globally on all themasked intensity maps. Through this procedure, we ob-tained that all the best-fit normalization parameters areof the order of unity, supporting a successful descriptionof the foreground emission. The normalized foregroundtemplates are then re-binned into the 11 energy binsused for the cross-correlation analyses, and subtractedfrom the corresponding intensity maps. The robustness of foreground removal and choice of foreground model arediscussed in Appendix A 2.The masked intensity map in the energy bin (1 . − .
3) GeV, after subtraction of the Galactic foregroundemission, is shown in the right panel of Fig. 2.
B. Galaxy catalogs
For our analysis we employ the 2MPZ catalog [49],which has been built by cross-matching 2MASS XSC,WISE and SuperCOSMOS all-sky samples. The cata-log contains ∼ galaxies and their photometric red-shifts have been reconstructed via an artificial neural net-work approach. All the 8 magnitudes (B, R, I, J, H, Ks,W1, and W2) measured in SuperCOSMOS, 2MASS andWISE are present. In order to perform our measurementwe use the mask described in Ref. [50], which avoids sys-tematics due to Galactic dust contamination or misiden-tification that derives from high stellar number densities.Our goal is to decipher the composition of the UGRBat low-z. The different γ -ray emitters considered inthis work—dark matter (DM), star forming galaxies(SFG), blazars (BLZ) and misaligned active galactic nu-clei (mAGN)—can show different levels of correlationswith different subsamples of the 2MPZ catalog that tracedifferent properties of galaxies. In fact, different γ -raysources can have different redshift behaviors and can behosted by different types of galaxies. Therefore, in an at-tempt to enhance the sensitivity of the cross-correlationanalysis to the different γ -ray type of sources, in some ofour analyses we subdivide the galaxy catalog in severalsubsamples, as described in the following subsection. FIG. 2. Left: Instance of a masked intensity map: the plots refer to the energy bin (1 . , .
3) GeV. Gray pixels define the mask,which covers Galactic latitudes b < ◦ and point sources. The left and right panels show the map without and with Galacticforeground removal.
1. Galaxy subsets
The full 2MPZ represents our reference catalog,shown in Fig. 1 (right). In addition, we consider thefollowing subsamples: • – The 2MASS Redshift Survey (2MRS,[51]) contains all the 2MASS sources for whicha spectroscopic redshift is available. The cata-log counts 50k objects with a mean redshift of z = 0 .
03, thus representing a low- z subsample of2MPZ. Since the DM signal is peaked at low red-shift, the 2MPZ has the potential to be more sen-sitive to the DM γ -ray emission. • Redshift bins – We perform a redshift slicing of2MPZ subdividing the catalog into three samples( z < .
07, 0 . < z < . z > . • B-luminosity bins – The B-band luminosity isa reasonable tracer of star formation activity (see,e.g., Ref. [52]), and thus, would be expected to cor-relate also with cosmic-ray induced astrophysical γ -ray emission [53]. We thus split the full 2MPZcatalog into three bins of absolute B-luminosity,with each sub-catalog containing again one thirdof the total number of galaxies. • K-luminosity bins – The K-band luminosity ofgalaxies are correlated with the stellar mass of thegalaxy which can be correlated with the halo massby, e.g., abundance matching [54]. Therefore, weconsider it a tracer of the object mass and we definethree sub-catalogs by slicing the full 2MPZ cataloginto three bins of absolute K-luminosity again eachone containing one third of the total number ofgalaxies. • High K–Low B – Objects with high-K and low-Bluminosities should have high mass and low level ofstar formation activity. Therefore they can be con-sidered as ideal targets for DM searches, since theymight have a reduced correlation with astrophys-ical γ -ray sources (having the emission driven bystar formation activity), whilst an enhanced corre-lation with γ -ray emission induced by DM (whichis related to the mass). In order to perform thisinvestigation, we select 10k objects in the cornerof the plane of K vs B absolute luminosity in the2MRS catalog (since the DM signal is peaked atlow- z ). We will report the results about this sam-ple only when focusing on the DM interpretationin Section IV C.Fig. 3 shows the redshift distributions of the full 2MPZcompared to those of the subsample catalogs. 2MRS isthe catalog peaking at the lowest redshift. The subsam-ples of the mid bins in both K and B luminosity havea redshift distribution close to the full 2MPZ, while thelow/high bins are peaked at lower/higher z .In addition to the subsamples listed above, we furtherdefine two selections of sources that aim at identifyingspecifically mAGN and BLZ in the 2MPZ catalog. Thisidentification will be useful to model the cross-correlationangular power spectrum of mAGN and BLZ, as describedin Section IV.Blazars are identified by cross-matching 2MPZ withthe WIBRALS catalog [55]. The latter is composed ofradio-loud WISE sources detected in all four WISE fil-ters, whose mid-infrared colors match typical colors ofconfirmed γ -ray emitting blazars. We select mAGNs bycross-matching 2MPZ with the AGN sample found inRef. [56]. The authors considered WISE and 2MASSdata and defined a statistical discriminator by compar-ing the measured infrared colors, producing a completesample of AGNs. This subset contains ∼ objects and FIG. 3. Redshift distributions of the different galaxy subsamples used in the cross-correlation analyses. Each distribution isnormalized to the total number of galaxies of its corresponding subsample. Low/Mid/High-B (K) refers to galaxies selectedaccording to their B (or K) luminosity. we remove blazars obtained from the WIBRALS catalog.
III. MEASUREMENTS
The cross-correlation APS is defined as: C ( ij ) (cid:96) = 12 (cid:96) + 1 (cid:88) m a ( i ) (cid:63)(cid:96)m a ( j ) (cid:96)m (2)where: a ( i ) (cid:96)m = (cid:90) d(cid:126)n δI ( i ) ( (cid:126)n ) Y (cid:96)m ( (cid:126)n ) (3)are the coefficients of the expansion of the fluctuations δI i ( (cid:126)n ) of the field I i ( (cid:126)n ) in terms of spherical harmon-ics Y (cid:96)m ( (cid:126)n ). In our case i and j correspond to the γ -ray and galaxy map fields. We determine the APS withPolSpice[57], a public code that computes both the two-point angular cross-correlation function in real space and the APS. PolSpice is based on the fast spherical harmonictransforms allowed by isolatitude pixelisations and it cor-rects for the effects introduced by masking following theapproach of Ref. [58]. PolSpice also provides an estimatefor the covariance matrix of the measurement, that willbe used for the statistical analysis discussed in the fol-lowing Sections.Before computing the APS, we remove the monopoleand dipole contributions from the input maps by ap-plying the HEALPix routine remove dipole , in orderto mitigate a possible leakage of these (large) terms tohigher multipoles (an effect due to multipole-mixing in-troduced by the masks).The finite angular resolution of the Fermi -LAT instru-ment suppresses the angular power spectrum at high mul-tipoles (the angular resolution of the galaxy surveys issignificantly better than the
Fermi -LAT one and the as-sociated suppression would show up only at higher multi-poles, in a range not considered in our analysis). In order
FIG. 4. Measured APS between the γ -ray map in the (1 . , . (cid:96) = 40 whilethe upper bound (cid:96) max is determined from the beam windowfunction and therefore depends depends on energy. All themeasured APS can be retrieved at this link [59]. to take this suppression into account, we correct the C (cid:96) with the beam window function: W (cid:96) ( E ) = 2 π (cid:90) − d cos θ P (cid:96) (cos θ ) PSF( θ, E ) (4)where P (cid:96) are the Legendre polynomials and PSF( θ, E )denotes the Fermi -LAT point spread function for the spe-cific IRF and energy, as provided by the
Fermi
ScienceTools. The energy-dependent beam window function isaveraged in each energy bin in accordance to the UGRBenergy spectrum E − α , where the spectral index is takenat α = 2 . (cid:104) W k(cid:96) (cid:105) = (cid:82) E max ,k E min ,k W (cid:96) ( E ) E − α dE (cid:82) E max ,k E min ,k E − α dE (5)The measured APS in the k -th energy bin is then de-fined as: C k(cid:96) = C k(cid:96), raw (cid:104) W k(cid:96) (cid:105) W pix , (6)where C k(cid:96), raw is the raw APS obtained from PolSpice inthe k -th energy bin and W pix is the pixel window functionassociated to the HEALPix pixeling.For the analyses discussed in the next Section, we re-bin the measured APS in 15 evenly-spaced logarithmicmultipole bins from 10 to 1000. Since the low multipoles C k(cid:96) can be affected by large-scale effects due to an imper-fect Galactic foreground removal and at large multipoles C k(cid:96) by an imperfect PSF correction, especially when the Bin E min [GeV] E max [GeV] (cid:96) min (cid:96) max beam window function starts deviating significantly from1, we must identify a suitable multipole range over whichwe perform our analyses: the lower limit is conservativelyset to (cid:96) min = 40; the upper limit (cid:96) max is defined from thecondition that the beam window function does not dropbelow a threshold corresponding approximately to the68% containment of the PSF in the specific k -th energybin: (cid:104) W k(cid:96) max (cid:105) = 0 . , (7)or l max = 1000, whichever is smaller. This conditionmakes (cid:96) max dependent on energy. The lower and upperbound of the multipole bins for each energy bin are shownin Table I.Fig. 4 shows an example of the measured APS, in the(1 . , .
3) GeV energy bin, for the cross-correlation withthe whole 2MPZ catalogue. The plot also shows the mul-tipole range ( l min , l max ) for this energy bin. Error barsare large at low multipoles because of cosmic variance,mask deconvolution and noise from Galactic foreground.They start becoming large also at multipoles above a fewhundreds because of the size of the Fermi -LAT PSF (andfinite statistics).All the measured APS can be retrieved at this link [59].
A. Amplitude and significance of the correlation
In order to provide a model-independent estimate ofthe amplitude and significance of the measured cross-correlations, we fit the APS in each energy bin with aterm which is multipole-independent (i.e., a constant)that we call C kp . This can be considered as the simplestmodel (i.e., a Poisson noise term) and provides an esti-mate of the amplitude which is similar to performing theaverage of the APS over the multipole range of interest.A more refined treatment, which involves modeling of the γ -ray components in the unresolved sky, is presented inthe next section. E [GeV] C p × E . / ∆ E [ G e V ( c m − s − s r − ) s r ] × − E [GeV] C p × E . / ∆ E [ G e V ( c m − s − s r − ) s r ] × − E [GeV] C p × E . / ∆ E [ G e V ( c m − s − s r − ) s r ] × − E [GeV] C p × E . / ∆ E [ G e V ( c m − s − s r − ) s r ] × − × FIG. 5. Multipole-independent APS C kp as a function of the energy, for the different galaxy subsamples considered in this work. Fig. 5 shows C kp as a function of the energy bin for thecross-correlation of the γ -ray flux maps with each of ourgalaxy subsamples. The plot indicates the presence ofa correlation signal between the galaxy distribution and γ -rays for all subsamples. In fact, the C kp are systemat-ically positive (i.e., they do not fluctuate around zero)and deviate from a null signal. To assess the significanceof the measurements, we compare the χ of a null signalwith the χ obtained from the C kp fit.We adopt a χ estimator defined as: χ = (cid:88) k =1 (cid:88) ∆ (cid:96), ∆ (cid:96) (cid:48) (8)( C k, mod∆ (cid:96) − C k, exp∆ (cid:96) )Γ − (cid:96), ∆ (cid:96) (cid:48) ,k ( C k, mod∆ (cid:96) (cid:48) − C k, exp∆ (cid:96) (cid:48) ) , where C k, exp∆ (cid:96) is the measured APS in the energy bin k andmultipole bin ∆ (cid:96) , C k, mod∆ (cid:96) is the APS model and Γ ∆ (cid:96), ∆ (cid:96) (cid:48) ,k is the covariance matrix, obtained from the PolSpice co- variance through multipole re-binning. We neglect thecovariance between different energy bins since the mainsource of error comes from the Poisson noise of the γ -raymaps, something which exhibits no correlation amongdifferent energy bins. Eq. 8 will be adopted throughoutthe paper for model comparison, including the analysis interms of γ -ray modeling as discussed in the next sections.We define a χ difference ∆ χ = χ − χ Cp , where χ is the null signal obtained from Eq. 8 by using C k ∆ (cid:96), mod = 0, and χ C p is obtained using C k ∆ (cid:96), mod = C kp .Table II shows the results for the different subsamples:for each case, ∆ χ > Subset χ χ C p ∆ χ z z χ results for the no-signal case and a multipole-independent C kp . The total num-ber of data-points considered for each sample (including en-ergy and multipole bins) is 114. the null hypothesis.Since in our correlation measurement we employ mapsof the integrated γ -ray flux, we expect the energy spec-trum to follow the integrated energy spectrum of theUGRB, namely I UGRB = (cid:82) ∆ E dE dI UGRB /dE . By multi-plying the vertical axis in Fig. 5 by E / ∆ E , we show (ap-proximately) the differential energy spectrum of the γ -rayemission responsible for the correlation signal, rescaledby E − . The statistical significance is not enough to de-rive firm conclusions on the energy dependence, but dif-ferent subsets seem to indicate a γ -ray population witha spectral index close to − γ -ray emit-ters, blazars typically show a hard spectrum (with indexabout − IV. RESULTS AND INTERPRETATIONA. Models
We model the clustering of galaxies in the samples pre-sented in Sec. II B 1 and of γ -ray emitters mentioned inthe introduction by making use of the halo model ap-proach. Galaxies are assumed to follow the matter powerspectrum with matter distributed in halos, and with thenumber of galaxies per halo defined by the so-called halooccupation distribution (HOD). The latter has been de-rived by fitting the auto-correlation APS of galaxies, asdescribed in the Appendix.The cross-correlation APS with γ -ray sources is com-puted as in Ref. [60] with two differences. The first in-volves the combination of flat spectrum radio quasars andBL Lacs into a single (effective) blazar class, as done inRef. [61]. The contribution to the UGRB of the γ -rayemitters considered in this work is shown in Fig. 6. Note E [GeV] -9 -8 -7 -6 E I [ G e V c m - s - s r - ] BLZ t e m p l a t e f itti ng z < 0.2SFGmAGN UG RB fr o m z < 4ann. DM FIG. 6. Energy spectrum of the UGRB as determined fromthe method described in Section II A (gray band) togetherwith the predicted contribution for the reference model ofthe different γ -ray emitters considered in this work. Withdashed lines, we show also their contribution to the UGRBfrom z < . from the figure that we expect the γ -ray emission we an-alyze in this paper (which is produced at z < .
2) toamount to a small fraction of the total UGRB, roughlyaround 10% at low energy. The three classes of emittersprovide comparable contributions, within a factor O (1).The second difference is related to the modeling of theshot-noise term. This contribution to the APS is the partat zero angular separation (i.e., (cid:96) -independent) of the 1-halo term. As recognized in Ref. [60], its modeling canbe very delicate. Here, we do not attempt to include itin our halo-modeling, whilst we follow two data-drivenapproaches. In the simplest approach, we just fit theshot-noise contribution (which is a constant term for eachAPS) by assuming a power-law energy dependence: dC ( j ) p dE = N j E − α j , (9)where the index j labels the galaxy sample. In this way,we introduce 2 parameters (normalization and power-lawindex) for each galaxy sample. This approach will becalled the “ free C p ” fit.The second approach, which is our reference one (called reference ), uses the fact that the shot-noise is given bythe average of the γ -ray flux of all the galaxies N gal ofa given sample j (let us remember that index k denotesthe energy bin): C ( jk ) p = 1 N j gal N j gal (cid:88) i =0 F γ ∆ E k ,j . (10)We adopt some empirical relations to predict F γ ∆ E k ,j fromthe optical/infrared magnitude of the each galaxy in thecatalog, as described in the Appendix. Furthermore,we need to identify blazars and misaligned AGNs in the2MPZ catalog, for which we use the procedure outlinedin Sec. II B 1. With these two ingredients, we are able toestimate the shot-noise contribution.The signal associated to annihilating DM is computed,again following the halo model approach, as described inRef. [33], with the “boost-factor” taken from Ref. [ ? ].The contribution depends on two parameters, the parti-cle DM mass M χ and the annihilation cross section σv .We will consider four different DM models, referring tofour specific annihilation final states endowed with differ-ent spectra and representative for a typical WIMP DM: b ¯ b , τ + τ − , W + W − , and µ + µ − .Summarizing, we fit the cross-correlation APS of thegalaxy samples presented in Sec. II B 1 and the Fermi -LAT γ -rays intensity maps with two approaches: • Reference model: C ( jk ) (cid:96), mod = C ( jk ) (cid:96), DM ( M χ , σv )+ N SFG × ( C ( jk ) (cid:96), SFG + C ( jk ) p, SFG )+ N BLZ × ( C ( jk ) (cid:96), BLZ + C ( jk ) p, BLZ )+ N mAGN × ( C ( jk ) (cid:96), mAGN + C ( jk ) p, mAGN ) . In this approach, the total number of free param-eters is 5, i.e., 3 normalizations ( N SFG , N BLZ and N mAGN ) for the astrophysical contributions and 2terms for the annihilating DM contribution ( M χ and σv ). The annihilation rate will be expressed interms of the “thermal” (or “natural scale”) value (cid:104) σv (cid:105) th = 3 × − cm / s by trading it for a dimen-sionless parameter N DM = σv/ (cid:104) σ a v (cid:105) th . • Free C p model: C ( jk ) (cid:96), mod = C ( jk ) (cid:96), DM ( M χ , σv ) + N SFG × C ( jk ) (cid:96), SFG + N BLZ × C ( jk ) (cid:96), BLZ + N mAGN × C ( jk ) (cid:96), mAGN + C ( jk ) p , where the last term C ( jk ) p = (cid:82) E k max E k min dE dC ( j ) p /dE ,with E k min and E k max being the energy boundariesof the k -th energy bin. With respect to the previouscase, this model adds 2 parameters for each sample j , associated with the C p term (see Eq. 9). B. Statistical analysis
Our fit is performed with the Monte Carlo parameterestimation code CosmoSIS [62]. Since the order of magni- tude of each parameter is unknown, we use a Metropolis-Hastings sampler with a flat prior in log-scale for eachparameter.The galaxy subsamples listed in Sec. II B 1 are analyzedseparately. For the cases involving three bins (redshift,B-luminosity and K-luminosity), we fit simultaneouslythe APS of the different bins, which are independent fromeach other (since the galaxy subsamples are not overlap-ping). For these samples, the number of parameters inthe fit is 5 (11) in the reference ( free C p ) model. For allthe other samples the number of parameter is 5 (7) inthe reference ( free C p ) model.As an example of the outcome, in Fig. 7, we show thetriangular plot obtained by fitting the 2MPZ split intoredshift bins. The vertical dashed and solid red (green)lines denote the 68% and the 95% CL upper (lower) limitsfound with the profile likelihood, respectively. In the2D panels, the 68% regions are identified in cyan whilethe 95% regions are in dark blue. In this example, theonly normalization which is not compatible with zero (at1 σ level) is for the mAGN population, see last panel ofsecond row.All the triangular plots for the various cases are avail-able at this link [59]. In the following, for the sake ofbrevity, we will focus our discussion on the 1D profilelikelihood distributions, except in the case of the DMparameters, for which we will discuss also the 2D planeshowing the bounds on the particle DM parameters inthe canonical annihilation rate vs DM mass space.Fig. 8 summarizes the results on the normalization pa-rameters of the astrophysical γ -ray sources for the ref-erence analysis. The upper three panels shows the 1Dlikelihood distributions for SFG, BLZ and mAGN ob-tained by organizing the galaxy data into the three dif-ferent subsamples that differentiate the galaxies in termsof redshift, K luminosity and B luminosity. The lowerthree panels show the results for SFG, BLZ and mAGNwhen the full 2MPZ catalog (blue) or the low-redshift2MRS catalog (yellow) are used.The corresponding DM results for the reference anal-ysis are shown in Fig. 9, for DM annihilating in the b ¯ b channel. The upper panels show the likelihood distri-butions for the annihilation rate for the different galaxysubsamples, while the lower panels show the correspond-ing 95% CL bounds on the annihilation rate as a func-tion of the DM mass for the same annihilation channel.The bounds for all the four annihilation channels consid-ered in this work ( b ¯ b , τ + τ − , W + W − , and µ + µ − ) and forthe analysis performed combining the three z -bins of the2MPZ catalog are shown in Fig. 10. This figure can beconsidered as the summary plot for what concerns thebounds on WIMP DM derived in this work.As a further investigation of the DM case, Fig. 11 con-siders galaxy samples for which the cross-correlation with γ -rays is expected to be enhanced, i.e., the low-redshift2MRS sample and its combination with the High K–LowB subsample of the 2MPZ catalog. Again, the left andright panels show the likelihood distribution for the anni-0 FIG. 7. Fit results for the 2MPZ redshift slicing subset for the reference analysis. All parameters are shown in log-scale. Thevertical dashed and solid red (green) lines denote the 68% and the 95% CL upper (lower) limits obtained from the profilelikelihood, respectively. In the 2D plots, the 68% regions are identified in cyan while the 95% regions are in dark blue. The 1Dprofile likelihood distributions on the diagonal are individually normalized to unity. hilation rate and the 95% CL bounds in the annihilationrate vs mass plane.Table III reports the best-fit values and the 68% upperand lower bounds (whenever present) for the astrophysi-cal and DM parameters, for the different galaxy samples.Discussion and interpretation of the results are presentedin the next section.For the free C p analysis, the results are shown inFigs. 12 and 13, that mirror the information in Figs. 8and 9, respectively. Table IV lists the best-fit values andthe 68% upper and lower bounds (whenever present) forthe astrophysical and DM parameters, for the different galaxy samples. In Table V we show the best fit resultsfor the C p normalizations and power-law indexes.Finally, the statistical significance of the reference and free C p models as compared to the null hypothesis ofabsence of signal are shown in Table VI in terms of the χ differences. C. Interpretation of the results
In this Section we discuss the interpretations of theresults presented in the previous section and conclude1 -2 -1 N SFG li k e li h oo d z binsK binsB bins10 -2 -1 N SFG li k e li h oo d -2 -1 N BLZ li k e li h oo d z binsK binsB bins10 -2 -1 N BLZ li k e li h oo d -2 -1 N mAGN li k e li h oo d z binsK binsB bins10 -2 -1 N mAGN li k e li h oo d FIG. 8. Profile likelihood distributions for the normalization parameters of the astrophysical γ -rays components, for the reference analysis. The upper panels show the results obtained for the different subsamples of the 2MPZ catalog. The lowerpanels show the results for the full 2MPZ and for the low-redshift 2MRS catalogs. The vertical solid (dashed) lines indicatethe 68% upper (lower) limits (whenever present in the plots).Sample N mAGN N SFG N BLZ N DM BF low up BF low up BF low up BF low up2MPZ (full) 0.02 - 3.24 0.76 0.13 1.29 1.95 - 3.47 190.55 - 575.442MRS (full) 0.35 - 0.74 0.06 - 67.61 0.02 - 3.16 7.59 0.56 25.70 z bins 2.45 0.85 3.47 0.07 - 0.23 0.02 - 1.95 181.97 - 478.63B bins 1.45 - 2.95 0.15 - 0.36 0.66 - 3.16 165.96 - 416.87K bins 2.09 0.27 3.16 0.14 - 0.31 0.03 - 2.14 165.96 - 426.58TABLE III. Best fit and 68% C.L. interval of the various parameters in the fit for the reference case. When the lower boundis not reported, it means that it is compatible with zero at the quoted CL. the consequences for the extragalactic γ -ray populationsconsidered in our modeling.
1. Star forming galaxies
Star-forming galaxies are poorly constrained by ouranalysis. We find no relevant peak in the 1D likelihooddistributions (i.e., in the left panels of Figs. 8 and 12) andupper bounds on its normalization are around 2–3 timesthe reference model. This implies that in order to providea significant contribution to the cross correlation APS measurements derived in this work, SFG would overshootthe total UGRB intensity level shown in Fig. 6. In otherwords, we found that SFG are a subdominant componentof the UGRB at low redshift.
2. Misaligned AGNs
Misaligned AGNs appear to be the population whichcan explain the bulk of the measured signal. They arethe only population that is singled out with statisticalconfidence in several datasets. It is interesting to note2 -2 -1 σv/ › σv fi th li k e li h oo d z binsK binsB bins 10 -2 -1 σv/ › σv fi th li k e li h oo d m DM [GeV] σ v / › σ v fi t h z binsK binsB bins 10 m DM [GeV] σ v / › σ v fi t h FIG. 9. Results for the DM case obtained with the the reference analysis and b ¯ b annihilation channel (with “boost-factor” fromRef. [ ? ]). The upper panels show the profile likelihood distribution for the annihilation rate. The lower panels show the 95%CL upper bounds for the annihilation rate vs the DM mass. The two panels in the first column refer to the analyses performedon different organization of the galaxy samples (redshift, K luminosity and B luminosity). The two panels in the second columnrefer to the analyses on the full 2MPZ catalog and on the low-redshift 2MRS catalog. the power of the “tomographic” approach. As shown inthe bottom right panel of Fig. 8, when considering thefull 2MPZ sample, no peak is present in the likelihooddistribution. However, the evidence appears when con-sidering the z and K -luminosity bins (top right panel).A preference for mAGNs is also found in the free C p case,as shown in Fig. 12. Note that in this case, the mAGNnormalization exhibits a lower limit already in the full 2MPZ sample.We see that the 2MRS catalog seems to set an upperbound for the mAGN normalization that excludes thebest-fits obtained with all the other samples (see bot-tom right panel of Fig. 8). On the other hand, this doesnot happen in the free C p case (bottom right panel ofFig. 12), where instead the upper limit is consistent withthe normalizations estimated from other galaxy samples.3 -2 -1 σv/ › σv fi th li k e li h oo d b ¯ b channel µ ¯ µ channel τ ¯ τ channel W ¯ W channel 10 m DM [GeV] σ v / › σ v fi t h b ¯ b channel µ ¯ µ channel τ ¯ τ channel W ¯ W channel FIG. 10. Left panel: profile likelihood for the DM annihilation rate for the four annihilation channels considered in thisanalysis. Right panel: 95% CL upper bounds on the DM annihilation cross-section as a function of the DM mass for the sameannihilation channels. The plot refers to the reference analysis performed combining the three z -bins of the 2MPZ catalog.Sample N mAGN N SFG N BLZ N DM BF low up BF low up BF low up BF low up2MPZ (full) 3.02 1.17 4.17 0.01 - 3.09 0.02 - 1.15 120.23 - 588.842MRS (full) 0.49 - 5.75 0.09 - 58.88 0.07 - 2.82 11.22 - 61.66 z bins 2.63 1.32 3.89 0.01 - 2.19 0.04 - 1.07 77.62 0.28 407.386B bins 2.45 1.05 3.39 0.05 - 1.91 0.02 - 0.79 25.70 - 371.54K bins 2.24 0.69 3.31 0.07 - 2.88 0.02 - 1.17 19.05 2.00 331.13TABLE IV. Best fit and 68% C.L. interval of the various parameters in the fit for the free C p case. When the lower bound isnot reported, it means that it is compatible with zero at the quoted CL. We remind that in the latter case, the C p are allowed tovary and are determined by the fit. These facts pointtoward a possible overestimate of the mAGN shot-noiseat very low- z (i.e., in the range covered by 2MRS ) inthe reference model. In fact, in this case the shot noisehas been derived from relations which show significantscatter (see Appendix): when applied to a very smallvolume like in the case of 2MRS, the shot-noise estimatemight be not very accurate. A dedicated analysis focus-ing on the low-redshift 2MRS catalog will be the subjectof future work.
3. Blazars
In the analysis of blazars, we can appreciate again howthe tomographic approach tightens the bounds in Fig. 8,pushing the normalization to lower values when goingfrom the lower to the upper panel. Taken at face value,the results of the reference case (reported also in Ta- ble III) would indicate that BLZ are constrained to be asubdominant component of the total UGRB (see Fig. 6where the BLZ component should be rescaled by a factor N BLZ ). Blazars are so constrained essentially becauseof their large shot-noise term that contributes in a non-negligible way to the cross APS signal we measure.On the other hand, in the free C p model, the boundsbecome weaker, actually suggesting the opposite picture,namely that BLZ are a subdominant component of thecross-correlation we measure and of the UGRB at lowredshift. Indeed, they need N BLZ to be quite larger than1 to become a relevant component in our measurement.With such values of N BLZ , blazars can provide the bulkof the total UGRB emission (a picture similar to the SFGcase).It is clear that to distinguish between the two inter-pretations, the model of the shot-noise term is crucial.Physically, this is because we have already observed andcataloged a significant fraction of the closest γ -ray emit-ting blazars, and thus the possible cross-correlation signal4 -2 -1 σv/ › σv fi th li k e li h oo d m DM [GeV] σ v / › σ v fi t h FIG. 11. Results for the DM case using galaxy catalog samples expected to be more sensitive to the DM γ -ray signal, i.e., thelow redshift catalog 2MRS (yellow line) and its combination with the High K–Low B subsample (blue line). The results referto the reference analysis. Left: profile likelihood distributions for the DM annihilation rate; the vertical solid (dashed) linesindicate the 68% upper (lower) limits (whenever present in the plots). Right: 95% upper bound on the annihilation rate vs theDM mass. Sample N α N α N α . × − . × − z bins 3 . × − − .
96 2 .
51 10 − − .
29 9 . × − . × − .
57 10 − − .
04 1 . × − . × − . × − . × − free C p case. for the unresolved part is generated by a relatively smallnumber of sources, providing a large shot-noise term. Asmentioned above, the model of the latter depends on pre-dicting the γ -ray luminosity of blazars from their IR lumi-nosity. If the relation obtained in Ref. [63] extends to theunresolved regime, the conclusion of the reference case islikely to hold. On the other hand, a lower γ -ray lumi-nosity for the corresponding IR luminosity would pointtowards the outcome of the free C p model. A future ded-icated cross-matching analysis of the Fermi -LAT FL8Ysource list with multi-wavelength data could help in clar-ifying the picture.
4. Dark Matter
We now discuss the implications for particle DM.Figs. 9 and 13 show the results for the reference and free C p methods, respectively, for the b ¯ b annihilation chan-nel. The different samples and methods provide compat- ible constraints, all excluding annihilation rates higherthan (about) the “thermal” rate for DM mass of 10 GeVand then increasing with a nearly linear trend for highermasses. The cases of τ + τ − and W + W − final states leadto similar results, whilst constraints for DM annihilatinginto µ + µ − are about one order of magnitude weaker, ascan be seen in Fig. 10.The 1D distributions of the annihilation rate reportedin Figs. 9 and 13 show a small peak (in all samples). Thepeak becomes enhanced and shifts to lower annihilationrates when the low-redshift 2MRS catalog is used, bothin the reference and free C p analyses. In order to un-derstand if it just a statistical fluctuation or it might berather a hint for a DM contribution, we deepen the inves-tigation by considering a further subsamples, tailored tothe expected behavior of a DM signal. Ideally, in order toemphasize the DM γ -ray contribution over the astrophys-ical ones, we need to focus the cross-correlation analysison a catalog with galaxies at low- z [64], for galaxies withhigh-mass and with the lowest possible level of star for-5 -2 -1 N SFG li k e li h oo d z binsK binsB bins10 -2 -1 N SFG li k e li h oo d -2 -1 N BLZ li k e li h oo d z binsK binsB bins10 -2 -1 N BLZ li k e li h oo d -2 -1 N mAGN li k e li h oo d z binsK binsB bins10 -2 -1 N mAGN li k e li h oo d FIG. 12. The same as Fig. 8, but for the free C p analysis.Subset χ χ χ C p ∆ χ − ref ∆ χ − free C p z bins 253.07 229.24 230.36 23.83 22.71B bins 263.21 233.82 236.18 29.39 27.03K bins 270.44 241.32 242.98 29.12 27.46TABLE VI. Comparison of the best-fit χ results for the absence of signal ( χ ), the reference analysis ( χ ), the free C p analysis ( χ C p ) and the relative differences of the two latter with the no-signal case. mation and AGN activity. To these ends, we select 10kgalaxies from the 2MRS catalog (low- z ), in the corner ofhigh K-luminosity (which corresponds to high mass) andlow B-luminosity (which corresponds to low star forma-tion rate). The results are shown in Fig. 11, where wefocus on the b ¯ b case and the reference analysis, for thesake of brevity.Interestingly, the peak in the likelihood distributionof the annihilation rate slightly increases in height andmoves its position towards lower values of the normaliza-tion parameter. Note that it is the most pronounced peakin the likelihoods of the DM annihilation rate amongthe different samples. Even though the statistical signifi-cance is too low to speculate on the possible presence of aDM contribution, we highlight that in Figs. 9, 11 and 13 the significance tends to increase when considering sam-ples with objects at lower redshift and with higher halomasses, just as expected for a dark matter origin. Thisstimulates to further pursue the particle DM quest ex-ploiting the cross-correlation approach with future dataand dedicated studies.The best-fit for the 2MRS/High-K/Low-B analysis oc-curs at M χ = 37 GeV and σv = 4 × (cid:104) σv (cid:105) th , thereforein slight excess over the “natural” scale. However, aswidely discussed in the literature (see, e.g., Ref. [33]), thenormalization of the DM signal can significantly vary de-pending on the modeling of the so-called “boost-factor”provided by the substructure contribution (because ofunknowns in the definition of the minimum halo mass,subhalo mass function and subhalo concentration param-6 -2 -1 σv/ › σv fi th li k e li h oo d z binsK binsB bins 10 -2 -1 σv/ › σv fi th li k e li h oo d m DM [GeV] σ v / › σ v fi t h z binsK binsB bins 10 m DM [GeV] σ v / › σ v fi t h FIG. 13. Same as Fig. 9 but for the free C p analysis. eter). Therefore, the normalization of the DM signal canbe easily modified by a factor of a few by introducing asubstructure modeling different from the one consideredhere. V. CONCLUSIONS
In this work, we have made an attempt to characterizethe unresolved γ -ray emission of the Local Universe. Tothis aim, we employed Fermi -LAT skymaps with detectedsources being masked and performed the measurement of their angular cross-correlation with the 2MPZ catalog.The latter contains about one million of galaxies with amedian redshift of 0.07. The cosmological volume probedby 2MPZ powers only about 10% of the total unresolved γ -ray background. Despite this small fraction, the tech-nique adopted here enables us to study the compositionof such emissions.The null hypothesis, i.e., the absence of correlation be-tween the two datasets, is excluded at a statistical confi-dence larger than 99 . γ -ray luminosity for faint blazarsis crucial.Finally, we evaluated the possible contribution of aparticle DM signal. The 95% C.L. bounds on the DM an-nihilation rates reach close to the “thermal” rate for DMmass of 10 GeV for b ¯ b , τ + τ − and W + W − annihilationchannels (while an order of magnitude weaker bound isfound for µ + µ − ) and then increasing with a nearly lineartrend for higher masses. Interestingly, when consideringsamples where the DM evidence is expected to increase(namely, correlation with objects at low- z , with high-mass, and low level of star formation), we see a slightlymore pronounced peak in the DM likelihood for the DMcontribution. Currently, the statistical significance of thiseffect is low, and it prevents us from deriving any firmconclusions on the presence of a DM signal. Neverthe-less, this result motivates to deepen the investigation ofcross-correlations between suitable galaxy catalogs (espe-cially low redshift ones, like 2MRS) and multiwavelengthobservations, to probe the potential contribution of DM. ACKNOWLEDGEMENTS
We would like to thank M. Bilicki, A. Cuoco, F. Mas-saro and M. Negro for discussions. This work is sup-ported by the “Departments of Excellence 2018 - 2022”Grant awarded by the Italian Ministry of Education, Uni-versity and Research (MIUR) (L. 232/2016). NF is sup-ported by the research grant “The Anisotropic Dark Uni-verse” Number CSTO161409, funded under the programCSP-UNITO “Research for the Territory 2016” by Com-pagnia di Sanpaolo and University of Torino. The work ofSH is supported by the U.S. Department of Energy underAward No. DE-SC0018327. MR acknowledges supportby the Excellent Young PI Grant: “The Particle Dark-matter Quest in the Extragalactic Sky” funded by theUniversity of Torino and Compagnia di San Paolo and by“Deciphering the high-energy sky via cross correlation”funded by Accordo Attuativo ASI-INAF n. 2017-14-H.0.SA, NF and MR acknowledge support from the project “Theoretical Astroparticle Physics (TAsP)” funded bythe INFN.
Appendix A: Validation and cross-checks
In this appendix, we present a series of tests performedin order to validate our analysis.
1. Theoretical estimation of the error
It is possible to provide a theoretical estimation of theerror δC l associated to the cross-correlation in each mul-tipole bin, assuming gaussian statistics: δC l = (cid:115) ( C ( γ, gal) l ) + C ( γ,γ ) l C (gal , gal) l (2 l + 1) f sky ∆ l , (A1)where f sky is the fraction of unmasked sky, ∆ l is themultipole bin size, and C ( γ, gal) l , C ( γ,γ ) l and C (gal , gal) l arethe cross-correlation, auto-correlation (including noise)of the γ -ray data and auto-correlation (including noise) ofthe galaxies, respectively. In the top left panel of Fig. 14we show the errors on the cross-correlation of the γ -raydata in the energy interval from 1 to 10 GeV with thewhole 2MPZ catalog. We find that the theoretical er-ror of Eq. A1 is similar to and typically slightly smallerthan the one estimated by PolSpice, which we then usethroughout our analyses.
2. Foreground dependence
In order to assess the independence of our analysisfrom Galactic γ -ray foreground subtraction, we performa cross-correlation analysis of the combined energy binsfrom 1 to 10 GeV (which contains about 60% of the totalphoton counts) using γ -ray data that have been cleanedup by the diffuse Galactic emission (as explained in sec-tion II A) and compare those results with a correspondinganalysis performed on the same data without foregroundremoval. The top central panel of Fig. 14 shows thatthe APS derived with and without foreground removalare in excellent agreement, which confirms the hypoth-esis that the cross-correlation of the γ -ray flux with ex-tragalactic tracers of the γ -rays emitters is not affectedby the Galactic γ -rays foreground. This suggests thatthe cross-correlation results are not strongly dependenton the specific foreground model used in foreground re-moval. To further confirm this point, we compute C kp asa function of the energy bin as in Fig. 5 employing differ-ent foreground models in the analysis. In addition to ourreference case, we introduce models A, B and C presentedin [4]. The differences in our results among the four casesare negligible, as can be seen in Fig. 14 (top-right) for theexample of the full 2MPZ catalog. Note however from8the central panel that the presence of an un-subtractedforeground emission results in a noisier dataset that isreflected in larger uncertainties, especially at lower mul-tipoles.
3. Lower multipoles removal
While in a full-sky APS analysis the different multi-poles are independent, the presence of masks results incouplings between different multipoles. PolSpice correctsfor this effect, but a residual contamination is still poten-tially present. Moreover, with the monopole being largelydominant (e.g., for γ -rays the total average intensity ismuch larger than its fluctuations), even a small residualcoupling can bias the measurement at higher multipoles.In our analysis we remove the monopole and the dipolebefore performing the APS measurement, and we con-sider the APS only for multipoles larger than 40, as dis-cusses in section III. We nevertheless performed a checkto verify that our measurement is not affected by lowermultipoles, by comparing the APS results in the mul-tipole window (cid:96) ≥
40 when we removed from the mapthe contribution of multipoles up to (cid:96) ≤ (cid:96) ≤ i) com-pute the spherical harmonic decomposition coefficients a lm of the maps with the HEALPix routine anafast ; ii) produce the corresponding skymaps containing only thestructure relative to multipoles up to (cid:96) = 5 ( (cid:96) = 10):this is obtained as a constrained realization by feedingthe HEALPix routine alm2map with the a lm ’s obtainedin point i) only up to (cid:96) = 5 ( (cid:96) = 10); iii) subtract themaps derived in point ii) from the original maps. Thisis an approximate way to subtract lower multipoles, be-cause the effect of the mask is not included in i) but itis useful to test the impact of possible leakages from lowto high multipoles. From these maps we then derive thecross-correlation APS and compare it with the APS de-termined with only the monopole or monopole and dipolesubtracted (the latter is what we do in our baseline anal-ysis). The results are shown in the bottom left panel ofFig. 14, which shows that all results are perfectly com-patible with each other and therefore there is no leakageof power from lower multipoles to the multipole windowof interest.
4. Correlation in real space
In order to test the robustness of our measurement, wealso compute the cross-correlation function in real space ξ ( θ ), which can then be transformed to the APS with theusual relation: ξ ( θ ) = (cid:88) l (2 l + 1)4 π C l P l (cos θ ) , (A2)where P l (cos θ ) are the Legendre polynomials and θ isthe physical angular scale. The correlation function ξ ( θ ) is determined by means of the following estimator: ξ ( θ ) = 1 (cid:80) a,b f ab ( θ ) (cid:88) a,b ( n γ − ¯ n γ ) ( n gal − ¯ n gal )¯ n gal f ab ( θ ) , (A3)where ( n γ − ¯ n γ ) and ( n gal − ¯ n gal ) represent the fluctua-tions of the γ -ray intensity flux and of the galaxy num-ber counts in every unmasked a -th and b -th pixel andthe function f ab ( θ ) assumes the value 1 when the angu-lar separation of the two pixels is θ and 0 otherwise. Wecompare the correlation function we measure by means ofEq. A3 with the corresponding ξ ( θ ) provided by PolSpice.The error associated to our estimator is computed with ajackknife re-sampling approach, dividing the sky into 20distinct patches and estimating the relative covariance.The bottom right panel in Fig. 14 shows the comparisonof the two methods: they nicely agree, with the jackknifemethod possibly underestimating the errors.
5. Comparison with previous measurement
Finally we compare our measurement with the resultsobtained in a previous analysis of cross-correlation be-tween γ -rays and the 2MPZ catalog, with a smaller pho-ton statistics [40]. The comparison is shown in Fig. 15and refers to the determination of the Poisson noise terms C kp defined in section III A and is performed for the full2MPZ sample. We see that our results and the results ofRef. [40] are in good agreement, and we can appreciatethe improvement in the statistical determination of thesignal with our new analysis. Appendix B: Halo occupation distribution ofgalaxies
In this work, we adopt the halo model to describe theclustering of structures. In order to estimate the angularcross-correlation of the unresolved γ -ray sky with sam-ples of galaxies, we need to describe how galaxies popu-late halos. To this end, we employ the halo occupationdistribution (HOD) formalism.We follow the approach described in Ref. [65] (for re-view on HOD, see also Refs. [66, 67]), where the HODis parameterized by distinguishing the contributions ofcentral and satellite galaxies: N = N cen + N sat , since dif-ferent formation histories typically imply different prop-erties for galaxies residing at the centers of halos withrespect to satellite galaxies. These can be modeled withthe following functional form for the central galaxies: (cid:104) N cen ( M ) (cid:105) = 12 (cid:20) (cid:18) log M − log M cut σ logM (cid:19)(cid:21) (B1)9 FIG. 14. Tests of stability of our results. Upper left panel: APS and its errors determined with the theoretical Gaussian estimate(red) and the error provided by PolSpice (blue). Upper central panel: APS obtained by using the masked photon maps with(blue) and without (pink) Galactic foreground subtraction. Upper right panel: Same as Fig. 5 bottom-right (blue-points), butemploying different models for the Galactic foreground [4]. Lower left panel: APS obtained by removing the monopole (darkblue), the monopole and the dipole (light blue), the first 5 multipoles (red), and the first 10 multipoles (pink). Lower rightpanel: Angular correlation function measured with the estimator of Eq. A3 (yellow) and with PolSpice (blue); the errors forthe former are determined by means of a jackknife technique. All results refer to a [1 ,
10] GeV energy bin (except for upperright panel) and to the full 2MPZ catalog. and with the following form for the satellite galaxies: (cid:104) N sat ( M ) (cid:105) = (cid:18) M − M cut M (cid:19) α for M > M cut (cid:104) N sat ( M ) (cid:105) = 0 for M ≤ M cut . (B2)With this formalism, we need four parameters for eachgalaxy population: M cut denotes the approximate halomass required to populate the halo with the consideredtype of galaxy, with the transition from 0 to 1 centralgalaxy modeled by means of Eq. (B1), and set by thewidth σ LogM . The satellite occupation is described by apower law (with index α and normalization set by themass parameter M .Eqs. (B1) and (B2) provide the number of galaxies ina halo of mass M . Concerning the spatial distribution,we treat central and satellite galaxies separately. Theformer is taken as a point-source located at the centerof the halo (the point-source approximation is expectedto break down only for (cid:96) (cid:38) ). Satellite galaxies areinstead described in an effective way with a spatial dis- tribution following the host-halo profile. In other words,we express the density field of galaxies with: g g ( x − x (cid:48) | M ) = (cid:104) N cen ( M ) (cid:105) δ ( x − x (cid:48) ) + (cid:104) N sat ( M ) (cid:105) ρ h ( x − x (cid:48) | M ) /M . (B3)Note that: (cid:90) d x g g ( x ) = (cid:104) N cen ( M ) (cid:105) + (cid:104) N sat ( M ) (cid:105) = (cid:104) N ( M ) (cid:105) . (B4)The value of the four HOD parameters of each sam-ple is derived by fitting the auto-correlation of the spe-cific catalog. We perform the measurement of the auto-correlation by employing the PolSpice tool, in the sameway as described in the main text for the cross correla-tion. The noise term is estimated with C ggN = 4 π f sky /N ,where N is the total number of galaxies outside the mask,and is subtracted from the measurement.The theoretical prediction for the 3D power spectrumis computed with the halo model approach (and assumingPoisson statistics) as:0 FIG. 15. Comparison with previous results for the Poissonnoise terms C kp as a function of the energy for the full 2MPZsample. Points refer to our analysis, the shaded regions showthe results obtained in the previous analysis of Ref. [40]. P hgg ( k, z ) = (cid:90) M max M min dM dndM (cid:104) N cen (cid:105) (cid:104) N sat (cid:105) ˜ v δ ( k | M ) + (cid:104) N sat (cid:105) ˜ v δ ( k | M ) ¯ n g (B5) P hgg ( k, z ) = (cid:34)(cid:90) M max M min dM dndM b h ( M ) (cid:104) N g (cid:105) ¯ n g ˜ v g ( k | M ) (cid:35) P lin ( k ) . (B6)The product (cid:104) N g (cid:105) ˜ v g ( k | M ) is the Fourier transform of (cid:104) N cen ( M ) (cid:105) δ ( x ) + (cid:104) N sat ( M ) (cid:105) ρ h ( x | M ) /M . Note that (cid:104) N g (cid:105) ˜ v g ( k = 0 | M ) = (cid:104) N g (cid:105) . The average numberof galaxies at a given redshift is given by ¯ n g ( z ) = (cid:82) dM dn/dM (cid:104) N g (cid:105) . Note that in Eq. B1, we do not in-clude the shot-noise term ∝ (cid:104) N (cid:105) since it has been sub-tracted from the data.The best-fit HOD parameters are reported in Ta-ble VII. A few examples of the comparison between the-oretical model and measured APS are shown in Fig. 16.It is clear from the plot that the models are stronglyconstrained by the measurements. Therefore, the uncer-tainty on the HOD parameters has negligible impact onthe cross-correlation observable and can be neglected inour analysis, where we consider only the best-fit values. Appendix C: Estimate of gamma-ray luminosityfrom other wavelengths
As mentioned in the main text, the Poisson noise termof the cross-correlation signal is given by the averagegamma-ray flux of objects in the catalogs. The computa-tion is performed in two steps. First, we derive a relationfor the (diffuse) gamma-ray production of all galaxiesgiven some tracer of the star formation rate. Then we
TABLE VII. Best fit values of the HOD parameters ofEqs.(B1) and (B2) for all the samples considered in this work.Catalog M cut σ Log M α M [10 M (cid:12) ] [10 M (cid:12) ]2MPZ (full) 1.8 0.32 1.15 2.82MRS (full) 1.6 0.22 1.0 2.02MPZ high- z z z add up emissions from blazars and misaligned AGN ifthe object has been classified as an host of these emit-ters.Here we describe how we derive the gamma-ray emis-1
50 100 500 1000
Multipole l -2 -1 l C l gg FIG. 16. Autocorrelation angular power spectrum for the2MPZ (full, mid-B, mid-K, mid- z ) and 2MRS catalogs. Pointsshow the measurements, while lines refer to the best-fit modelderived as described in the text. sion of AGNs and star-forming galaxies starting from agiven magnitude in the optical/infrared. Note that suchrelations suffer from significant uncertainty. If the latteris due just to random scatter around the reported rela-tions, the impact of these uncertainty in our analysis issubdominant. In fact, in order to compute the Poissonnoise term, we add up the flux of a very large number ofobjects. On the other hand, if the adopted relations arebiased, this could in principle affect our conclusion. Toovercome this issue, we introduce also a model in whichthe Poisson noise term is not modeled but left free andfitted.
1. Blazars
The gamma-ray flux of blazars is computed using therelation between the infrared magnitude at 12 µ m and the energy flux between 0.1 GeV and 100 GeV foundin Ref. [63]. From their Fig. 2, one obtains a F Eγ = A [W3] − β with A = 10 − . ± . erg cm − s − and β =4 . ± .
17. We employ the W3 magnitude measured bythe WISE survey and provided in the catalog.
2. Misaligned AGNs
Predictions for the gamma-ray flux of misalignedAGNs are typically derived from their radio emission[16, 68], with the best-fit relation found to be: L γ =10 − . ( L RC / erg s) . , where L γ is the luminosity be-tween 0.1 GeV and 100 GeV and L RC is the 5GHz radiocore luminosity. Ref. [69] shows a correlation between 1.4GHz luminosity and the 12 µ m luminosity of WISE AGNs(see their Fig. 13). These two relations allow us to pre-dict the Poisson noise term of mAGNs starting from theW3 magnitude of the 2MPZ catalog. The predicted av-erage gamma-ray flux agrees well with a more direct esti-mate we obtained on a smaller sample obtained by cross-matching the 2MPZ sources with the FIRST catalog [70],to directly extract radio fluxes (then linked to gamma-rayfluxes using again the relation of Refs. [16, 68]).
3. Star-forming galaxies
Star formation is expected to trigger gamma-ray pro-duction in galaxies. Indeed, galaxies detected in γ -raysshow a tight correlation between the luminosity in therange (0 . − L γ = (1 . ± . × (SFR / M (cid:12) yr) . ± . erg/s [11].In turn, the star formation rate is correlated with theB-band magnitude (see, e.g., Fig. 5 in Ref. [52]). InRef. [52], they found L B = 13 . × SFR / M (cid:12) yr with ascatter within one dex. We estimate the average gamma-ray flux of star forming galaxies starting from the B-bandmagnitude reported in the 2MPZ catalog and using theabove two relations. [1] M. Fornasa and M. A. S´anchez-Conde, Phys. Rept. ,1 (2015), arXiv:1502.02866 [astro-ph.CO].[2] F. Acero et al. (Fermi-LAT), Astrophys. J. Suppl. ,23 (2015), arXiv:1501.02003 [astro-ph.HE].[3] M. Ackermann et al. (Fermi-LAT), Phys. Rev. Lett. ,151105 (2016), arXiv:1511.00693 [astro-ph.CO].[4] M. Ackermann et al. (Fermi-LAT), Astrophys. J. , 86(2015), arXiv:1410.3696 [astro-ph.HE]. [5] M. Ackermann et al. , Astrophys. J. , 30 (2011),arXiv:1108.0501 [astro-ph.CO].[6] M. Ajello et al. , Astrophys. J. , 108 (2012),arXiv:1110.3787 [astro-ph.CO].[7] M. Ajello et al. , Astrophys. J. , 73 (2014),arXiv:1310.0006 [astro-ph.CO].[8] M. Di Mauro, F. Calore, F. Donato, M. Ajello, and L. La-tronico, Astrophys. J. , 161 (2014), arXiv:1304.0908[astro-ph.HE]. [9] M. Di Mauro, F. Donato, G. Lamanna, D. A. Sanchez,and P. D. Serpico, Astrophys. J. , 129 (2014),arXiv:1311.5708 [astro-ph.HE].[10] S. Ando and E. Komatsu, Phys. Rev. D73 , 023521(2006), arXiv:astro-ph/0512217 [astro-ph].[11] M. Ackermann et al. (Fermi-LAT), Phys. Rev.
D85 ,083007 (2012), arXiv:1202.2856 [astro-ph.HE].[12] J. P. Harding and K. N. Abazajian, JCAP , 026(2012), arXiv:1206.4734 [astro-ph.HE].[13] A. Cuoco, E. Komatsu, and J. M. Siegal-Gaskins,Phys. Rev.
D86 , 063004 (2012), arXiv:1202.5309 [astro-ph.CO].[14] 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].[15] S. Ando and E. Komatsu, Phys. Rev.
D87 , 123539(2013), arXiv:1301.5901 [astro-ph.CO].[16] M. Di Mauro, A. Cuoco, F. Donato, and J. M.Siegal-Gaskins, JCAP , 021 (2014), arXiv:1407.3275[astro-ph.HE].[17] M. Fornasa et al. , Phys. Rev.
D94 , 123005 (2016),arXiv:1608.07289 [astro-ph.HE].[18] S. Ando, M. Fornasa, N. Fornengo, M. Regis,and H.-S. Zechlin, Phys. Rev.
D95 , 123006 (2017),arXiv:1701.06988 [astro-ph.HE].[19] S. Dodelson, A. V. Belikov, D. Hooper, and P. Serpico,Phys. Rev.
D80 , 083504 (2009), arXiv:0903.2829 [astro-ph.CO].[20] D. Malyshev and D. W. Hogg, Astrophys. J. , 181(2011), arXiv:1104.0010 [astro-ph.CO].[21] M. Lisanti, S. Mishra-Sharma, L. Necib, and B. R. Safdi,Astrophys. J. , 117 (2016), arXiv:1606.04101 [astro-ph.HE].[22] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo,and M. Regis, Astrophys. J. , L31 (2016),arXiv:1605.04256 [astro-ph.HE].[23] H.-S. Zechlin, A. Cuoco, F. Donato, N. Fornengo,and A. Vittino, Astrophys. J. Suppl. , 18 (2016),arXiv:1512.07190 [astro-ph.HE].[24] M. Di Mauro, S. Manconi, H. S. Zechlin, M. Ajello,E. Charles, and F. Donato, Astrophys. J. , 106(2018), arXiv:1711.03111 [astro-ph.HE].[25] J.-Q. Xia, A. Cuoco, E. Branchini, M. Fornasa, andM. Viel, Mon. Not. Roy. Astron. Soc. , 2247 (2011),arXiv:1103.4861 [astro-ph.CO].[26] S. Camera, M. Fornasa, N. Fornengo, and M. Regis,Astrophys. J. , L5 (2013), arXiv:1212.5018 [astro-ph.CO].[27] S. Ando, A. Benoit-L´evy, and E. Komatsu, Phys. Rev.
D90 , 023514 (2014), arXiv:1312.4403 [astro-ph.CO].[28] S. Ando, JCAP , 061 (2014), arXiv:1407.8502[astro-ph.CO].[29] N. Fornengo, L. Perotto, M. Regis, and S. Camera,Astrophys. J. , L1 (2015), arXiv:1410.4997 [astro-ph.CO].[30] M. Shirasaki, S. Horiuchi, and N. Yoshida, Phys. Rev.
D90 , 063502 (2014), arXiv:1404.5503 [astro-ph.CO].[31] S. Camera, M. Fornasa, N. Fornengo, and M. Regis,JCAP , 029 (2015), arXiv:1411.4651 [astro-ph.CO].[32] A. Cuoco, J.-Q. Xia, M. Regis, E. Branchini, N. For-nengo, and M. Viel, Astrophys. J. Suppl. , 29 (2015),arXiv:1506.01030 [astro-ph.HE]. [33] M. Regis, J.-Q. Xia, A. Cuoco, E. Branchini, N. For-nengo, and M. Viel, Phys. Rev. Lett. , 241301 (2015),arXiv:1503.05922 [astro-ph.CO].[34] J.-Q. Xia, A. Cuoco, E. Branchini, and M. Viel, Astro-phys. J. Suppl. , 15 (2015), arXiv:1503.05918 [astro-ph.CO].[35] M. Shirasaki, S. Horiuchi, and N. Yoshida, Phys. Rev.
D92 , 123540 (2015), arXiv:1511.07092 [astro-ph.CO].[36] C. Feng, A. Cooray, and B. Keating, Astrophys. J. ,127 (2017), arXiv:1608.04351 [astro-ph.CO].[37] T. Troster et al. , (2016), 10.1093/mnras/stx365,[Mon. Not. Roy. Astron. Soc.467,no.3,2706(2017)],arXiv:1611.03554 [astro-ph.CO].[38] M. Shirasaki, O. Macias, S. Horiuchi, S. Shirai,and N. Yoshida, Phys. Rev.
D94 , 063522 (2016),arXiv:1607.02187 [astro-ph.CO].[39] E. Branchini, S. Camera, A. Cuoco, N. Fornengo,M. Regis, M. Viel, and J.-Q. Xia, Astrophys. J. Suppl. , 8 (2017), arXiv:1612.05788 [astro-ph.CO].[40] A. Cuoco, M. Bilicki, J.-Q. Xia, and E. Bran-chini, (2017), 10.3847/1538-4365/aa8553, [Astrophys. J.Suppl.232,10(2017)], arXiv:1709.01940 [astro-ph.HE].[41] M. Shirasaki, O. Macias, S. Horiuchi, N. Yoshida, C.-H. Lee, and A. J. Nishizawa, Phys. Rev.
D97 , 123015(2018), arXiv:1802.10257 [astro-ph.CO].[42] D. Hashimoto, A. J. Nishizawa, M. Shirasaki, O. Ma-cias, S. Horiuchi, H. Tashiro, and M. Oguri, (2018),arXiv:1805.08139 [astro-ph.CO].[43] P. A. R. Ade et al. (Planck), Astron. Astrophys. et al. (Fermi-LAT), Astrophys. J. Suppl. ,18 (2017), arXiv:1702.00664 [astro-ph.HE].[47] See: https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/for further details.[48] https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html.[49] M. Bilicki, T. H. Jarrett, J. A. Peacock, M. E. Clu-ver, and L. Steward, Astrophys. J. Suppl. , 9 (2014),arXiv:1311.5246 [astro-ph.CO].[50] D. Alonso, A. I. Salvador, F. J. S´anchez, M. Bilicki,J. Garc´ıa-Bellido, and E. S´anchez, Mon. Not. Roy. As-tron. Soc. , 670 (2015), arXiv:1412.5151.[51] J. P. Huchra, L. M. Macri, K. L. Masters, T. H. Jar-rett, P. Berlind, M. Calkins, A. C. Crook, R. Cutri,P. Erdoˇgdu, E. Falco, T. George, C. M. Hutcheson,O. Lahav, J. Mader, J. D. Mink, N. Martimbeau,S. Schneider, M. Skrutskie, S. Tokarz, and M. Westover,Astrophys. J. Suppl. , 26 (2012), arXiv:1108.0669.[52] M. S. Bothwell, R. C. Kennicutt, and J. C. Lee, Mon.Not. Roy. Astron. Soc. , 154 (2009), arXiv:0908.1122[astro-ph.CO].[53] M. Ackermann et al. (Fermi-LAT), Astrophys. J. ,164 (2012), arXiv:1206.1346 [astro-ph.HE].[54] P. S. Behroozi, C. Conroy, and R. H. Wechsler, Astro-phys. J. , 379 (2010), arXiv:1001.0015 [astro-ph.CO].[55] R. D’Abrusco, F. Massaro, A. Paggi, H. A. Smith,N. Masetti, M. Landoni, and G. Tosti, Astrophys. J.Suppl. , 14 (2014), arXiv:1410.0029 [astro-ph.HE].[56] R. Edelson and M. Malkan, Astrophys. J. , 52 (2012),arXiv:1203.1942. , 29 (2015),arXiv:1506.01030 [astro-ph.HE].[61] M. Ajello, D. Gasparrini, M. Sanchez-Conde, G. Zahari-jas, M. Gustafsson, et al. , Astrophys.J. , L27 (2015),arXiv:1501.05301 [astro-ph.HE].[62] J. Zuntz, M. Paterno, E. Jennings, D. Rudd, A. Man-zotti, S. Dodelson, S. Bridle, S. Sehrish, andJ. Kowalkowski, Astron. Comput. , 45 (2015),arXiv:1409.3409 [astro-ph.CO].[63] F. Massaro and R. D’Abrusco, Astrophys. J. , 67(2016), arXiv:1609.08615 [astro-ph.HE].[64] N. Fornengo and M. Regis, Front. Physics , 6 (2014),arXiv:1312.4835 [astro-ph.CO]. [65] Z. Zheng, A. A. Berlind, D. H. Weinberg, A. J. Ben-son, C. M. Baugh, et al. , Astrophys.J. , 791 (2005),arXiv:astro-ph/0408564 [astro-ph].[66] A. A. Berlind and D. H. Weinberg, Astrophys.J. , 587(2002), arXiv:astro-ph/0109001 [astro-ph].[67] A. Cooray and R. K. Sheth, Phys.Rept. , 1 (2002),arXiv:astro-ph/0206508 [astro-ph].[68] D. Hooper, T. Linden, and A. Lopez, JCAP , 019(2016), arXiv:1604.08505 [astro-ph.HE].[69] B. Mingo, M. G. Watson, S. R. Rosen, M. J. Hardcastle,A. Ruiz, A. Blain, F. J. Carrera, S. Mateos, F.-X. Pineau,and G. C. Stewart, Mon. Not. Roy. Astron. Soc. ,2631 (2016), arXiv:1607.06471.[70] R. H. Becker, R. L. White, and D. J. Helfand, Astrophys.J.450