Limits on dark matter annihilation in prompt cusps from the isotropic gamma-ray background

Recent studies indicate that thermally produced dark matter will form highly concentrated, low-mass cusps in the early universe that often survive until the present. While these cusps contain a small fraction of the dark matter, their high density significantly increases the expected gamma-ray flux from dark matter annihilation, particularly in searches of large angular regions. We utilize 14 years of Fermi-LAT data to set strong constraints on dark matter annihilation through a detailed study of the isotropic gamma-ray background, excluding with 95% confidence dark matter annihilation to $b\bar{b}$ final states for dark matter masses below 120 GeV.


I. INTRODUCTION
While the particle properties of dark matter are not understood, its gravitational properties are strongly constrained by observations across a wide range of distance scales from the cosmic microwave background [1] to the smallest dwarf galaxies [2].These constraints require only that dark matter has a small scattering cross section and small thermal motions in the low-redshift universe.Under these conditions, both theoretical arguments and simulations imply that dark matter structures grow hierarchically from the "bottom up", producing quasi-equilibrium dark matter halos with density profiles that diverge approximately as r −1 in their inner region [3].
Recently, several studies have provided an important extension to theoretical expectations for the small-scale density distributions of gravitationally interacting dark matter.They have shown that, so long as the dark matter is initially smooth on the smallest scales, a pronounced cusp forms at the moment of collapse of each peak in the initial mass density field [4][5][6].Unlike the dark matter profiles produced subsequently by accretion and mergers, these "prompt cusps" have strongly peaked dark matter densities which rise as ρ ∝ r −1.5 outside a very small central core set by the initial phase-space density of the dark matter [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18].Contrary to previous expectations, these cusps can survive until the present era unless they are disrupted by tidal encounters with individual stars [19], a process that is only important for the small fraction of cusps that reside near the centers of galaxies [19,20].
While prompt cusps only provide O(1%) of the total dark matter content and are too low-mass and diffuse to produce measurable dynamical or gravitational lensing effects in cold dark matter scenarios, their extremely high densities greatly enhance the dark matter annihilation signal in models where such interactions are possible.Within the paradigm of weakly interacting massive particles (WIMPs) [22], prompt cusps 1. Upper limits on the WIMP annihilation cross-section, accounting for prompt cusps that form at early times and persist throughout the bulk of both the Milky Way and extragalactic halos.The dotted line shows the expected thermal WIMP annihilation cross-section [21].For dark matter with the thermal cross section, constraints from the isotropic γ-ray background obtained from our analysis exclude at 95% confidence masses below 120 GeV annihilating to b b and masses around 90 GeV annihilating to W + W − .Constraints lie near the thermal cross-section for τ + τ − .dominate the total dark matter annihilation signal from all but the densest regions [20], significantly changing the flux and morphology of the dark matter annihilation signal.
Most importantly for γ-ray searches, the existence of a large number of prompt cusps implies that on astrophysical scales the dark matter annihilation signal traces the distribution of cusp number density, rather than the square of the (smooth) dark matter halo density.This implies that the highest sensitivity is achieved by searches of large regions containing the bulk of the dark matter, rather than by searches focused on regions of high local dark matter density.Note that individual prompt cusps are unlikely to be resolved [23].
In this paper, we make a new estimate of the isotropic γ-ray background (IGRB) based on 14 years of Fermi-LAT data.Modeling the γ-ray contribution from blazar and non-blazar source classes, as well as the simultaneous contribution from extragalactic star-forming (SF) activity, we find no evidence for an excess and strongly constrain any contribution from dark matter, as seen in Fig. 1.Comparing our results to the flux predicted from prompt cusps in the Milky Way halo and extragalactic environments, we place upper limits on the dark matter annihilation cross-section, excluding with 95% confidence the thermal annihilation cross-section for dark matter particle masses below 120 GeV annihilating to b b.Notably, this limit is in tension with models where dark matter annihilation produces the Galactic Center γ-ray excess [24].

II. PROMPT CUSPS IN WIMP SCENARIOS
Prompt cusps have ρ = Ar −1.5 density profiles that extend out to a limiting radius r cusp .Each prompt cusp formed from the collapse of a local maximum in the initial linear density field, and Ref. [5] showed that the density normalization A and the radius r cusp for each cusp are directly connected to the properties of its progenitor density peak.Specifically, numerical simulations suggest that is accurate at the 10 percent level for individual prompt cusps.
Here ρ0 is the comoving dark matter density, a coll is the cosmic scale factor when the cusp's progenitor peak collapses (according to the ellipsoidal collapse approximation in Ref. [25]), and R ≡ (δ/|∇ 2 δ|) 1/2 is the peak's characteristic comoving size, defined in terms of its height δ ≡ (ρ − ρ)/ρ and curvature ∇ 2 δ.We adopt ρ0 = 33.1 M ⊙ kpc −3 [1].Furthermore, the relation is accurate to within a factor of a few for individual prompt cusps, with uncertainty dominated by the unclear operational definition of a prompt cusp's edge.Our results are only logarithmically sensitive to r cusp .As Ref. [5] noted, Eqs. ( 1) and ( 2) are demanded by dimensional considerations; only the numerical factors are tuned to match simulations.Since the annihilation rate within a ρ ∝ r −1.5 density profile diverges logarithmically at small radii, it is also necessary to quantify the minimum radius at which the profile is valid.As noted by Ref. [5], thermal motion in the initial conditions imposes a maximum phase-space density f max that cannot be exceeded by any gravitationally evolved structure.This leads to the requirement that the ρ ∝ r −1.5 profile transitions to a finite-density core within the core radius [5] where G is the gravitational constant.The density of a cored prompt cusp changes from ρ = Ar −1.5 for r ≫ r core to ρ ≃ Ar −1.5 core for r ≪ r core .Here we adopt the cored cusp description of Ref. [19], which is derived from an Ansatz for the phase-space density and for which the volume integral of the squared density (out to r = r cusp ) is The annihilation rate inside a prompt cusp is proportional to J. Note that typically r cusp /r core ∼ 500 for WIMP models [20], so the quantity in brackets in Eq. ( 4) is approximately 7.  Prompt cusps accrete dark matter halos around them, which also contribute to the annihilation signal.However, Ref. [20] noted that the annihilation rate in prompt cusps exceeds that predicted by previous halo and subhalo models (e.g.[26]) by at least an order of magnitude in the cosmological average, which suggests that halos only contribute to annihilation around the 10 percent level.We neglect their contribution.

