Modelling the flux distribution function of the extragalactic gamma-ray background from dark matter annihilation
PPrepared for submission to JCAP
Modelling the flux distributionfunction of the extragalacticgamma-ray background from darkmatter annihilation
Michael R. Feyereisen, a Shin’ichiro Ando a and Samuel K. Lee b,c a GRAPPA Institute, University of Amsterdam, Science Park 904, 1098 XH Amsterdam,Netherlands b Princeton Center for Theoretical Science, Princeton University, Princeton, NJ 08544, USA c Broad Institute, 75 Ames Street, Cambridge, MA 02142, USAE-mail: [email protected], [email protected], [email protected]
Abstract.
The one-point function (i.e., the isotropic flux distribution) is a complementarymethod to (anisotropic) two-point correlations in searches for a gamma-ray dark matterannihilation signature. Using analytical models of structure formation and dark matter haloproperties, we compute the gamma-ray flux distribution due to annihilations in extragalacticdark matter halos, as it would be observed by the Fermi Large Area Telescope. Combiningthe central limit theorem and Monte Carlo sampling, we show that the flux distribution takesthe form of a narrow Gaussian of ‘diffuse’ light, with an ‘unresolved point source’ power-lawtail as a result of bright halos. We argue that this background due to dark matter constitutesan irreducible and significant background component for point-source annihilation searcheswith galaxy clusters and dwarf spheroidal galaxies, modifying the predicted signal-to-noiseratio. A study of astrophysical backgrounds to this signal reveals that the shape of the totalgamma-ray flux distribution is very sensitive to the contribution of a dark matter component,allowing us to forecast promising one-point upper limits on the annihilation cross section. Weshow that by using the flux distribution at only one energy bin, one can probe the canonicalcross section required for explaining the relic density, for dark matter of masses around tensof GeV.
Keywords:
Dark Matter, Fermi-LAT, Structure Formation, Substructure Boost, Blazar a r X i v : . [ a s t r o - ph . C O ] A ug ontents P ( F ) 94.1.1 Dominant effects 94.1.2 Subleading effects 94.2 Computing P ( F ) 114.2.1 Monte Carlo method combined with the central limit theorem 114.2.2 Flux distribution and instrumental sensitivity 12 P ( C ) 205.4 Caveats 22 B.1 Rationale of our method 25B.2 Deriving the form of P k ( F ) 26B.3 Cross checks 27 C Angular resolution limit 27 – 1 –
Introduction
The Large Area Telescope (LAT) onboard the Fermi satellite [1] measured the energy spec-trum [2] and angular anisotropies [3] of the diffuse extragalactic background of gamma rays.Components contributing to this background include blazars [4, 5], star-forming and star-burst galaxies [6], and misaligned active galaxies [7]. The combination of these sources givesreasonably good fit to the spectral data [5, 8], while the anisotropies are consistent with theblazar component alone [3, 9]. Independently of this, dark matter has emerged as the pre-ferred explanation of many astrophysical and cosmological features through gravity (galacticrotation curves, Ω m (cid:29) Ω b , lensing by galaxy clusters, etc.). If particle dark matter producesgamma rays (e.g., by self-annihilation) as in the case of weakly interacting massive particles(WIMPs) motivated by popular particle-physics models [10, 11], then it could also contributeto this diffuse signal (in some unknown proportion) [12]. Given that the known astrophysicalsources yield reasonable fit to the spectrum of the gamma-ray background, the dark mattercomponent started to be tightly constrained only through the spectral data (e.g., [5, 8, 13]).Recently, new analysis techniques beyond the energy spectrum and angular two-pointcorrelations were proposed and investigated extensively. Among them is to take cross correla-tions of gamma-ray data with local galaxy catalogs [14, 15] and matter distribution throughlensing data [16, 17]. Although recent measurements of the cross correlations [18–21] areconsistent with the hypothesis of no dark matter signal, they yield tight constraints thereof(e.g., [22]).Complementary to the studies on these two-point functions, the one-point function(i.e. the photon-count or flux distribution) would leverage the isotropic component of thediffuse signal. For example, the flux distribution of Milky-Way subhalos has been used toconstrain particle dark matter properties in light of Fermi unidentified sources [23]. Theone-point function of Fermi-LAT data has been experimentally fit to a combination of thediffuse background and blazar-like sources by Ref. [24]. Theoretically, Ref. [25] studiedthe one-point probability density function (PDF) of the gamma-ray flux due to Galacticsubhalos, and showed that it features power-law tail at high-flux end (see also Refs. [26, 27]).Understanding the one-point PDF for all the relevant sources will be important also forpossible detections of dwarfs or galaxy clusters with gamma rays (e.g., [28]).In this paper, we extend the theoretical framework of Ref. [25] to include the contributionof dark matter annihilation in the extragalactic halos. We model the gamma-ray flux fromthe population of dark matter halos using the mass and luminosity functions predicted for thestructure formation scenario in the Universe with cold dark matter and cosmological constant(ΛCDM). Since there are many such halos, the total flux observed by Fermi is predicted using‘large N ’ statistical tools. By combining the central limit theorem (CLT) in the low-fluxregime and a Monte Carlo method in the high-flux regime, we find that the differential fluxdistribution has a roughly Gaussian peak, as result of the diffuse emission of a large numberof very faint sources, but with a power-law high-flux tail due to the rare occurrence of anexceptionally bright halo. The all-sky flux in our fiducial model lies around the sensitivitylimit of the Fermi, well-below the measured gamma-ray background [2], and at roughly thesame level as the dark matter flux expected from the Fornax galaxy cluster. We find thatthe detectability of a dark matter signature from galaxy clusters over the extragalactic darkmatter background decreases, when the luminosity boost due to halo substructure increases.We also illustrate how to disentangle a dark matter signal from astrophysical backgrounds inthe presence of the photon shot noise, and forecast the upper limit on the annihilation cross– 2 –ection one might expect to obtain using the one-point function alone (although actuallyperforming this analysis with its due rigour is beyond the scope of this theoretical paper).Given a fiducial model for the dark matter halo substructure boost, the particle dark mattermass, etc., we find a 5 σ upper bound roughly a factor of two above the thermal cross section.This paper is organised as follows. In Sec. 2, we construct the model of the flux dis-tribution observed at the Fermi-LAT. In Sec. 3, we detail our specific model choices, andin Sec. 4, we present our main results, including a sensitivity analysis of our distribution tothe model choices, and a probabilistic method for summing the fluxes from a (quite literally)astronomically large number O (10 ) of halos. In Sec. 5, we discuss consequences of thisstudy for indirect DM searches, by comparing the predicted distribution to the gamma-rayfluxes of galaxy clusters, dwarf spheroidals, and blazars. We conclude the paper in Sec. 6. The goal of the present work is to theoretically predict the PDF P ( F ), which gives theprobability of observing a total gamma-ray flux F arising from dark matter annihilation inextragalactic halos in a Fermi pixel of a particular size. In this section, we present a formalismfor constructing P ( F ) given models for the cosmology, halo properties, and annihilationprocess; we proceed in a completely general manner, postponing specific choices for thesemodels until Sec. 3.We construct a Bayesian hierarchical model, to predict unknown gamma-ray observablesfrom well-constrained ΛCDM parameters and fitted models of N-body simulations. In thehierarchical Bayesian approach, uncertainties in the parameters of probability distributionsare modelled with their own distributions (and recursively). This allows us to systematicallycombine the uncertainties on physics at widely differing scales, and thereby to perform asensitivity analysis of our model (Sec. 4.1). Throughout the paper, F represents the differential flux, i.e., a number of photons receivedper unit area, unit time, and unit energy range [ F ( E ) = d N γ /dAdtdE ]. The PDF P ( F )for observing a total differential flux F from all of the halos in a pixel depends on the PDF P ( F ) for observing F from any individual halo. We thus proceed by first deriving the latterquantity.Because the differential flux F from an individual halo is completely determined by itsrest-frame differential luminosity L = d N γ /dtdE and its redshift z , we can write P ( F ) = (cid:90) dL dz P ( F | L, z ) P ( L, z )= (cid:90) dL dz δ [ F − F ( L, z )] P ( L | z ) P ( z ) . (2.1)Here, the usual relation for the differential flux, F ( E ; L, z ) = e − τ ( E,z ) (1 + z ) L [(1 + z ) E ]4 πd L ( z ) , (2.2) Throughout this paper, we denote probability distributions by P ( · · · ) and distinguish them using therandom variables that they describe, along with subscripts if necessary. Conditional and parameterised dis-tributions are denoted as P ( ·|· ). Exceptions to this convention are Poisson and Normal distributions, denoted P ( ·|· ) and G ( ·|· , · ) respectively. – 3 –epends on the luminosity distance d L ( z ) and the pair-production optical depth τ ( E, z )for gamma-ray photons, and also accounts for the redshift of photons emitted with rest-frame energies E (1 + z ) to observed energies E . We can interpret P ( L | z ) = dN/dL ( z ) asthe redshift-dependent halo differential-luminosity function. Assuming that the halos areisotropically distributed across the Universe, the number of halos at redshift z is proportionalto the comoving volume δV ( z ) of the corresponding redshift slice δz , therefore we also have P ( z ) = dN/dz ∝ dV /dz .Alternatively, we can rewrite Eq. (2.1) in terms of the halo mass M to obtain P ( F ) = (cid:90) dL dM dz δ [ F − F ( L, z )] P ( L | M, z ) P ( M | z ) P ( z ) , (2.3)where we can similarly interpret P ( M | z ) = dN/dM ( z ) as the redshift-dependent halo massfunction.In principle, the distribution P ( L | M, z ) in Eq. (2.3) captures the scatter in the relationbetween the differential luminosity and the mass of a halo, which also depends on redshift.This is because the halo luminosity is determined not only by the properties of the darkmatter particle and the details of the annihilation process, but also by the density profile ρ ofthe halo, which usually shows scatter for any given M . The halo profiles can be completelycharacterised by some parameters θ h (such as ρ s , r s , r vir . . . in the case of the NFW profile [29])so that (for any given particle dark matter model) we have L = L ( θ h ). If we further assumethat the distribution of halo profiles can be described by a halo model that gives P ( θ h | M, z ),we can write P ( L | M, z ) = (cid:90) d θ h P ( L | θ h ) P ( θ h | M, z )= (cid:90) d θ h δ [ L − L ( θ h )] P ( θ h | M, z ) . (2.4)We can then use this expression to simplify Eq. (2.1), giving P ( F ) = (cid:90) dM dz d θ h δ [ F − F ( θ h , z )] P ( θ h | M, z ) P ( M | z ) P ( z ) , (2.5)where the flux relation is now written in terms of θ h .In order to make the numerical calculation of Eq. (2.5) more tractable, we shall neglectthe scatter in the distribution P ( θ h | M, z ) in this work. That is, we take the distribution ofthe halo-profile parameters θ h = { θ h, , . . . , θ h,n } to be given by P ( θ h | M, z ) = n (cid:89) i =1 δ [ θ h,i − ¯ θ h,i ( M, z )] , (2.6)where the functions ¯ θ h,i ( M, z ) give the mean values for the parameters. With this assumption,we can perform the integrals over θ h and M in Eq. (2.5), leaving only an integral over z : P ( F ) = (cid:90) dz (cid:12)(cid:12)(cid:12)(cid:12) ∂F∂M (cid:12)(cid:12)(cid:12)(cid:12) − dNdM dVdz . (2.7) Here and elsewhere, equalities between probability densities and number densities is meant modulo anormalisation. – 4 –ere, functions of M in the integrand are evaluated at the value of M defined implicitly bythe flux relation F = F ( ¯ θ h ( M, z ) , z ) for the corresponding values of F and z .Eq. (2.7) then gives the PDF P ( F ) for the differential flux from individual halos. Toreiterate, this is given in terms of (i) the cosmological model, which affects the differentialflux F ( θ h , z ), the volume P ( z ), and the mass function P ( M | z ); (ii) the halo model, whichgives the mean values ¯ θ h ( M, z ) of the halo-profile parameters; and (iii) the optical depth andthe details of the annihilation process, which affect the normalisation of F ( θ h , z ). In Sec. 3,informed by observations and simulations, we shall make specific, fiducial choices for thesequantities and calculate the resulting P ( F ). However, before doing so, we shall completeour discussion of our general formalism by reviewing how P ( F ) can be used to find P ( F ),the PDF of the total flux F from all of the halos in a pixel. We assume that all dark matter gamma-ray sources may be treated as point sources (see alsoSec. 3.1.1 and Sec. 5.4 below), allowing us to equate the differential flux F and the differentialflux per pixel. The dark matter differential flux F arriving at any given pixel of the Fermisky map, is the summed flux F = (cid:80) i F i of any number of individual halo point sources [30],where each differential flux F i is an independent and identically distributed (i.i.d.) randomvariable with the distribution P ( F ). The distribution of a sum of random variables is theconvolution of all the original distributions [31]; Since the F i are i.i.d., the distribution of thetotal differential flux per pixel is the autoconvolution [25] P k ( F ) = P ( F ) (cid:63) P ( F ) (cid:63) · · · (cid:63) P ( F ) = ( P ) (cid:63)k , (2.8)where k is the number of halos contributing to this flux. Since furthermore we do not knowhow many halos are thus stacked in a pixel, the number k of fluxes in the sum is itself arandom variable. If we assume this number k of halos per pixel is Poisson-distributed overthe sky with some mean N (cid:48) , we can model the total differential flux per pixel as P ( F ) = (cid:90) dN (cid:48) P ( N (cid:48) ) (cid:88) k P (cid:0) k | N (cid:48) (cid:1) P k ( F ) , (2.9)where the uncertainties in k and N (cid:48) are marginalised away. Since the numbers k and N (cid:48) ofextragalactic halos are very large, both P ( N (cid:48) ) and P ( k | N (cid:48) ) are thin enough to be approx-imated by delta functions, so that P ( F ) = P k ≈ N (cid:48) ( F ). Thus, the only additional physicalinput required to compute P ( F ) from P ( F ) is N (cid:48) , which is discussed below. Cosmology directly determines the redshift distribution of halos (via isotropy P ( z ) = dN/dz )and their mass distribution (via the gravitational collapse of inhomogeneities that yields dN/dM ). The normalisation of these number densities clearly corresponds to the total num-ber N of halos in the Universe. The number of halos per Fermi pixel is then N (cid:48) = Ω pix π N = Ω pix π (cid:90) dNdM dVdz dM dz , (3.1)– 5 –here Ω pix is the pixel size (expressed in units of a solid angle) and the mass function dN/dM is described in Sec. 3.1.2. The parameters of the ΛCDM cosmology relevant to our modelare (Planck+WMAP from [32]) { Ω b h , Ω c h , Ω Λ , H , σ , n S } = { . , . , . , . , . , . } . (3.2)Furthermore, the integration limits for Eq. (3.1) (and also Eq. (2.7)) are chosen as follows:we assume that our dark matter candidate forms structures down to 10 − M (cid:12) [33–37], andallow for virialised dark matter structures up to 10 M (cid:12) , 100 times more massive than galaxyclusters. We assume that we can measure luminosity from structures that form between z = 5and 10 − , the latter of which corresponds to a distance of roughly 45 kpc, well outside of thebaryonic content of the Milky Way in any direction.We model the energy-dependent angular resolution of Fermi-LAT [1] as follows: θ ( E ) = (cid:40) . ◦ ( E/ GeV) − . for 0 . < E/ GeV < . ◦ for 20 < E/ GeV , (3.3)such that θ ≈ . ◦ at our fiducial observing energy (justified in Sec. 4.1) of E = 1 GeV.This is slightly larger than 0 . ◦ quoted in [1], which is valid for normally incident photonsonly. Since we adopt this angular resolution as a size of each pixel, Ω pix ≈ πθ ( E ), we entirelyneglect the instrumental point-spread function. This is also justified because the number ofsources per pixel is found very large ( N (cid:48) ∼ × ) and the flux is diffuse.The largest possible mass of a point source at a given redshift M Ext ( z ) may be deter-mined from the critical virial radius R Ext ≈ d L ( z ) tan( θ/
2) that fits in a pixel. We can usethis as an integration limit in Eq. (3.1) to cross-check our assumption (Sec. 2.2) that all darkmatter gamma-ray sources are point-like. We find less than 0.28 extended sources per pixel.This is much more than twenty orders of magnitude smaller than the total number of sources,but still represents a non-negligible absolute amount given the large number of pixels. Wediscuss this issue more thoroughly in 5.4.
The halo mass function, first addressed heuristically by Press and Schechter [38] and subse-quently formalised in, e.g., Ref. [39], is computed as dndM = ¯ ρM f ( ν ) dνdM , ν = (cid:18) δ c σ (cid:19) , (3.4)where δ c = 1 . /D ( z ) is the (linear) critical overdensity, D ( z ) is the linear growth factor, and σ ( M ) is the rms deviation of primordial density fluctuations, smoothed to scale M [39]. Thefunctional form of σ ( M ) (required to calculate dν/dM ) is determined from the literature[40] with normalisation set by the cosmological parameter σ [32]. The function f ( ν ) isderived from the excursions of these density fluctuations above a ‘barrier’ [39] that encodesthe physics of halo collapse (including δ c ). For an approachable presentation of the formalism,see Ref. [41].In addition to ellipsoidal collapse, our fiducial mass function incorporates a virialisedhalo’s angular momentum and the cosmological constant into its barrier δ c . It has a self-similar f ( ν ) well-fit by the following function (Eq. (163) in [42]): νf ( ν ) ∝ (cid:18) . aν ) . + 0 . aν ) . (cid:19) (cid:114) aν π exp (cid:32) − . aν (cid:20) . aν ) . + 0 . aν ) . (cid:21) (cid:33) , – 6 –here a = 0 . The differential luminosity from annihilation in a dark matter halo of mass M is given by theproduct of a particle physics term and an astrophysical J -factor, i.e. the line-of-sight integralof the dark matter density squared, boosted by the annihilations in halo substructures. Thedark matter density can be parameterised using the same profiles that fit N-body simulationswell. For an NFW profile, the J factor of a point-like source can be analytically recast as J = (1 + B ) a ( c vir ) ρ s M vir , (3.5)in terms of the substructure boost B , the virial concentration c vir , the scale density ρ s of theprofile, and the analytic function (e.g., [45]) a ( c ) = (cid:18) − c ) (cid:19) (cid:18) ln(1 + c ) − c c (cid:19) − . The concentration parameter c of an NFW halo is related to the background density atthe time that the halo forms: small mass halos are more concentrated than high-mass halosbecause they form earlier (hierarchical halo formation). The concentration is also the linkbetween scale parameters ( ρ s , r s ) of the halo profile and the halo’s mass content [46]. The presence of halo substructures enhances the luminosity of the halo as a whole.Furthermore, halo substructure is expected to be denser than a host halo of the same mass(e.g., [45]). This higher density entails a larger number of annihilations, further enhancing theluminosity. The boost factor B parameterises this substructure luminosity as a proportionof the host luminosity. Since the J factor increases as the density squared, substructure isexpected to contribute between twice to twenty times as much luminosity as host structuresof mass 10 < M/M (cid:12) < [48].There are a number of N-body fitted models for NFW concentration and substructureboost to choose from in the literature. The models for B used in this study are listed inTable 1.The optimistic model is based on fit to the numerical simulations [49], which is wellmotivated for mass scales larger than the current resolution limit. For smaller scales, on theother hand, it heavily relies on validity of the phenomenological extrapolation. Reference [48]pointed out that such a power-law extrapolation was unphysical, and came up with morephysically motivated model by adoping the flattening of the concentration-mass relationfound for field halos. We convert between the virial and ( · · · ) conventions found in the literature using the prescription ofRef. [47]. This allows us to interchangeably convert between c vir and c at a given mass, and to convertbetween M vir and M given a choice of concentration-mass relation. – 7 –oost models [ref] FormulaNo boost B = 0Fiducial [48] Eq. (3.6)Optimistic [49] B = 1 . × − ( M /M (cid:12) ) . Table 1 : Substructure boost models considered in this study, ordered from least to mostoptimistic in terms of annihilation signal detection prospects. The nonzero models are derivedfrom fits to N-body simulations.The boost factor for our fiducial model is given as [48]:log B = (cid:88) i =0 B i log ( M ) i , (3.6) B i = (cid:8) − . , . , − . , . × − , . × − , − . × − (cid:9) . (3.7)The underlying concentration model is very close to other N-body-motivated models [50, 51].The third model represents the most conservative (but unexpected) situation in which thereis no substructure boost. In this study, we consider all the boost models in order to bracketuncertainties on subhalos, although the fiducial model is preferred over the others. The particle physics component K ( E ) = (cid:104) σv (cid:105) ( dN/dE ) /m χ of this model ( L ( E ) = J K ( E ))is assumed here to be a standard WIMP model, with annihilations into gamma rays via b ¯ b . (cid:104) σv (cid:105) is taken to be the thermal cross section 3 × − cm s − , and the WIMP mass is takento be m χ = 85 GeV. Our parameterisation of the photon number per energy per annihilatingparticle is [52]: dNdE = 0 .
42 exp( − x ) m χ ( x . + 0 . , x = Em χ . (3.8)With the values quoted above, we have K = 1 . × − cm s − MeV − at 1 GeV. Sinceour method is independent of any specific particle physics model, and depends only linearlyon K , it is trivial to rescale this study’s results to accommodate other particle physics models(as demonstrated in Sec. 5.3). By restricting our analysis to small enough redshifts and energies, we do not need to considerphotoionisation or pair production [53]. We can then use a very rough parameterisation [52]of gamma-ray absorbtion: e − τ ( E,z ) = exp (cid:34) − z . (cid:18) E (cid:19) . (cid:35) . (3.9)That this is very simplistic does not matter too much, since it is not a very important effectat low energies and redshifts anyway — as can be justified quantitatively by a sensitivityanalysis. – 8 – Results
In this section, we calculate P ( F ) and P ( F ) given the physical inputs of Sec. 3. We find that P ( F ) is roughly power-law like with a log-slope of approximately − P ( F ) is sensitive. More importantly,we find that P ( F ) takes the shape of a Gaussian at low flux, connecting smoothly with apower-law tail at high flux. This power-law tail, which has the same slope as P ( F ), is theregime where the flux is dominated by a singularly bright source. P ( F )The probability [Eq. (2.5)] that any single halo produces a differential flux observed with avalue F is then given by marginalising away the uncertainty in its mass and redshift. Thismarginalisation is clearly sensitive to the model choices we have presented above. We nowdiscuss the sensitivity of P ( F ) to the model choices. We find that the two most significant effects on P ( F ) (from amongst the effects we consid-ered; see below) are the choice of observing energy at the Fermi-LAT, and the substructureboost model B ( M ). These are represented in Fig. 1, which illustrates not only that P ( F )is more complex than just a power law P ( F ) ∝ F − (even in the absence of a substructureboost), but also that the substructure drastically influences the shape of P ( F ).The choice of the energy at which we study our differential flux represents a tradeoffbetween the amount of flux we expect to observe (which decreases with energy; Fig. 1), andthe instrumental capacity to actually observe it (which increases with energy; Ref. [1]). Theobserving energy of 1 GeV adopted in this study represents the lowest energy at which wecan leverage Fermi’s best angular resolution and effective area. The instrument’s energyresolution of ∆ E/E = 9% at 1 GeV is optimum for normally incident photons [1]. Anysystematic error due to misidentified energies is therefore minimised by this choice of energy(although we do not account for this instrumental effect in our analysis).Consequently, the effect of a gamma-ray optical depth (which predominantly affectsphotons from distant sources) only affects the low-flux region of P ( F ). This simply reflectsthe fact that at equal luminosity, distant sources produce less flux. At our low energy of1 GeV, the attenuation factor of a source at z = 3 is e − τ ∼ .
86: the net effect is definitelysmaller than a few percent after marginalising over the nearby halos, for which the attenuationis truly negligible. A better parameterisation than [52] was therefore not deemed necessaryfor this exploratory analysis.Since our differential flux is roughly F ∼ L/d L , there is an extent to which ‘bright’halos are massive and nearby, while ‘faint’ halos are light and distant. Therefore the largestmodeling factor in the high-flux tail is the choice of boost model. This is reflected in the rightpanel of Fig. 1: the optimistic model [49] gives very large high-flux tail, in comparison to thefiducial model. The mass scale at which the fiducial and optimistic boost models intersect( ∼ M (cid:12) [48]) is present in the flux distribution. We find that P ( F ) is less sensitive toother modeling choices. We have used various analytical models for the power spectrum transfer functions T ( k )[54, 55], and for c vir fits to N-body simulations [48, 50, 51]. The resulting distributions P ( F )– 9 –ntensity F/ Ω pix I P ( I ) [ c m − s − s r − M e V − ] − . − . . . . . − − − − − − − − [cm − s − sr − MeV − ]No BoostFiducialOptimistic10 − − − − − − − − Figure 1 : The flux/intensity PDF for a single dark matter source, P ( F ), with its depen-dence on photon energy (left) and boost models (right). In the right panel, the blue, black,and red curves represent respectively the pessimistic, fiducial, and optimistic models of thesubhalo boost. The choice of boost model clearly and significantly affects the functional formof the one-point function. The log-slopes of the fiducial model are offset (black dashed) andquantified for convenience. The flux F and intensity I of the gamma-ray background arerelated via the instrument’s pixel size: F = I Ω pix , where Ω pix ≈ . × − sr for E = 1 GeVphotons.were found to be robust to changes in these inputs. Indeed, the function a ( c ) in the J factor [Eq. (3.5)] is very smoothly decreasing for c >
1, so the choice of concentration modeldoes not influence the final result much. Similarly, σ ( M ) only changes by about an orderof magnitude over many orders of magnitude of halo masses, so varying models of T ( k ) (oreven using an N-body fit [56]) does not significantly change P ( F ) either.Using such a simple, ‘self-similar’ mass function, is not without shortcomings: theseMarkov-process models do not follow individual halo histories, and cannot account for dy-namical effects such as dynamical friction or tidal stripping. This, amongst other concerns,underlines the danger of our uncontrolled extrapolation of the mass function down to halomasses 10 − M (cid:12) . Nevertheless, mass function fits to N-body simulations are often very good;the fits even favour ellipsoidal collapse models [43, 44] (which we adopt in this study) overthe spherical collapse models of the seminal papers. The ellipsoidal collapse model is accu-rate to a few % over the large halo mass range (while spherical collapse [38] underestimatesthe number of FOF halos [57]) and gives an excellent fit in the range of 10 –10 M (cid:12) [58].Ellipsoidal collapse is also suitable for halos as small as 10 M (cid:12) and as early as z = 15 [59].We studied dependence of the mass function on Λ, in addition to ellipsoidal collapse [42].Without significantly changing P ( F ),we find that the delayed structure formation gives aroughly three times larger total flux (by increasing the number of halos N ) than if we hadused a mass function for which the cosmological constant is ignored [43, 44].We will need, when computing the total P ( F ), to compute the first few moments of– 10 – ( F ) at an intermediate step: E P ( F ) = (cid:90) dF (cid:48) P ( F (cid:48) ) F (cid:48) , (4.1) V P ( F ) = (cid:90) dF (cid:48) P ( F (cid:48) )( F (cid:48) − E P ( F ) ) . (4.2)After multiplying Eq. (4.1) by the mean number of halos [Eq. (3.1)] and dividing by thepixel size, one obtains the mean intensity of the gamma-ray background from dark matterannihilation [60]. Similarly the variance [Eq. (4.2)] is related to the angular power spectrumafter similar corrections [61].Since Eq. (2.7) entails an integration over redshift for each F , and since we find P ( F )to be relatively smooth, we calculate 250 logarithmically equidistant points over the ∼ θ ( M, z ). Of these, the largest quantifiable uncertainty is the scatter about theconcentration c , estimated at 15% [51]. P ( F ) Once P ( F ) and the number of halos per pixel N (cid:48) is specified, the calculation of P ( F ) requiresno additional physical assumptions. However, the large number of halos per pixel k ≈ N (cid:48) makes any exact calculation of the autoconvolution [Eq. (2.8)] prohibitively expensive, evenusing Fourier methods [25, 30, 62]. One might attempt to use the central limit theorem(CLT) to approximate P k ( F ) by a Gaussian. But, although the CLT guarantees convergencein distribution for k → ∞ , at finite k there will be deviations from a Normal distribution,particularly in the tails [63]. This is especially true for power-law-like distributions such as P ( F ), since the stable distributions of sums of power laws may be non-Gaussian [31]. Seealso Appendix A.For the purpose of solving the autoconvolution [Eq. (2.8)], we can choose to split P ( F )into two physically interesting contributions: the multitude of low-mass, faint halos con-tributing less than some cutoff flux F ∗ ; and the rare, large point sources brighter than F ∗ .This scheme is illustrated in Fig. 2. The former will contribute an isotropic backgroundwith Gaussian statistics, while the latter will add a power-law contribution to the high-fluxtail, which smoothly matches with P ( F ) at high fluxes (where flux is dominated by a sin-gle source). The noise at high flux in our Monte Carlo generated P ( F ) can therefore beconfidently ignored as a numerical artefact (cf. Fig. 10 in Appendix B).Despite the large value of k , our cutoff F ∗ may be chosen such that only a few thousandof these high-flux sources remain in each Fermi pixel. We can model the contribution of thesesources by Monte Carlo, drawing number of these rare sources from a Poisson distributionand their flux from P ( F ). We obtain finally the flux from k halos as the sum of these twocontributions; The flux distribution from k sources is given by the following convolution: P k ( F ) = G F 1. The angle arctan( R/d A ) then determines the number of pixelsover which the flux is averaged into an intensity. This corresponds to flux dilutions overroughly 10 pixels for Coma and 60 pixels for Fornax, explaining why intensities from Comaand Fornax appear inverted in lower two panels of Fig. 4: the total flux increases whenconsidering substructure, but flux per solid angle decreases more for Fornax than for Coma.The fact that the intensity from Fornax appears to decrease from the top panel to the lowerpanels of Fig. 4 is then just a manifestation of the difference between seeing Fornax as apoint source in the top panel or as an extended source in the lower panels.For the optimistic boost model, Coma stands out in the tail of P ( F ), while Fornaxis only barely more visible than if it (pessimistically) had no substructure. Although ourtreatment of source extension is somewhat naive, the diffuse gamma-ray background would– 15 –e a limiting factor in cluster analysis for the fiducial boost model, even if the effects ofextension were favourably revised by a factor of three or four (see middle panel of Fig. 4).The intrinsically poor signal-to-noise in cluster searches is also independent of the anni-hilation cross section or mass: changing these particle physics parameters would not changethe signal-to-noise ratio, since both the target cluster and the gamma-ray background arerescaled by the same factor. Dwarf spheroidal galaxies are another strategic choice for gamma-ray point-source searches,to which dark matter annihilations in the Milky Way substructures constitute a knownbackground [28]. We now discuss the background due to extragalactic sources.We consider Draco and Segue 1, with J -factors from Ref. [69] since we do not expect ourvirialised halo model (Sec. 3.2) to apply to them. Again, these sources have larger J -factorsthan other known dwarfs, and should thereby set the strongest constraints on non-detection.These sources lie well above the isotropic background component, even if we assume norelevant substructure boost for the dwarfs ( B = 0). Consequently, the more substructureboost there is in extragalactic dark matter halos, the worse the signal-to-noise for dwarfs willbe. Since we know the full distribution P ( F ) for the extragalactic background, the p -value ofan excess signal at 1 GeV (due to a dwarf spheroidal) can readily be estimated. Even thoughthe mean intensities yield poor signal-to-noise ratios, a pixel as bright as Segue 1 would berelatively uncommon (though not absent from the Fermi skymap, see Fig. 3). In this section we consider a number of known astrophysical backgrounds that would maskthe signal of our fiducial model. Following the parameterisation of the blazar source count from Ref. [4], we assume thefollowing power-law for P ( F ) of unresolved blazars: d NdSd Ω = 1 . × − (cid:18) S cm − s − (cid:19) − . cm s deg − , (5.3) P ( F ) = 2 . × − (cid:18) F cm − s − MeV − (cid:19) − . cm s MeV , (5.4)where S is the gamma-ray flux integrated above 100 MeV. We extrapolate this relation downto a lower flux limit of S ≥ . × − cm − s − , and perturb around this fiducial value to test the model’s sensitivity to this extrapolation (see Table 3). We also enforce an upperlimit of S ≤ × − cm − s − , in order to maintain a low blazar detection efficiency [4]without altering the power-law form of the blazar P ( F ). Assuming an E − . spectrum, thecorresponding range in differential flux at 1 GeV is then 2 . × − ≤ F/ (cm − s − MeV − ) ≤ . × − .Using d N/dSd Ω, we find ∼ . P ( F ) by Monte Carlo simulation. Inpractice, we split the blazar P ( F ) into two contributions: the delta-function of zero sources This value is ten times smaller than the ‘medium band’ from Ref. [4], corresponding to the faintest observed source out of the mix of FSRQs and BL Lacs that contributes at 1 GeV. – 16 –ntensity F/ Ω pix [cm − s − sr − MeV − ] F P ( F ) − − − − − Figure 5 : One-point function P ( F ) for the three dark matter models (with boosts color-coded as previously), alongside the P ( F ) of the diffuse contribution of blazars (green). Thedashed red band represents the measurement of the unresolved EGB from the Fermi data at1 GeV [65], while the dashed green line is the mean of the blazar PDF. S min N blz/pix (cid:104) F (cid:105) / Ω pix (% EGB)0.72 1.25 0.89 (14.8%)0.36 1.97 0.921 (15.3%)0.18 3.08 0.946 (15.8%) Table 3 : Sensitivity summary of our unresolved blazar model (with fiducial values in thecentral row). The first two columns pertain to Eqn. (5.3): S min (in units of 10 − cm − s − )is the lowest flux to which we extrapolate the source count distribution, from which a numberof (faint, unresolved) blazars per pixel may be derived. The next two columns summarisethe corresponding mean intensity of the blazar P ( F ), both as an absolute value in units of10 − cm − s − sr − MeV − , and as a proportion of the unresolved EGB [65].per pixel (treated analytically) and the contribution of more than zero sources per pixel(Monte Carlo): P Blazar ( F ) = ∞ (cid:88) k =0 P ( k | . P k ( F ) = e − . δ ( F ) + ∞ (cid:88) k =1 P ( k | . P k ( F ) . (5.5)The (nonzero) blazar flux distribution is plotted in Fig. 5. Like the flux distribution for thedark matter, it has a non-negligible skew and approximates the single-source distribution at– 17 –igh fluxes. P Blazar ( F ) also follows the single-source powerlaw between the minimal intensityof a single source I min = F min / Ω pix and the minimal intensity of two sources, with a discon-tinuous derivative at this threshold. We see from Fig. 5 and Table 3 that unresolved blazarscontribute roughly 15% of the unresolved EGB at 1 GeV, reproducing Ref. [4]: even a singleblazar is at least as bright as the entire dark matter contribution of our fiducial model, andthe mean blazar flux is between one and two orders of magnitude brighter than the darkmatter, depending on the boost model. Fortunately, using the entire distributions (insteadof just their means) allows the thin, peaked dark matter to be statistically extracted fromthe broad, powerlawlike blazars (see Sec. 5.3). There are a number of other backgrounds to consider, such as cosmic rays or starburstgalaxies. Since we expect all contributions of the EGB to add up to the experimentallyobserved flux, we must convolve the distributions above with the P ( F ) of these other isotropicbackgrounds. However, a complete model of these backgrounds is beyond the scope of thisexploratory analysis. The one-point function of these other isotropic components is (forconvenience) assumed Gaussian with a small, arbitrarily chosen, but non-negligible variance( µ/σ ∼ ) and mean determined by requiring that the mean intensity due to P EGB ( F ) (afteradding the dark matter and blazar components) be equal to the experimentally determinedvalue 6 × − cm − s − sr − MeV − [65]. This Gaussian could be thought of as the centrallimit theorem approximation to the flux distribution of a large number of sources, since µ/σ scales like √ N .We can also consider not convolving a dark matter component into the gamma-raybackground, to study (in Sec. 5.3) the distinguishability of an alternate hypothesis modelwith dark matter, P Alt ( F ), from a null hypothesis model without dark matter P Null ( F )(Both of these flux distributions are represented in Fig. 6). We then have P EGB ( F ) = P Blazar ( F ) (cid:63) G rest ( F ) (cid:63) P DM ( F ) , (5.6)which can be solved using the blazar P ( F ) model [Eq. (5.5)] above, as P EGB ( F ) = (cid:2) e − . δ ( F ) + (1 − e − . ) P MCBlazar ( F ) (cid:3) (cid:63) G rest ( F ) (cid:63) P DM ( F ) , (5.7)where the analytical factor of (1 − e − . ) captures the normalisation of the Monte Carlo. Itshall be crucial in what follows that the mean of P EGB ( F ) has been fixed as an experimentalinput, rather than allowed to vary freely as a model parameter. The purpose of this section is to describe the influence of modelling choices on the shape ofthe flux distribution of the total unresolved EGB, as represented in Fig. 6.The thin peaks and power-law-tailed peaks at low and high flux in Fig. 6, correspondrespectively to whether or not a blazar is present in the associated pixel. The relative peakheights are determined (to a first approximation) by the number of unresolved blazars perpixel, via the mixture coefficient e − . . Broadening the peaks (by convolving the astrophysicalcomponents with a dark matter component) introduces a correction to this determinationof height. The location of the mean of the low-flux peak is set via (i) Eqn. (5.6), (ii) the(assumed) mean of the blazar flux distribution (cf. Table 3), and (iii) the experimentallydetermined mean of the total flux distribution [65]. The difference between the dark matter– 18 –ntensity F/ Ω pix × [cm − s − sr − MeV − ] F P ( F ) − − Figure 6 : Predicted flux distribution P EGB ( F ) of the extragalactic gamma-ray background,with (black) and without (green) a contribution from dark matter annihilations. The distri-butions have two peaks, based on whether or not a blazar is present in the associated pixel.The mean EGB derived from Fermi [65] is represented by the vertical line (red, dashed). Across-section twice the canonical value was used to visually enhance the differences betweenthese distributions.distribution’s mean and most probable fluxes (cf. Table 2) introduces an additional percent-level shift between the locations of the low-flux peaks with and without dark matter.The width and depth of the ‘gap’ between the high and low flux peaks is related, inthe absence of dark matter, to the lower flux limit on blazars assumed in Sec. 5.2.1. Thisrelation relies on the interplay between two effects: on one hand, extrapolating to lowerfluxes increases the mean blazar contribution, shifting the location of the low-flux peak by anequal and opposite amount to even lower fluxes and widening the gap. On the other hand,extrapolating the faintest blazar contribution to lower fluxes also fills the gap with the fluxesfrom these faint blazars. The latter effect more than compensates for the former, such thatoverall the gap closes as it shifts to lower fluxes with an increasing blazar contribution.However, the contribution of the dark matter component to this gap is just as dramatic.Since the distribution is skewed, it broadens both peaks preferentially to higher fluxes; how-ever, the location of the low-flux peak is fixed by the constraint imposed on the mean of P EGB ( F ), while the high-flux peak is free to shift under this broadening. This widens thegap with increasing brightness of the dark matter annihilation signal.In addition to these shifts, the slope of the flux-tail due to dark matter annihilation( − . 5) is quite different from the one for blazars ( − . ± × − cm − s − sr − MeV − , Ref. [65]), and may be open to exploitation byone-point function methods.As the sum of a diffuse background and a few unresolved sources, the gamma-ray back-ground in Eq. (5.6) has the same origin as the dark matter background. By adopting awider variance for the Gaussian in Eq. (5.6), even the ‘Gaussian with a power-law tail’ formof the dark matter component can even be reproduced, corroborating our interpretation ofthese features of the dark matter P ( F ) and justifying our faint/bright F ∗ analysis. Onemight worry that this ‘truly diffuse plus unresolved point sources’ form also undermines theprospects of an unambiguous dark matter detection above such a background. Indeed, model-fitting one-point techniques such as that described in Ref. [24] cannot separate degeneratemodels, especially given the angular resolution limitation (illustrated in Fig. 3) that preventsthe dark matter P ( F ) power-law tail from contributing significantly to the observed flux.However, one need not worry too much: the asymmetry-induced shift of the peaks of P ( F )may be a sufficiently distinctive feature to extract a dark matter signal from the Fermi datanonetheless. P ( C )The observable given by Fermi is not the gamma-ray flux F , but the discrete number of photoncounts per pixel C . Photon arrival may then be modelled as a Poisson rate with a meandetermined by the differential gamma-ray flux and the exposure (cid:15) = (time) × (detector area) × (photon energy). For a five-year Fermi mission, correcting for the field of view, we have anexposure of (cid:15) ≈ . × cm s MeV sr pixel − . Marginalising over the uncertain fluxdistribution then gives P ( C ) = (cid:90) P EGB ( F ) P ( C | (cid:15)F ) dF. (5.8)This Poisson arrival uncertainty substantially smooths away the differences between the nulland alternate flux models, as evidenced by Fig. 7. However, the percent-level shift betweenthe low-flux peaks due to the dark matter distribution’s skewness survives, since a percentdifference with C ∼ O (100) is still a few photons. There is also a larger, opposite shift in thepoint-source-driven high-flux tail due to our imposed value of the distribution’s mean.We can define, given our number of pixels N pix , the test statistic [25] χ = (cid:88) C (cid:32) N pix [ P Null ( C ) − P Alt ( C )] (cid:112) N pix P Null ( C ) (cid:33) , (5.9)The choice of bounds for the sum over count bins is somewhat arbitrary; a formal optimisationof this test statistic would be beyond the scope of this analysis. We choose to focus on thepeak 65 < C < 165 of the distribution, in which we anticipate sufficiently many pixels per‘count bin’ to trust a χ test. The lower panel of Fig. 7 illustrates the terms of the sum inthis test statistic, from which one may obtain the p -value at which data with dark matterfollowing P Alt ( C ) exactly would reject our dark-matter-free null hypothesis.– 20 –C N p i x P ( C ) [ p i x e l ] 60 80 100 120 140 160-100-50050 Photon Count C [pixel − ] N p i x ∆ P ( C ) [ p i x e l ] Figure 7 : Predicted count distributions of EBG photons with (black) and without (green)a dark matter component. The green bands represent the Poisson errors σ ∝ (cid:112) N P ( C ) onthe dark-matter-free model. The lower panel shows difference between the two models. Across-section twice the canonical value was used to visually enhance the differences betweenthese distributions.For the fiducial dark matter model, the percent-level shift between these peaks is justsmall enough that the null cannot be rejected by the data. Since the dark matter cross section (cid:104) σv (cid:105) enters in our model only as a proportionality factor for the dark matter halo luminosity,we can rescale our P DM ( F ) to inexpensively repeat this detectability study for higher valuesof the cross section. We can then forecast for which values of the cross section a dark mattercomponent would become distinguishable from the background using the one-point functionalone. This is summarised in Fig. 8, which shows that (given a perfect understanding of thebackgrounds) the one-point function could probe a dark matter annihilation signal with across section roughly a factor two times larger than the canonical value 3 × − cm s − .Fig. 8 also shows that our fiducial choice of lower blazar flux extrapolation is fortuitouslyclose to a detectability optimum, but our forecast does not deteriorate much upon rescalingthis value.This result, complementary to two-point function analyses, could even be strengthened– 21 –.8 1.9 2.0 2.1 2.2 2.3 2.4Cross Section (cid:104) σv (cid:105) / × − cm s − p v a l u e S min = 0 . 36 0 . . 18 2 σ σ σ σ − − − − − − − Figure 8 : Predicted statistical significance of a (hypothetical) one-point-function-only de-tection of a dark matter annihilation signal above perfectly characterised astrophysical back-grounds, as a function of the dark matter cross section. Curves are labelled by the flux S min down to which the blazar distribution is extrapolated (see Table 3). Horizontal lines (blue,dashed) represent some common choices of confidence level. Including the energy-dependenceof the flux distributions would improve these results, at the cost of a greater dependence onthe annihilation spectrum.by including the energy dependence of the differential flux to break the degeneracy with theastrophysical backgrounds [28]. Such a study would remain sensitive to (but would allow aquantitative analysis of) the assumptions and uncertainties of the astrophysical backgroundmodel. Yet, even without this spectral input, our forecasted one-point upper limit on thecross-section is on par with the most recent (spectral) constraints [5] based on the mean valuealone.In addition to the extragalactic dark matter flux, there will be a component due toGalactic substructures. The one-point distribution of such a Galactic component has beenpredicted [25], and similarly features a power-law high-flux tail. Due to the energy spec-trum Eq. (3.8), if the mean intensity from subhalos at the anticenter integrated above10 GeV is ∼ − cm − s − sr − [25], then the mean differential intensity at 1 GeV is (cid:104) I (cid:105) ∼ − cm − s − sr − MeV − . This is of the same order of magnitude as the extragalac-tic component discussed above, and with the same high-flux F − . power-law tail. Thus,including the Galactic component would further enhance the expected signal-to-noise for po-tential detection. We finally note that the Galactic component will show a dipole feature,with more flux from the Galactic center than the anticenter, which can in principle be usedto discriminate it from the isotropic extragalactic component. There are a number of caveats on the results presented in this study. Firstly, a large numberof assumptions were used to simplify the hierarchical model without a proper sensitivityanalysis. We assumed there is no scatter in the halo parameters [Eq. (2.6)], or any uncertainty– 22 –n the average number N (cid:48) of halos in each pixel. An NFW profile was assumed for thehalos despite the fact that Einasto and uncusped profiles would probably give less flux. Theuncertainty incurred by extrapolating the mass function dN/dM and the boost models B ( M )down to M min , is compounded by our ignorance of the actual value of M min .Secondly, we have not studied how our results depend on pixel size, particularly theeffects of source extension. We have merely assumed (Sec. 2.2) that all sources are point-like,since there are on average only 0.28 extended dark matter sources per pixel, a negligiblefraction of all N (cid:48) = 7 × halos (Sec. 3.1.1). However, extended sources must be eithermassive or nearby (Sec. 4.1), and therefore would tend to have large fluxes, affecting thedistribution P F >F ∗ ( F ). Our point-source-based P ( F ) is therefore not applicable to theseobjects at high flux. That this compromises the analysis should be obvious from Fig. 4: theone-point functions do not account for the clusters’ extension, and the clusters’ fluxes clearlydo not live in the domain of the PDF they should be drawn from. Thus, extended sourcesshould be dealt with in a complete analysis. Some of the elements of such an analysis (suchas an energy-dependent angular resolution, a substructure-boost-dependent mass-to-solid-angle conversion, or a redshift-dependent mass threshold M Ext ( z )) have been presented inthe main text. Note, however, that the distribution of faint and distant sources P F We constructed a hierarchical model that predicts, using analytical models of ΛCDM struc-ture formation, the flux distribution of gamma rays from extragalactic dark matter annihi-lation in unresolved point sources. The uncertainties on this flux subject to the modelingchoices we studied are typically percent-level; in the case of the substructure boost function,they remain smaller than a factor of three. We then compute, without requiring any addi-tional physical assumptions, the flux distribution per pixel P ( F ), which has the characteristicform of an isotropic diffuse Gaussian matched at high flux to the point-source distributionwith a power-law slope of − . 5. This distribution is non-Gaussian and asymmetric; however– 23 –he most likely flux and the mean flux are comparable at the percent-level in all but theoptimistic boost model, salvaging previous ‘mean intensity’ constraints on the dark matterproperties from this potential systematic effect.The fluxes predicted for our fiducial model lie just within the reach of the Fermi-LAT,and should be observable by the tenth year of the mission. We also showed that the dis-tinctive features of the power-law-tailed Gaussian distribution all live above Fermi’s angularresolution. Therefore, the extragalactic gamma-ray emission due to dark matter annihilationconstitutes an irreducible and significant background for point-source annihilation searcheswith clusters or dwarf spheroidals. Ironically, an optimistic boost model would be detrimen-tal to these searches, by deteriorating the signal-to-noise of these point sources (to unity orworse for galaxy clusters).We also discussed the astrophysical backgrounds from which a dark matter annihila-tion signal would need to be extracted. These include unresolved blazars (which contributean order of magnitude more flux than the fiducial dark matter model) and other diffusecomponents, which were all convolved together into a total model for the gamma-ray back-ground. The scarcity of unresolved blazars make this distribution quite rich in features; mostprominently, it has two distinct peaks of most probable fluxes, the inter-peak gap being verysensitive to the dark matter component.Even accounting for the Poisson noise of photon arrivals that come with such low fluxes,a contribution to the gamma-ray background of the order of a vanilla WIMP model may be de-tectable above well-characterised astrophysical backgrounds using the flux distribution alone.Using the energy-dependence of the flux distribution should further break the degeneracy be-tween the components of the gamma-ray background, and should allow one-point functionmethods to complement and strengthen existing constraints set by two-point-function anal-yses. Acknowledgments This work was supported by the Netherlands Organisation for Scientific Research throughVidi grant (MRF and SA). A Failure of the central limit theorem When constructing P ( F ) from many individual sources (each an independent and identicallydistributed realisation of P ( F )), it is easy to show that we expect significant deviations froma Normal distribution, even for the large number N (cid:48) ≈ × of sources per pixel. Theratio between the mean flux and the standard deviation of the fiducial model is approximately( µ/σ ) N (cid:48) = 0 . 23, which (using the Gaussian’s associated cumulative distribution) would firmlyplace 40% of a Gaussian P ( F ) at negative fluxes. This is not only mathematically impossible(given that the sum of positive random variables must be positive), it is also physicallynonsensical. Note that using distributions with widths instead of delta functions for P ( N (cid:48) )or P ( k | N (cid:48) ) makes this problem even worse by broadening the unphysical ‘Gaussian’ P k .We can also show that we expect deviations from a Gaussian with slightly more rigour:The k -autoconvolution definition of P k ( F ) = ( P ) (cid:63)k ( F ) prompts us to work in the Fourierspace of probability distributions, i.e., with the characteristic functions φ ( t ) and φ k = ( φ ) k .– 24 –arginalising P k ( F ) with the Poisson distribution P ( k | N (cid:48) ) gives [30]: P ( F | N (cid:48) ) = ∞ (cid:88) k =0 e − N (cid:48) N (cid:48) k k ! F − (cid:104) φ k (cid:105) (A.1)= F − (cid:34) e − N (cid:48) ∞ (cid:88) k =0 ( µφ ) k k ! (cid:35) = F − (cid:104) e N (cid:48) ( φ − (cid:105) , where we have used linearity of the (inverse) Fourier transform and have recognised the powerseries expansion of the exponential function. By Taylor expanding the characteristic function φ , we generate (by construction) the first moments of P : φ = 1 + E ( P )( it ) + V ( P ) ( it ) − it ) , (A.2)such that the characteristic function associated to P ( F | µ ) becomes (to second order in t ) φ P ( F | N (cid:48) ) = exp (cid:18)(cid:2) N (cid:48) E ( P ) (cid:3) ( it ) − (cid:2) N (cid:48) V ( P ) (cid:3) t (cid:19) e N (cid:48) o( − it ) . (A.3)We recognise this first exponential as the characteristic function of a Normal distribution, so P ( F | N (cid:48) ) ≈ G (cid:0) F | N (cid:48) E ( P ) , N (cid:48) V ( P ) (cid:1) (cid:63) F − (cid:104) e N (cid:48) o( − it ) (cid:105) . (A.4)If P ( N (cid:48) ) is a δ -function we recover upon marginalisation the near-Gaussian form of P ( F )anticipated in the main text, with explicit deviations due to higher moments convolved in.Proving that these deviations lead to a power-law tail at high flux is of less interest thanthe observation, in the main text, that in this regime the flux is almost completely due to asingle bright source. This leads, in the next section, to a Monte Monte CarloCarlo analysisof these corrections.The analysis above does not assume N (cid:48) , k → ∞ . This reflects the fact that we have not derived the Central limit Theorem (CLT), the ‘Gaussian (cid:63) Corrections’ form was derived fora finite number of sources. However, the CLT is still valid asymptotically as k → ∞ since P does have a finite variance despite its power-law-like behaviour. B Our method B.1 Rationale of our method The deviations from Gaussianity in the sum P k of random variables are bounded by theBerry-Esseen (BE) theorem [31, 63], which relates the Kolmogorov-Smirnov distance (thelargest deviation at any point), to the ‘absolute skewness’ of P ( F ) and the large but finitenumber of dark matter halos k :sup F | CDF k ( F ) − Φ( F ) | ≤ Cρ √ k , ρ = E ( | X − µ | ) σ , (B.1) Formally, ρ is the sum of the partial skewnesses above and below the mean. We simply call ρ theskewness. – 25 –here C ∼ . 5, and where CDF k and Φ( F ) are the cumulative distribution for P k and theNormal distribution, respectively. We note that the absolute value in ρ distinguishes it fromthe conventional skewness γ , and gives the inequality ρ ≥ γ . For our power-law-like P ( F ),the skewness is huge and the BE bound Cρ k − / is uninformative: the Kolmogorov-Smirnovdistance (a quantity definitionally less than one), is ‘constrained’ by the first few momentsof our fiducial P ( F ) to be less that about 10 .This perspective suggests splitting P ( F ) into low-flux and high-flux contributions, toreduce the skewness of each contribution. In the following, we (i) derive and (ii) cross-checkthe Monte Carlo method presented in the main text. B.2 Deriving the form of P k ( F )We derive the behaviour for P k as follows: marginalise P ( F ) into high and low flux contri-butions, P ( F ) = (1 − (cid:15) ) P F We should, of course, check that our choice of F ∗ is physically sensible. If we want to split P k ( F ) into an isotropic background and candidate point sources, we must verify that ourfaint sources do not contribute more than zero or one photons each. Recall that a singleGeV photon per pixel in 5 years of Fermi data corresponds to a differential flux/intensity of F = 6 × − m − yr − MeV − pixel − . Choosing F ∗ ≈ . × − m − yr − MeV − pixel − ,over two orders of magnitude smaller, guarantees that faint sources are indeed faint. Thisvalue of F ∗ was in fact chosen algorithmically, such that the BE bound for P F We can compute the probability of seeing a dark matter signal above a given flux, by lookingat the exceedance (complementary cumulative) distribution Ψ( F ) associated to P ( F ):Ψ( F max ) = (cid:90) ∞ F max P ( F ) dF. The probability of not realising this high flux tail in N pixel trials is then given by the binomialdistribution, and we want to solve the following for F max : B (0 | N pixel , Ψ( F max )) ≤ α , α = { . , . , . } . (C.1)– 27 – ( F ) F ∗ = 1 . × − F ∗ = 4 . × − F ∗ = 9 . × − F/ Ω pix × [cm − s − sr − MeV − ] r e l a t i v ee rr o r ( % ) Figure 9 : Local sensitivity analysis of our Monte Carlo on the cutoff flux F ∗ , with 10 samples of P F >F ∗ k ∗ for each F ∗ . The values of the cutoff F ∗ quoted in the legend are in unitsof m − yr − MeV − pixel − . The error in the lower panel is relative to the average of thethree Monte Carlo simulations.For k = 0 successes and p = Ψ( F max ) (cid:28) 1, we can expand this in a binomial series that wecan truncate at first order when Ψ( F max ) (cid:28) N − :1 − N pixel Ψ( F max ) + O ( N pixel Ψ( F max )) ≤ α. (C.2)When F max lies in the power-law tail of P ( F ), we can use the high-flux-tail equivalence(Fig. 10) of P ( F ) and P ( F ) to find, approximately, (cid:90) ∞ F max P ( F ) dF ≤ − αN pixel . (C.3)For a power-law-like P ( F ) ≈ AF γ , γ < (cid:90) ∞ F max AF γ dF = − A ( F max ) γ +1 / ( γ + 1) = − F max P ( F max ) / ( γ + 1) . (C.4)– 28 –ntensity F/ Ω pix [cm − s − sr − MeV − ] P r o b a b ili t y d e n s i t y ( unn o r m a li s e d ) − − − − − − − − − Figure 10 : Unnormalised P ( F ) and P ( F ) (solid, black) for the fiducial model. The scalingfactor to achieve the matching between the two distributions is not a free parameter, it isthe number of halos k − k ∗ = 1864 contributing to the high-flux tail in the fiducial model.We show also the (unnormalised) Gaussian with equivalent mean and variance predicted bya naive CLT (green, dashed).In the limit that ( α, γ ) → (0 , − . 5) we reproduce Eq. (4.4) from the main text. Moregenerally, for α → γ = − . δγ , we have F max P ( F max ) ≤ . N pixel − δγN pixel . (C.5)Since (as visible in Fig 1) the highest fluxes in the tail of P ( F ) have a steeper log-slope than − . δγ < References [1] Fermi-LAT Collaboration, W. Atwood et al., The Large Area Telescope on the FermiGamma-ray Space Telescope Mission , Astrophys. J. (2009) 1071–1102, [ arXiv:0902.1089 ].[2] Fermi-LAT Collaboration, M. Ackermann et al., The spectrum of isotropic diffuse gamma-rayemission between 100 MeV and 820 GeV , Astrophys. J. (2015), no. 1 86,[ arXiv:1410.3696 ].[3] Fermi-LAT Collaboration, M. Ackermann et al., Anisotropies in the diffuse gamma-raybackground measured by the Fermi-LAT , Phys. Rev. D (2012) 083007, [ arXiv:1202.2856 ]. – 29 – Fermi-LAT Collaboration, The Fermi-LAT high-latitude Survey: Source Count Distributionsand the Origin of the Extragalactic Diffuse Background , Astrophys. J. (2010) 435–453,[ arXiv:1003.0895 ].[5] M. Ajello, D. Gasparrini, M. S´anchez-Conde, G. Zaharijas, M. Gustafsson, et al., The Origin ofthe Extragalactic Gamma-Ray Background and Implications for Dark-Matter Annihilation , Astrophys. J. (2015), no. 2 L27, [ arXiv:1501.05301 ].[6] I. Tamborra, S. Ando, and K. Murase, Star-forming galaxies as the origin of diffusehigh-energy backgrounds: Gamma-ray and neutrino connections, and implications for starbursthistory , JCAP (2014), no. 09 043, [ arXiv:1404.1189 ].[7] M. Di Mauro, F. Calore, F. Donato, M. Ajello, and L. Latronico, Diffuse γ -ray emission frommisaligned active galactic nuclei , Astrophys. J. (2014) 161, [ arXiv:1304.0908 ].[8] M. Di Mauro and F. Donato, The composition of the Fermi-LAT IGRB intensity: emissionfrom extragalactic point sources and dark matter annihilations , arXiv:1501.05316 .[9] A. Cuoco, E. Komatsu, and J. Siegal-Gaskins, Joint anisotropy and source count constraints onthe contribution of blazars to the diffuse gamma-ray background , Phys. Rev. D86 (2012)063004, [ arXiv:1202.5309 ].[10] G. Jungman, M. Kamionkowski, and K. Griest, Supersymmetric dark matter , Phys. Rept. (1996) 195–373, [ hep-ph/9506380 ].[11] D. Hooper and S. Profumo, Dark matter and collider phenomenology of universal extradimensions , Phys. Rept. (2007) 29–115, [ hep-ph/0701197 ].[12] M. Fornasa and M. A. Sanchez-Conde, The nature of the Diffuse Gamma-Ray Background , arXiv:1502.02866 .[13] S. Ando and K. Ishiwata, Constraints on decaying dark matter from the extragalacticgamma-ray background , arXiv:1502.02007 .[14] S. Ando, A. Benoit-L´evy, and E. Komatsu, Mapping dark matter in the gamma-ray sky withgalaxy catalogs , Phys. Rev. D (2014), no. 2 023514, [ arXiv:1312.4403 ].[15] S. Ando, Power spectrum tomography of dark matter annihilation with local galaxy distribution , JCAP (2014), no. 10 061, [ arXiv:1407.8502 ].[16] S. Camera, M. Fornasa, N. Fornengo, and M. Regis, A Novel Approach in the WeaklyInteracting Massive Particle Quest: Cross-correlation of Gamma-Ray Anisotropies and CosmicShear , Astrophys. J. (2013) L5, [ arXiv:1212.5018 ].[17] S. Camera, M. Fornasa, N. Fornengo, and M. Regis, Tomographic-spectral approach for darkmatter detection in the cross-correlation between cosmic shear and diffuse gamma-ray emission , arXiv:1411.4651 .[18] J.-Q. Xia, A. Cuoco, E. Branchini, M. Fornasa, and M. Viel, A cross-correlation study of theFermi-LAT γ -ray diffuse extragalactic signal , Mon.Not.Roy.Astron.Soc. (2011) 2247–2264,[ arXiv:1103.4861 ].[19] M. Shirasaki, S. Horiuchi, and N. Yoshida, Cross-Correlation of Cosmic Shear andExtragalactic Gamma-ray Background: Constraints on the Dark Matter AnnihilationCross-Section , Phys. Rev. D (2014), no. 6 063502, [ arXiv:1404.5503 ].[20] N. Fornengo, L. Perotto, M. Regis, and S. Camera, Evidence of Cross-correlation between theCMB Lensing and the Γ -ray sky , Astrophys. J. (2015), no. 1 L1, [ arXiv:1410.4997 ].[21] J.-Q. Xia, A. Cuoco, E. Branchini, and M. Viel, Tomography of the Fermi-LAT gamma-raydiffuse extragalactic signal via cross-correlations with galaxy catalogs , Astrophys. J. Suppl. (2015), no. 1 15, [ arXiv:1503.05918 ]. – 30 – 22] M. Regis, J.-Q. Xia, A. Cuoco, E. Branchini, N. Fornengo, et al., Particle dark matter searchesoutside the Local neighborhood , arXiv:1503.05922 .[23] A. Berlin and D. Hooper, Stringent Constraints on the Dark Matter Annihilation Cross SectionFrom Subhalo Searches with the Fermi Gamma-Ray Space Telescope , Phys. Rev. D (2014),no. 1 016014, [ arXiv:1309.0525 ].[24] D. Malyshev and D. W. Hogg, Statistics of gamma-ray point sources below the Fermi detectionlimit , Astrophys. J. (2011) 181, [ arXiv:1104.0010 ].[25] S. K. Lee, S. Ando, and M. Kamionkowski, The Gamma-Ray-Flux Probability DistributionFunction from Galactic Halo Substructure , JCAP (2009) 007, [ arXiv:0810.1284 ].[26] S. Dodelson, A. V. Belikov, D. Hooper, and P. Serpico, Identifying Dark Matter AnnihilationProducts In The Diffuse Gamma Ray Background , Phys. Rev. D (2009) 083504,[ arXiv:0903.2829 ].[27] E. J. Baxter, S. Dodelson, S. M. Koushiappas, and L. E. Strigari, Constraining Dark Matter inGalactic Substructure , Phys. Rev. D (2010) 123511, [ arXiv:1006.2399 ].[28] E. Carlson, D. Hooper, and T. Linden, Improving the Sensitivity of Gamma-Ray Telescopes toDark Matter Annihilation in Dwarf Spheroidal Galaxies , Phys.Rev. D91 (2015), no. 6 061302,[ arXiv:1409.1572 ].[29] J. F. Navarro, C. S. Frenk, and S. D. White, A Universal density profile from hierarchicalclustering , Astrophys. J. (1997) 493–508, [ astro-ph/9611107 ].[30] P. A. Scheuer, A statistical method for analysing observations of faint radio stars , in Mathematical Proceedings of the Cambridge Philosophical Society , vol. 53, pp. 764–773,Cambridge Univ Press, 1957.[31] V. V. Petrov, Sums of independent random variables , vol. 82. Springer, 1975.[32] Planck Collaboration, P. Ade et al., Planck 2013 results. XVI. Cosmological parameters , Astron. Astrophys. (2014) A16, [ arXiv:1303.5076 ].[33] A. M. Green, S. Hofmann, and D. J. Schwarz, The First wimpy halos , JCAP (2005) 003,[ astro-ph/0503387 ].[34] A. Loeb and M. Zaldarriaga, The Small-scale power spectrum of cold dark matter , Phys.Rev. D71 (2005) 103520, [ astro-ph/0504112 ].[35] E. Bertschinger, The Effects of Cold Dark Matter Decoupling and Pair Annihilation onCosmological Perturbations , Phys.Rev. D74 (2006) 063509, [ astro-ph/0607319 ].[36] S. Hofmann, D. J. Schwarz, and H. Stoecker, Damping scales of neutralino cold dark matter , Phys.Rev. D64 (2001) 083507, [ astro-ph/0104173 ].[37] S. Profumo, K. Sigurdson, and M. Kamionkowski, What mass are the smallest protohalos? , Phys.Rev.Lett. (2006) 031301, [ astro-ph/0603373 ].[38] W. H. Press and P. Schechter, Formation of galaxies and clusters of galaxies by selfsimilargravitational condensation , Astrophys. J. (1974) 425–438.[39] J. Bond, S. Cole, G. Efstathiou, and N. Kaiser, Excursion set mass functions for hierarchicalGaussian fluctuations , Astrophys. J. (1991) 440.[40] J. M. Bardeen, J. Bond, N. Kaiser, and A. Szalay, The statistics of peaks of gaussian randomfields , The Astrophysical Journal (1986) 15–61.[41] A. Lapi, P. Salucci, and L. Danese, Statistics of Dark Matter Halos from the Excursion SetApproach , Astrophys. J. (2013) 85, [ arXiv:1305.7382 ].[42] A. Del Popolo, On the cosmological mass function theory , Astron. Rep. (2007) 709–734,[ astro-ph/0609166 ]. – 31 – 43] R. K. Sheth, H. Mo, and G. Tormen, Ellipsoidal collapse and an improved model for thenumber and spatial distribution of dark matter haloes , Mon. Not. Roy. Astron. Soc. (2001)1, [ astro-ph/9907024 ].[44] R. K. Sheth and G. Tormen, An Excursion set model of hierarchical clustering : Ellipsoidalcollapse and the moving barrier , Mon. Not. Roy. Astron. Soc. (2002) 61,[ astro-ph/0105113 ].[45] S. Ando, Gamma-ray background anisotropy from galactic dark matter substructure , Phys. Rev.D (2009) 023520, [ arXiv:0903.4685 ].[46] J. S. Bullock, T. S. Kolatt, Y. Sigad, R. S. Somerville, A. V. Kravtsov, et al., Profiles of darkhaloes. Evolution, scatter, and environment , Mon. Not. Roy. Astron. Soc. (2001) 559–575,[ astro-ph/9908159 ].[47] W. Hu and A. V. Kravtsov, Sample variance considerations for cluster surveys , Astrophys. J. (2003) 702–715, [ astro-ph/0203169 ].[48] M. A. Sanchez-Conde and F. Prada, The flattening of the concentration-mass relation towardslow halo masses and its implications for the annihilation signal boost , Mon. Not. Roy. Astron.Soc. (2014) 2271, [ arXiv:1312.1729 ].[49] L. Gao, C. Frenk, A. Jenkins, V. Springel, and S. White, Where will supersymmetric darkmatter first be seen? , Mon. Not. Roy. Astron. Soc. (2012) 1721, [ arXiv:1107.1916 ].[50] F. Prada, A. A. Klypin, A. J. Cuesta, J. E. Betancort-Rijo, and J. Primack, Haloconcentrations in the standard LCDM cosmology , Mon. Not. Roy. Astron. Soc. (2012)3018–3030, [ arXiv:1104.5130 ].[51] B. Diemer and A. V. Kravtsov, A universal model for halo concentrations , Astrophys. J. (2015), no. 1 108, [ arXiv:1407.4730 ].[52] L. Bergstrom, J. Edsjo, and P. Ullio, Spectral gamma-ray signatures of cosmological darkmatter annihilation , Phys. Rev. Lett. (2001) 251301, [ astro-ph/0105048 ].[53] A. A. Zdziarski and R. Svensson, Absorption of x-rays and gamma rays at cosmologicaldistances , The Astrophysical Journal (1989) 551–566.[54] D. J. Eisenstein and W. Hu, Baryonic features in the matter transfer function , TheAstrophysical Journal (1998), no. 2 605.[55] G. Efstathiou and J. R. Bond, Isocurvature cold dark matter fluctuations , Monthly Notices ofthe Royal Astronomical Society (1986), no. 1 103–121.[56] A. A. Klypin, S. Trujillo-Gomez, and J. Primack, Dark matter halos in the standardcosmological model: Results from the bolshoi simulation , The Astrophysical Journal (2011),no. 2 102.[57] M. S. Warren, K. Abazajian, D. E. Holz, and L. Teodoro, Precision determination of the massfunction of dark matter halos , The Astrophysical Journal (2006), no. 2 881.[58] P. Col´ın, A. Klypin, O. Valenzuela, and S. Gottl¨ober, Dwarf dark matter halos , TheAstrophysical Journal (2004), no. 1 50.[59] M. Sasaki, P. C. Clark, V. Springel, R. S. Klessen, and S. C. Glover, Statistical properties ofdark matter mini-haloes at z¿= 15 , arXiv preprint arXiv:1405.4018 (2014).[60] S. Ando, Can dark matter annihilation dominate the extragalactic gamma-ray background? , Phys. Rev. Lett. (2005) 171303, [ astro-ph/0503006 ].[61] S. Ando and E. Komatsu, Anisotropy of the cosmic gamma-ray background from dark matterannihilation , Phys. Rev. D (2006) 023521, [ astro-ph/0512217 ].[62] G. Haines and A. G. Jones, Logarithmic fourier transformation , Geophysical JournalInternational (1988), no. 1 171–178. – 32 – 63] V. V. Petrov, “Limit theorems of probability theory. 1995.”[64] Fermi-LAT Collaboration, M. Ackermann et al., Limits on Dark Matter Annihilation Signalsfrom the Fermi LAT 4-year Measurement of the Isotropic Gamma-Ray Background , arXiv:1501.05464 .[65] Fermi-LAT Collaboration, A. Abdo et al., The Spectrum of the Isotropic Diffuse Gamma-RayEmission Derived From First-Year Fermi Large Area Telescope Data , Phys. Rev. Lett. (2010) 101101, [ arXiv:1002.3603 ].[66] S. Ando and D. Nagai, Fermi-lat constraints on dark matter annihilation cross section fromobservations of the fornax cluster , Journal of Cosmology and Astroparticle Physics (2012), no. 07 017.[67] F. Zandanel and S. Ando, Constraints on diffuse gamma-ray emission from structure formationprocesses in the coma cluster , Monthly Notices of the Royal Astronomical Society (2014),no. 1 663–671, [ http://mnras.oxfordjournals.org/content/440/1/663.full.pdf+html ].[68] I. Arad, A. Dekel, and A. Klypin, Phase-space structure of dark matter haloes: scale-invariantprobability density function driven by substructure , Monthly Notices of the Royal AstronomicalSociety (2004), no. 1 15–29.[69] M. Ackermann et al., Constraining dark matter models from a combined analysis of milky waysatellites with the fermi large area telescope , Physical Review Letters (2011), no. 24.(2011), no. 24.