A. The prompt cusp distribution
We evaluate the abundance of prompt cusps and the distribution of their structural properties as in Ref. [20].Due to Eqs.
(1) and ( 2), these follow directly from the properties of the peaks in the initial linear density field, and thus can be calculated explicitly from the linear power spectrum of dark matter density fluctuations [27].Adopting Planck 2018 cosmological parameters [1], we use the CLASS code [28] to evaluate the linear dark matter power spectrum at redshift z = 31.The power spectrum at k > 2 × 10 4 Mpc −1 is then evaluated analytically according to Ref. [29] and matched in amplitude to the CLASS output.The resulting power spectrum is shown as the black curve in Fig. 2 and is valid for dark matter models with no initial thermal motion.Such models give rise to density variations on arbitrarily small scale.
However, WIMP dark matter was once in equilibrium with the Standard Model plasma, and the resulting thermal motion smooths out density variations below some characteristic freestreaming scale k −1 fs .The smoothing scale in the initial dark matter distribution is of fundamental importance to the distribution of prompt cusps [30].We evaluate k fs via Eq.( 48) of Ref. [31] and multiply the idealized power spectrum in Fig. 2 by the Gaussian function exp(−k2 /k 2 fs ).The colored curves show the resulting power spectra for a sample of WIMP scenarios.Particles of lower mass, m DM , have higher thermal velocities, so they stream over greater distances, leading to lower values for k fs .The temperature T d at which the dark matter kinetically decoupled from the Standard Model plasma is also relevant; later decoupling (lower T d ) also leads to larger streaming distances and hence lower k fs .
For each dark matter power spectrum, we use the procedure described in Appendix C of Ref. [4] to evaluate the number density n peaks of initial density peaks and to generate a Monte Carlo sample of 10 7 such peaks.These peaks are distributed in δ, ∇ 2 δ, and the ellipticity e and prolateness p of the tidal field.We note that this procedure is exact, aside from the distribution of e and p, which is still very nearly exact for the power spectra that we consider [20].Each peak's collapse time follows as a coll = [δ ec (e, p)/δ]1/g a, where a = 1/32 is the scale factor at which the power spectrum was evaluated, δ ec (e, p) is the threshold for ellipsoidal collapse given by Ref. [25], 1 and g ≃ 0.901 [29] due to the baryon-driven suppression of dark matter density perturbations at small scales (an effect responsible for the feature in Fig. 2 around k ∼ 10 2.5 Mpc −1 ).The random sample of peaks is then translated into a sample of prompt cusps using Eqs.(1-3), with f max determined as in Ref. [20].

B. Annihilation in prompt cusps
We are interested the annihilation signal from a mass M of dark matter that is large enough to contain many prompt cusps.The volume integral of the squared density contributed by prompt cusps in this case, per M , is = n peaks ⟨J⟩ sample /ρ 0 (6) on average, where the sum is over the values of J (Eq. 4) for the cusps inside the mass M and the angle brackets denote the average over the Monte Carlo sample of cusps obtained earlier. 2 Note that (per Eq. 6) ρ eff does not depend on M .It also has negligible variance between regions as long as M is much larger than the mass scale of individual cusps.ρ eff is dimensionally a mass density; in this respect, it is the prompt cusps' contribution to the mass-weighted average dark matter density, since ρ eff = 1 M cusps ρ dM .It is a convenient quantity because the average annihilation rate Γ per dark matter mass M , due to prompt cusps, is where ⟨σv⟩ is the annihilation cross section.Here we assume the dark matter is its own antiparticle and that σv is velocityindependent at lowest order.We insert the subscript '0' to note that this relationship does not account for the further evolution of the cusps after their formation, which we will discuss below.
Figure 3 shows how ρ eff varies with the mass m DM and kinetic decoupling temperature T d of the dark matter particle.It is approximately logarithmically sensitive to both of these parameters. 3We will assume T d = 30 MeV for the remainder of this article, a common choice (e.g.Refs.[32,33]) that is relatively conservative.Models with lower T d couple more strongly to Standard Model particles, which can put them at odds with nondetection results in terrestrial experiments [34].Higher T d , on the other hand, would only increase ρ eff , boosting the annihilation signal further and strengthening any inferred constraints on the annihilation cross section.
Not all prompt cusps survive.Halos frequently merge together, and if they have comparable masses, it is expected that their central cusps also merge.This process effectively removes one prompt cusp while only slightly altering the other [5,35].Also, a substantial number of the lower initial density peaks do not even form prompt cusps because their mass accretes onto another object before they collapse.We model these effects by scaling ρ eff by a factor f glob < 1.We adopt f glob = 0.5 based on recent simulations [20], but this quantity is uncertain at the ∼ 0.1 (∼ 20 percent) level.Note that f glob propagates cleanly through our analysis, so that our constraints on ⟨σv⟩ can be scaled by 0.5/f glob if further research refines the predicted value of f glob .
In a merger between halos of highly unequal mass, the smaller object generally survives as a subhalo of the larger 3 Reference [20] found that ρ eff ≃ 0.08 log object.Today, almost all prompt cusps are expected to be subhalos of larger systems.However, even the smallest observationally resolved halos (e.g.those of the Milky Way's faintest satellites) are still expected to contain a very large number of prompt cusps (∼ 10 10 , perhaps).Subhalos are gradually stripped by tidal forces (e.g.[36][37][38][39]), suppressing their annihilation radiation [40,41].However, this only significantly affects prompt cusps in the innermost regions of present-day halos where the dominant disruption effect is, instead, encounters with individual stars (e.g.[19,42,43]).In a cosmological average, these effects are generally minor.For example, Ref. [20] used the tidal stripping model of Ref. [39] to estimate that tidal stripping suppresses the observable annihilation signal from external galactic halos at the 10 percent level.Similarly, Ref. [44] estimated that under 10 percent of dark matter resides in regions where stars are abundant, such that disruption by stellar encounters would be relevant.Since these effects are smaller than the uncertainty in f glob , we neglect them.However, if further research yields more precise modeling, the global effects of mergers, tidal stripping, and stellar encounters, as well as any small corrections to the prompt cusp density profiles (Eqs.1-4), can be absorbed at first order into f glob , scaling our constraints on ⟨σv⟩ by 0.5/f glob as discussed above.
Tidal stripping and encounters with stars are, however, highly relevant for the prompt cusp signal from our own Galactic halo as a result of our position near its center.In this context, we will denote by f gal (x) < 1 the scaling factor for the prompt cusp annihilation rate due to tidal stripping and stellar encounters, which is a function of position x.

III. THE ISOTROPIC ANNIHILATION SIGNAL
We now quantify the contribution of annihilating dark matter to the IGRB resulting from the prompt cusp distributions derived above.Let dN γ /dE be the average differential photon count, per energy E, produced by each annihilation event.We will consider b b, τ + τ − , and W + W − annihilation channels, adopting in each case the photon spectra dN γ /dE tabulated by Ref. [45] (which use the results of Ref. [46]).Note that the differential photon rate emitted per mass M of dark matter is then ⟨Γ/M ⟩dN γ /dE, where ⟨Γ/M ⟩ is the average annihilation rate per dark matter mass.We will neglect the contribution of dark matter not in prompt cusps.Even though only a percent-level fraction of the dark matter is in cusps, Ref. [20] used comparisons with previous halo-based predictions of the annihilation rate to show that the prompt cusps dominate the average signal.

A. Galactic contribution
A significant fraction of the dark matter contribution to the IGRB comes from our own Galactic halo (e.g.[47]).Our halo has a dark matter distribution that we approximate with a spherically symmetric density profile ρ gal (r).According to the above considerations, the average annihilation rate per dark matter mass at radius r due to prompt cusps is obtained by scaling Eq. ( 7) by the suppression factors discussed in Sec.II.As noted above, we only include the contribution of prompt cusps to the annihilation signal.While the smooth Galactic halo dominates the signal from the Galactic Center, it contributes only at the percent level beyond ∼20 o therefrom.
We use the CUSP-ENCOUNTERS code [19] to evaluate f gal (r)ρ eff .This code adopts the description of the Galaxy's baryonic components used by Ref. [48], which includes axisymmetric stellar and gas disks and a stellar bulge, the parameters of which are consistent with observational constraints in Refs.[49,50].The dark matter halo is taken to initially have a density profile of the Navarro-Frenk-White form [3] with concentration c ≡ R 200c /r s = 8.7 and mass M 200c = 10 12 M ⊙ , but it is subjected to adiabatic contraction due to the presence of the baryons according to the prescription of Ref. [51].We will test the impact of the choice of c and M 200c .For the initial sample of prompt cusps derived in Sec.II, orbits are randomly drawn from the distribution function of the Galactic halo (assuming velocity isotropy for simplicity).The prescriptions of Refs.[19,39] are then applied to modify the cusps' density profiles due to tidal stripping from the Galactic potential and impulsive encounters with stars along their orbits.
The photon flux that we receive from annihilation in cusps in the Galactic halo, per area and solid angle, is then where ρ(r) is the density profile of the Galactic halo and r = r 2 0 + ℓ 2 − 2rℓ cos θ is the Galactocentric radius at distance ℓ along the chosen line-of-sight, which is at angle θ from the Galactic Center.Here r 0 = 8.2 kpc is the distance to the Galactic Center.
Equation ( 9) is a function of both the sky angle θ and the photon energy E. To illustrate the angular dependence, we consider the integrated energy flux, dEEd 2 Φ/dΩdE.For the example model of a 100 GeV WIMP that decouples at a temperature of 30 MeV and annihilates into b b with cross section ⟨σv⟩ = 10 −26 cm 3 s −1 , we plot the integrated flux from Galactic cusps as the green curve in Fig. 4. In our analysis we will consider only Galactic latitudes |b| > 20 deg. and longitudes |l| > 80 deg., as we will discuss in Section IV.Within this region of the sky, the angle θ from the Galactic Center ranges from 80.6 deg. to 160 deg., and the annihilation signal from Galactic cusps differs from its average value (dotted line) by up to 30 percent.While the density of the Galactic halo rises strongly toward the Galactic Center, this trend is not reflected as strongly in the annihilation rate because prompt cusps nearer to the Galactic Center are suppressed more by stellar encounters and (to a lesser extent) by tidal forces, as discussed in Ref. [19].For dark matter models of higher mass, the prompt cusps form slightly earlier and thus have higher density, reducing their susceptibility to suppression by these effects.However, this effect is small; for a 100 TeV WIMP, the annihilation signal from Galactic cusps differs from its average value by a maximum of 38 percent.The parameters of the Galactic halo are not exactly known.Figure 5 shows how the sky-averaged flux from annihilation in Galactic cusps changes as the halo mass M 200c and initial NFW concentration c are varied across the uncertainty range of the observational inference in Ref. [51].More massive and more concentrated halos yield stronger signals.However, within the 68 percent confidence range, and even most of the 95 percent confidence range, the flux from Galactic cusps varies by under 10 percent from the fiducial value (which adopts c = 8.7 and M 200c = 10 12 M ⊙ ).We will neglect this uncertainty.As we will see next, Galactic cusps only source about a third of the total annihilation signal, so the total effect of the uncertainty in Galactic parameters is much smaller than 10 percent.

B. Extragalactic contribution
The average annihilation rate per dark matter mass in the extragalactic field is The photon flux from annihilation in extragalactic cusps is As we discussed in Sec.II B, tidal stripping and stellar encounters (incorporated via f gal in Eq. 9) are small effects in the cosmological average, so we do not explicitly include them here.Due to the cosmological redshift, extragalactic cusps contribute with a different energy spectrum from Galactic cusps.This difference is illustrated in Fig. 6 for the same example 100 GeV WIMP model.Extragalactic cusps (orange curve) contribute about twice the energy flux of Galactic cusps (green curve), and their contribution is skewed toward slightly lower energies.The integrated energy flux from extragalactic cusps is shown in Fig. 4 (orange line) for the same model.It is evi-dent again that extragalactic cusps contribute at about twice the level of Galactic cusps.We also show the total cusp signal (blue curve), summing Galactic and extragalactic contributions.Within the region of the sky under consideration, at Galactic latitudes |b| > 20 deg. and longitudes |l| > 80 deg., the total signal differs from its average value (dashed line) by well under 10 percent.Such a small variation would be indistinguishable from an isotropic signal, and thus for the remainder of the analysis, we assume that the dark matter signal is isotropic and consider only its average value.

IV. ISOLATING THE OBSERVED ISOTROPIC BACKGROUND FLUX
With the predicted dark matter contribution to the IGRB established, we now begin our analysis of the observed IGRB using 14 years of Pass 8 Fermi-LAT data.We first separate the diffuse isotropic emission from anisotropic foreground emission components using a morphological analysis.The γ-ray foreground includes contributions from Galactic and extragalactic point sources, as well as diffuse Galactic γ-ray emission.To separate these components, we use a model very similar to "Model A" in Ref. [52], with a few modifications.In particular, we use different angular cuts: |l| > 80 • and |b| > 20 • , which excludes slightly more of the Galactic disk and Galactic Center.The excluded regions are more strongly affected by foregrounds and hence provide lower signal-tonoise when analyzing the IGRB.
The diffuse Galactic emission comes from a combination of the inverse Compton scattering (ICS) of starlight by relativistic electrons, as well as the bremsstrahlung and hadronic interactions of relativistic electrons and protons, respectively, with Galactic gas.The morphology of this emission closely traces the Galactic gas density, which is traced in turn by a combination of several key gas tracers, including: • COuse the radio observations of carbon monoxide (CO) from the Planck survey [53] as a tracer of cold molecular hydrogen gas.The component of CO-traced H 2 has a small scale height, and thus, due to our cut on Galactic latitude, gives only a minuscule contribution.
• Dust -We use a three-dimensional dust map based on Gaia data [54] to trace the cold neutral medium (CNM).
• 21 cm -We use the two-dimensional full-sky 21 cm map called HI4PI [55], as a tracer of atomic hydrogen.
We refer to Ref. [52] for a more detailed account of these γray components.We infer the γ-ray source normalizations in a Bayesian framework, using Hamiltonian Monte Carlo, where each energy bin is fitted separately.We use strictly positive wide flat box priors for all γ-ray components.Furthermore, we perform a jackknife binning of the sky area to quantify systematic uncertainties.We divide the sky into ten jackknife regions, which all cover half of the non-masked sky.These regions are defined as follows.For each jackknife region and energy bin, we sample the posterior density using Markov chain Monte Carlo (MCMC).After a thorough burn-in phase, where the mode is located and the step-size is tuned, the MCMC is run for 80 000 steps.The variance of the full chain is roughly equal to the variance between realizations separated by 20 steps.We thin the chain by this factor, giving a total of 4000 MCMC realizations.The uncertainty given by the difference between jackknife regions is dominant and significantly larger than the statistical uncertainty associated with the individual MCMC chains.The different jackknife regions not only allow us to quantify systematic uncertainties, but they also enable us to estimate possible correlations between the respective energy bins.If the fitted IGRB in different energy bins has the same upward or downward fluctuations in the same jackknife regions, that implies that they are subject to a similar fitting error; in other words, there is a positive correlation between those energy bins.We calculate the fitted IGRB correlation matrix over the different energy bins, by taking the correlation factor of the respective jackknife region's MCMC realizations. 4The correlation matrix is visible in Fig. 7.Although this approach gives a handle on correlations, it is not complete; in principle, there could also be systematic correlations between energy bins that 4 In exact terms, the correlation factor between two energy bins, labeled by indices i and j, is calculated according to: Corr. factor (i,j) = where k denotes the jackknife region, n denotes the MCMC realization index, and x i,k,n is an MCMC realization of the IGRB normalization factor.The quantities xi and Std.(x i ) represent the mean and standard deviation of the MCMC chains of the ith energy bin, all jackknife regions included.The numerical factor in the denominator is equal to the product of the number of jackknife regions and number of MCMC iterations.are not dependent on sky angle, which we would not be sensitive to with this procedure.
The black points in the upper panel of Fig. 8 show the isotropic energy flux that we obtain by this procedure.We note that this is not exactly the IGRB, because it includes likely contamination by misidentified cosmic rays, which we will discuss and account for in Sec.VI.

V. ASTROPHYSICAL CONTRIBUTIONS TO THE IGRB
The majority of the IGRB is known to originate from astrophysical sources (e.g.[47,56]).In order to model these contributions, we consider spectral components from blazar and non-blazar sources.Our non-blazar component includes extragalactic star-forming activity (SF) and misaligned AGN (mAGN).Our blazar component includes BL Lacertae objects (BL Lacs) and flat spectrum radio quasars (FSRQs).We adopt the IGRB model in Ref. [56] and allow the relative contribution of each source class to vary according to the lognormal distributions found therein.This model accounts for the potential IGRB contribution from any galaxy that hosts a non-blazar supermassive black hole since it allows for misaligned AGN gamma-ray emission to be considered independently from the gamma-ray emission originating from its starforming activity.
This model also accounts for the electromagnetic cascades that originate from γ-rays above about 100 GeV as they generate high-energy electrons when scattering with the extragalactic background light (EBL+CMB) [57,58].In turn, these electrons produce gamma rays after inverse-Compton scattering with the same radiation fields.It is important to account for these effects since the spectral shape of any sufficiently energetic source class is significantly altered by these cascades, especially above 100 GeV.While the blazar model that we adopt does not account for cascade production, we note that FIG. 8.The isotropic emission (black points) and the results of our fit to it (black curve).This emission includes contributions from cosmic rays that have been misidentified as γ rays by the Fermi-LAT instrument, γ-ray contributions from extragalactic star-forming activity, blazar and non-blazar AGN source classes, as well as potential contributions from dark matter annihilation.We find that our combined model, including astrophysical γ-ray emission (blue and orange curves) and misidentified cosmic rays (green curve), fits the observed data without any evidence for an excess from dark matter annihilation.Here we show the 95 percent confidence limit on the emission from dark matter with a mass of 100 GeV annihilating into b quarks (dashed curve).
(1) the contribution from unresolved blazars cannot be significantly enhanced by their cascades [59] and (2) the blazar contribution is expected to be subdominant to SF and mAGN at all relevant gamma-ray energies [56].
Since there is a considerable degeneracy in the spectral shape of the SF, mAGN, and BL Lac contributions, their spectral shape is considered to be the same while the normalization is fit using their joint log-normal distribution as a prior.We note that the summed log-normal distribution can be directly calculated by using the full covariance matrix for the individual uncertainties calculated in Table II of Ref. [56].Because the spectrum of the FSRQ contribution differs from the remaining astrophysical components, the FSRQ source class is considered independently.As shown in Fig. 8, our model can account for the entirety of the measured isotropic gamma rays, allowing us to constrain the addition of a potential dark matter signal.

VI. DARK MATTER LIMITS
Dark matter annihilation and astrophysical sources would both contribute to the diffuse isotropic emission that we ob- tained in Sec.IV.To obtain a limit on the dark matter annihilation cross section, we now employ a combination of spectral information and multiwavelength constraints to differentiate the contributions from each component.In particular, we utilize a template fit approach that incorporates the astrophysical sources discussed in Sec.V, the contamination from cosmic ray (CR) backgrounds, and the Galactic and extragalactic dark matter signal outlined in Sec.III.
To account for uncertainties in the model of astrophysical sources and the CR background contamination, we introduce nuisance parameters.Subsequently, we perform a template fit to the measurement of the isotropic flux using four templates: SF+mAGN+BlLac, FSRQ, CR background contamination, and a dark matter model.The astrophysical contribution of SF, mAGN, and BlLac have similar spectral shapes and exhibit significant degeneracy in normalization, whereas the sum of their contributions is more tightly constrained.Hence, we opt to consider a single template that encompasses all three populations.In our default simulations, the CR background is set to the results of Ref. [60] and represents instrumental misidentification of charged CR events as a γ-ray flux.We note that Ref. [60] was calculated using Pass 7 data, and we discuss in detail systematic tests of uncertainties in the CR background below.
The free parameters involved in the fit are: • The normalization, A SF+mAGN+BlLac , of the astrophysical contribution from with respect to the model from Ref. [56].We choose a Gaussian prior in log 10 A SF+mAGN+BlLac centered at µ = 0 with width σ determined by the uncertainty from Ref. [56].
• The deviation of the spectral slope, δ SF+mAGN+BlLac , of the SF+mAGN+BlLac contribution with respect to the reference model.We apply a flat prior for between −0.2 and 0.2.
• The normalization, A FSRQ , of the astrophysical FSRQ contribution.The prior is again Gaussian in log 10 A FSRQ with µ = 0 and σ taken from Ref. [56].
We perform 62 fits, one without dark matter (i.e.⟨σv⟩ fixed to 0) and the remaining 61 fits using discrete values of the dark matter mass approximately logarithmically spaced between 5 GeV and 100 TeV.For each fit we minimize the negative loglikelihood defined as where C −1 is the matrix inverse of Here ϕ (d) is the measured isotropic flux with the covariance matrix V , ϕ (m) is the sum of the model templates, and σ (CR,bkg) is the uncertainty of the CR background contamination extracted from Ref. [60].
For each fit, we perform a comprehensive MultiNest parameter scan to obtain the full posterior distribution.First we analyze the fit without dark matter.By utilizing solely the astrophysical source and CR contamination model, we achieve a satisfactory fit to the data.Moreover, the posterior distributions of the nuisance parameters, introduced for the astrophysical source populations, align well with the reference model normalizations, and the deviation of the slope is consistent with zero.The marginalized posteriors for each parameter are visualized by the dashed black lines in Fig. 9.The χ 2 value amounts to 17.5 for 15 data points, allowing us to proceed with deriving dark matter limits.
Considering an example involving dark matter, the blue lines in Fig. 9 show the posteriors for m DM = 100 GeV and annihilation into b quarks.In this section, we utilize annihilation into b quarks as an example to investigate systematic uncertainties.However, we will present the limits for annihilation into τ leptons and W gauge bosons in the concluding section.At 100 GeV there is no preference for dark matter, as can be seen from multiple angles.First, the χ 2 of the fit does not improve.Second, the posterior distribution of the parameters describing the astrophysical templates are not affected by the additional dark matter template.Third, the posterior distribution in ⟨σv⟩ is flat at low cross sections.We also show the posterior distribution in more detail in Fig. 10.Notice in particular that there is no degeneracy between the dark matter cross section and astrophysical parameters.
FIG. 10. Posterior distribution of the three parameters of the astrophysical model and the dark matter annihilation cross section, for a dark matter mass of 100 GeV.The inner and outer contours enclose 68 and 95 percent of the distribution, respectively.
The dark matter limit is determined from the difference of the marginalized χ 2 between the cases with and without dark matter annihilation, which is given by Here L(⟨σv⟩) is the likelihood defined in Eq. ( 13) that is marginalized over the vector of astrophysical parameters (denoted θ), namely: where π(θ) is the prior distribution in θ.Finally, the dark matter limit at 95% confidence level is obtained at ∆ χ2 (⟨σv⟩) = 3.84. 5Technically, L(⟨σv⟩) is obtained from the MultiNest scans with dark matter annihilation, while L(⟨σv⟩ = 0) is taken from the scan without it. 6The rightmost panel in Fig. 9 shows the marginalized χ 2 and the dark matter limit for the example dark matter masses of 100 and 300 GeV.In Fig. 8 we compare our model with the isotropic flux data (black points).We show fits for the astrophysical templates that correspond to the median parameters at the peak of the marginalized posteriors, while the dark matter template is shown for ⟨σv⟩ at the 95% limit.The black line shows the sum of the astrophysical templates (i.e.without dark matter).The residuals, defined as (data-model)/model, in the lower panel are also evaluated for a model that does not include the dark matter template.
The cross-section limit for dark matter masses ranging from 5 GeV to 5 TeV is represented by the solid blue line in Fig. 11 As expected, stronger limits are found for smaller dark matter masses.The precise increase of the upper limit results from the interplay between the decreasing number density of dark matter particles (scaling with m −1 DM ) and the decrease of the astrophysical background flux at higher energies.At 10 GeV, we obtain a limit of about 6×10 −27 cm 3 /s for dark matter annihilation to a b b final state, which then rises to approximately 4 × 10 −25 cm 3 /s at 1 TeV.
We note that around 500 GeV the dark matter limit appears relatively weak.This is due to a slight preference for the dark matter template.However, this preference is not statistically significant.To assess this, we conduct an additional fit where the dark matter mass is treated as a free parameter.The resulting χ 2 difference between the fit with and without dark matter is 9.8, corresponding to a local significance of less than 3σ for two additional fit parameters employing the frequentist interpretation.Furthermore, the Bayes factor, which is the ratio of the fit with and without dark matter, is 2.2, leads to the same conclusion.The best χ 2 and evidences of all fits are summarized in Tab.I.The significance of the dark matter preference is not robust against the systematic uncertainties (see discussion below) and, therefore, should not be taken seriously.
While our default approach involves incorporating correlated uncertainties for the isotropic flux measurement, we also evaluate our constraints in models where we do not include correlations between different energy bins.This approach is TABLE I. Results of the MultiNest scans, comparing fits with and without the inclusion of dark matter.The fits with dark matter (DM) are based on two additional free parameters: the dark matter mass and the annihilation cross section.The provided values include the best-fit χ 2 for each case, along with the corresponding log-evidences.We also provide the dark matter mass and the annihilation cross section for the best fit.While the inclusion of dark matter shows slight improvements in the fits, these enhancements are not statistically significant (see text for details).more conservative since it allows for greater flexibility in both the astrophysical model and dark matter, as the shape is less constrained by the correlation.This is evident from the minimal χ 2 values, which decrease to 4.6 in the fit without dark matter and 1.9 in the fit with dark matter.Consequently, the χ 2 per degree of freedom is significantly lower than one, providing further evidence that our estimation of the covariance matrix is reasonable.The dark matter limit obtained for this case is comparable to the default scenario.Only for M DM ranging from 20 to 80 GeV does this approach provides a slightly weaker limit by about a factor of 2 as shown in Fig. 11.

Fit
Another systematic uncertainty concerns the CR background contamination taken from Ref. [60].As detailed above, in the benchmark analysis we assume that the CR contamination rate in our analysis is the same.While this is a reasonable assumption, we note that there are some differences between our analysis and Ref. [60] that might impact the CR contamination.Ref. [60] uses specific additional cuts to the Pass 7 data with ULTRACLEAN veto in their analysis to decrease the CR contamination systematically, separately at low and high energies.We employ the updated Pass 8 data, which are supposed to feature a better background rejections but on the other hand we do not use additional cuts to suppress the contamination.
To address our incomplete knowledge of the CR contamination we performed two additional tests.First, we fixed the CR background template to 50% of our benchmark and, second, we allowed the normalization of the CR background (but not its spectrum) to vary as a free parameter in the fits.When the CR background template is fixed at 50%, we do not achieve a satisfactory fit to the isotropic flux using only the astrophysical model.The best-fit χ 2 value increases to 43, indicating a poor fit.Furthermore, the nuisance parameter δ SF+mAGN+BlLac converges to the lower boundary at −0.2.Consequently, we obtain a stronger dark matter limit, particularly at a few tens of GeV.However, this limit should be treated with caution.Conversely, when the normalization of the CR contamination template is allowed to vary as a free parameter, it converges to approximately 175%.In this scenario, the limit becomes more conservative, with the most significant impact observed once again at a few tens of GeV.We note that while neither of these techniques are optimal, the Fermi-LAT data on cosmic-ray contamination for the Pass 8 dataset has not been publicly released, making it difficult to self-consistently analyze the rate of cosmic-ray contamination.The approach shown here is effective in bracketing the possible levels of cosmic-ray contamination.
The summary presented in Table I shows that the preference for a dark matter signature is not robust against different systematic assumptions.While the local significance is measured at 2.6σ for our benchmark case, it diminishes to 1.9σ when the cosmic ray (CR) background contamination is allowed to vary in the fit, and further drops to 1.1σ in the case without correlations.Only in the scenario with CR background contamination fixed at 50% the significance rises to approximately 4.4σ.However, this particular case fails to provide a satisfactory fit to the overall data.This discussion highlights the critical importance of accurately estimating the CR background contamination for future dark matter searches.The importance of the CR background can also be seen from Fig. 8: The bump-like spectral shape of the CR background resembles that of a dark matter signature leading to a degeneracy with a potential dark matter signal.
As a final cross-check, we directly use the official IGRB measurement conducted by Fermi [60] in 2014 to derive the dark matter limit.In this particular case, the CR background contamination is already subtracted by the collaboration, thus eliminating its influence.Additionally, it is important to note that we cannot consider correlations in this case, since they are not provided by the collaboration.The limit obtained using this dataset is comparable, but it exhibits fewer discernible features and a smaller preference for the dark matter template.This can be attributed to the larger uncertainties associated with the IGRB measurement itself.

VII. DISCUSSION AND CONCLUSIONS
We have examined the consequences of highly concentrated prompt cusps of dark matter that form at the onset of nonlinear structure formation and survive throughout the bulk of present-day dark matter halos.Due to the contribution of these systems, we find that strong WIMP annihilation constraints come from the isotropic γ-ray background.These constraints can be stronger than those from targeted studies of regions with high central density, such as dwarf spheroidal galaxies.12.Comparison with other results.This figure plots the upper limit on ⟨σv⟩ as a function of the dark matter particle mass mDM for the case of annihilation into b quarks.The solid curve shows our constraint using prompt cusps, while the dashed curve shows the Fermi Collaboration's constraint using a halo model [33]; both of these constraints are based on measurements of the IGRB.The dot-dashed curve shows the Fermi Collaboration's limit from dwarf spheroidal galaxies [61].The shaded region is compatible with an annihilation origin for the Galactic Center γ-ray excess [62].
Utilizing 14 years of Fermi-LAT data, we produced a stateof-the-art model for the isotropic γ-ray background.We fit this model with a combination of astrophysical components produced by star-forming galaxies, misaligned AGN, blazars, and residual cosmic-ray contamination.We found no evidence for a remaining excess that might be due to dark matter annihilation, and so we are able to set strong limits on the dark matter annihilation cross-section.These limits are shown in Fig. 1 for dark matter annihilating into b quarks, τ leptons, or W bosons.
In Figure 12, we compare the results of our analysis with the Fermi collaboration's bound on the dark matter annihilation cross section [33] using the IGRB measurement from Pass 7 Fermi-LAT data [60].This result used substructure models that did not account for prompt cusps.Our IGRBbased limits are stronger by more than an order of magnitude throughout most of the mass range, allowing us to exclude with 95% confidence thermally annihilating WIMP models with particle masses below 120 GeV (for a b b final state).As noted by Ref. [20], the annihilation rate in prompt cusps exceeds the annihilation rate assumed by Ref. [33] by about a factor of 5. Thus, prompt cusps are responsible for the majority of the improvement, while the additional years of data and the IGRB astrophysical model that we adopt (Sec.V) are likely responsible for the rest.
We additionally compare our results to Fermi-LAT limits on annihilation in dwarf spheroidal galaxies from Ref. [61] (which differ only marginally from a recent update [63]) as well as dark matter fits to the Galactic Center γ-ray excess [62].The existence of prompt cusps is not expected to significantly change these estimates, because these searches target regions that are of such high average density that the smooth dark matter component should dominate the annihi-lation rate.We find that our IGRB constraints are stronger than the dwarf constraints throughout most of the parameter space, indicating that if present treatments of prompt cusps are accurate, then Fermi-LAT observations of the IGRB provide among the strongest indirect constraints on velocityindependent annihilation of O(100) GeV dark matter. 7Additionally, our constraints are in tension with the cross sections necessary for dark matter to produce the Galactic Center excess.
Finally, we note that the annihilation limits that we set in this work depend strongly on recent advances in our understanding of the astrophysical composition of the IGRB.Future γ ray and multiwavelength studies capable of decreasing the systematic uncertainty in our understanding of extragalactic star formation activity and mAGN contributions to the IGRB thus offer the potential to unlock extremely sensitive searches for dark matter annihilation with current and future γ-ray observations.

FIG. 2 .
FIG. 2. Power spectra of dark matter density variations, evaluated in linear theory at z = 31 and shown in the dimensionless form, P(k) ≡ [k 3 /(2π 2 )]P (k).The black curve shows the prediction for idealized cold dark matter.Density variations in WIMP models are smoothed by thermal motion in the dark matter, imposing a smallscale cutoff in P(k) that depends on the particle mass mDM and the temperature T d at which it decoupled from the Standard Model plasma.The colored curves show P(k) for such scenarios; different colors correspond to different mDM while different line styles correspond to different values of T d .

FIG. 3 .
FIG.3.The effective dark matter density ρ eff contributed by prompt cusps, relevant to the dark matter annihilation rate (see Eq. 7), as a function of the particle mass mDM.ρ eff is the annihilation J factor (proportional to the annihilation rate) per mass of dark matter, or equivalently, the mass-weighted average dark matter density.Different colors correspond to different kinetic decoupling temperatures.

FIG. 4 .
FIG. 4. Integrated energy flux from annihilation in prompt cusps, as a function of angle from the Galactic Center.We adopt here a 100 GeV WIMP that decouples at 30 MeV and annihilates into b b with a cross section of ⟨σv⟩ = 10 −26 cm 3 s −1 .

FIG. 5 .FIG. 6 .
FIG.5.Sensitivity of the annihilation signal from Galactic cusps to the mass M200c and initial NFW concentration c of the Milky Way halo.We consider a 100 GeV WIMP that decouples at 30 MeV.The color scale (with contour lines every 10 percent) indicates the fractional change in the sky-averaged gamma-ray flux from Galactic cusps, compared to the fiducial parameters c = 8.7 and M200c = 10 12 M⊙.For comparison, the dashed and dotted curves mark the 68 and 95 percent confidence ranges, respectively, of the observational inference in Ref.[51].Although the annihilation signal is stronger for higher M200c and higher c, its variation remains at the 10 percent level across the range of likely parameters.

FIG. 7 .
FIG. 7. Correlation coefficients for the 15 energy bins of the isotropic component.The diagonal is valued unity by definition.Otherwise nearby energy bins typically have strong positive correlations.

FIG. 9 .
FIG.9.The three panels on the left display the marginalized posterior distribution for the fit parameters of the astrophysical model.Dashed black lines show the fit without dark matter while the solid blue lines show the result of the fit for a fixed dark matter mass of 100 GeV and annihilation into b quarks.The rightmost panel displays the marginalized χ 2 as a function of the dark matter annihilation cross section.The blue (yellow) line corresponds to a dark matter mass of 100 (300) GeV.The vertical dotted lines indicate the 95% upper limit on the dark matter annihilation cross section.

FIG. 11 .
FIG.11.Upper limits on dark matter annihilation into b quarks.The solid curve is our fiducial result, while the other blue curves test the sensitivity of the result to different assumptions (see the text for more detail).The dotted black line shows the annihilation cross section expected for a thermal WIMP[21].
FIG.12.Comparison with other results.This figure plots the upper limit on ⟨σv⟩ as a function of the dark matter particle mass mDM for the case of annihilation into b quarks.The solid curve shows our constraint using prompt cusps, while the dashed curve shows the Fermi Collaboration's constraint using a halo model[33]; both of these constraints are based on measurements of the IGRB.The dot-dashed curve shows the Fermi Collaboration's limit from dwarf spheroidal galaxies[61].The shaded region is compatible with an annihilation origin for the Galactic Center γ-ray excess[62].