Primordial black holes as silver bullets for new physics at the weak scale

Observational constraints on gamma rays produced by the annihilation of weakly interacting massive particles around primordial black holes (PBHs) imply that these two classes of dark matter candidates cannot coexist. We show here that the successful detection of one or more PBHs by radio searches (with the Square Kilometer Array) and gravitational wave searches (with LIGO/Virgo and the upcoming Einstein Telescope) would set extraordinarily stringent constraints on virtually all weak-scale extensions of the standard model with stable relics, including those predicting a WIMP abundance much smaller than that of dark matter. Upcoming PBH searches have in particular the potential to rule out almost the entire parameter space of popular theories such as the minimal supersymmetric standard model and scalar singlet dark matter.


I. INTRODUCTION
The formation and growth of black holes (BHs) inevitably modifies the dark matter (DM) distribution around them.If DM is in the form of weakly interacting massive particles (WIMPs) which self-annihilate, the increase in DM density can significantly boost the annihilation rate.This process has been discussed in the context of supermassive BHs at the centers of galaxies [1][2][3][4][5][6][7] and intermediate-mass BHs [8][9][10].The argument has been more recently extended to the case of primordial black holes (PBHs), which can form before big bang nucleosynthesis [11,12] and could constitute a significant, yet subdominant, fraction of DM.In this case, the WIMP annihilation rate around PBHs would lead to a gamma-ray background exceeding the one observed by the Fermi Large Area Telescope (LAT), leading to stringent constraints on the relative PBH abundance, f PBH ¼ Ω PBH =Ω DM [13][14][15][16].
In this paper, we explore the compatibility of PBHs and WIMP DM from the opposite viewpoint.We focus on the prospects for discovering PBHs with upcoming radio and gravitational wave (GW) searches, and on the implications such a discovery would have on even a small relic density of WIMPs in the Universe.Specifically, we consider three discovery scenarios: (i) The detection of GWs produced by the merger of BHs with mass M ≲ 1 M ⊙ with LIGO/Virgo.(ii) The detection of GWs produced by the merger of Oð10 M ⊙ Þ BHs at redshift z > 40 with the Einstein Telescope.(iii) The detection of the radio emission produced by the accretion of gas onto 1-1000 M ⊙ BHs with the planned Square Kilometer Array (SKA).The scenarios we consider are summarized in Table I.
We estimate the abundance of PBHs in each scenario, given a number of detections (Fig. 1) and, from that, we calculate the gamma-ray luminosity of WIMP overdensities around PBHs in the Universe.By comparing this with the observed diffuse extragalactic gamma-ray flux and with unidentified gamma-ray point sources in the 3FGL TABLE I. PBH detection scenarios.For each of the three detection scenarios we consider, we show the benchmark PBH mass M ⊙ which we assume, as well as N min and N max , the smallest and largest number of detections we consider.Note that N max is the largest number of detections which would be compatible with current PBH constraints.

Sub-solar (O3)
High redshift (ET) SKA M PBH [M ⊙ ] 0.5 10 100 N min 1 1 1 0 N max 80 24000 80 Fermi-LAT catalogue, we show that a positive detection of even a small number of PBHs in any of the above scenarios would set extraordinarily stringent constraints on weakscale extensions of the standard model with stable relics, including those predicting a WIMP abundance much smaller than that of the total dark matter.
Our code is available on GitHub [19], and a link below each figure [20] provides the Python code with which it was generated.An archived version is available on Zenodo [21].

II. DISCOVERING PBHs WITH GRAVITATIONAL WAVES
Upcoming GW observations have the potential to unequivocally detect PBHs.We consider two main detection scenarios based on this channel.The first case is the observation of nearby mergers of sub-solar-mass BHs.The lightest known BH is in the x-ray binary XTE J1650-500, with a mass of 5.1 M ⊙ [22] or perhaps larger [23].In standard theories of stellar evolution, compact objects do not form below the Chandrasekhar mass of ∼1.4 M ⊙ [24][25][26].As a result, the detection of sub-solar-mass BHs would strongly suggest a primordial origin (though see Refs.[27,28] for alternative exotic formation mechanisms).The second case is the observation of BH-BH mergers at high redshift with proposed detectors such as the Einstein Telescope (ET) [29] and Cosmic Explorer (CE) [30].
At low redshift, mergers of astrophysical BHs are commonplace, but at earlier times BHs have typically not had sufficient time to collapse and form binaries [31].Above around z ∼ 40, the merger rate of astrophysical BHs should be negligible [32].High-redshift mergers therefore point towards a primordial BH population which formed much earlier [33].

III. PBH MERGER RATES: SETTING THE STAGE
To calculate the PBH merger rate as a function of the PBH abundance, we follow Refs.[34,35] (though see also Refs.[36][37][38][39][40][41][42]).In this scenario, pairs of PBHs decouple from the Hubble flow and form binaries before matterradiation equality [43][44][45].For a given value of (f PBH , M PBH ), we calculate the comoving PBH merger rate (in the source frame) as a function of redshift as Here, n PBH ¼ f PBH Ω DM ρ c;0 =M PBH is the comoving number density of PBHs in terms of the critical density today, ρ c;0 .The merger time distribution is Pðt merge Þ, which we evaluate at t merge ¼ t½z, the age of the Universe at redshift z.The merger time depends on the semimajor axis a and eccentricity e of the binary orbit [46].The distribution Pðt merge Þ is thus obtained from the distribution of Pða; ejM PBH Þ, calculated assuming a uniform initial distribution of PBHs and accounting for torques due to all other PBHs and smooth density perturbations [34].We also account for the impact of local DM halos on the binary properties [35], for which we assume a DM density profile with slope γ ¼ 9=4 [16], leading to a roughly 30% enhancement in the merger rate.
To obtain projected constraints on R (or equivalently f PBH ) from future GW observation campaigns, we follow LIGO/Virgo [47,48] and perform a Bayesian analysis. 1We assume that the number of observed mergers is Poisson distributed, PðNjΛÞ ¼ PoissðΛÞ, with the mean dz fðzÞ is the differential timevolume sensitivity for an observation time T [47,48].We calculate the comoving volume as a function of redshift dV c =dz [49], assuming the Planck 2018 best-fit cosmological parameters [50].The selection function fðzÞ encodes the detection probability for a given source at redshift z.
FIG. 1. Constraints on the PBH fraction from future observations of PBHs.Median value (solid line) and symmetric 90% credible intervals [17] (shaded region) of f PBH ¼ Ω PBH =Ω DM , assuming the observation of N PBH candidates.In blue, we assume the observation of BH mergers with component masses of 0.5 M ⊙ during LIGO O3.In orange, we assume BH mergers with component masses of 10 M ⊙ are observed at redshift z ≥ 40 during one year of operation of the Einstein Telescope.In green, we assume the observation of 100 M ⊙ PBHs in the Milky Way in radio and x-ray searches.The gray hatched boundaries show the current 95% upper limit on f PBH for each PBH mass [18]. 1 We will assume that the threshold for detection of mergers is high enough that there is no contamination from terrestrial "events" (i.e., p astro ≃ 1 for all observations).
As in Refs.[47,48], we assume a Jeffreys prior2 on Λ (though see Sec. 3 of the Appendix for a more detailed discussion of the prior dependence of our results).Using Bayes' theorem [52], the resulting posterior on Λ, given N observations, is The corresponding posterior on the redshift-averaged merger rate R ¼ Λ=hVTi is then We also incorporate a log-normal prior on the total timevolume sensitivity hVTi, with a 30% relative uncertainty in order to account for possible calibration uncertainties [48].
We marginalize over hVTi to obtain the marginal posterior for R, which can be transformed into a posterior on f PBH : We now apply this formalism to the specific cases of subsolar-mass PBHs observed with LIGO and high-redshift mergers observed with ET.

IV. SUB-SOLAR-MASS MERGERS
Let us assume that during a one-year campaign, LIGO's third observing run (O3) observes N events which are consistent with merging sub-solar-mass BHs.Let us also assume that the PBH mass function is monochromatic, and that the component masses of the merging BHs are recovered exactly.In this case, we choose as a benchmark a PBH mass of M PBH ¼ 0.5 M ⊙ .
LIGO O3 is sensitive to such mergers only in a small range of redshifts, z ≲ 0.02.We therefore assume that the PBH merger rate is constant over this range, replacing it with the redshift-averaged value R. We estimate the sensitive timevolume using the approximation of Ref. [53]: Here, T is the analyzable live time of the detectors, and D avg is the average detector range.This is obtained from D max , the horizon distance of the detector, by averaging the detector response over both orientation and location of the binary, which gives D avg ≈ D max =2.26.We assume a duration of one year for O3 with a 70% duty cycle [54].We also assume a further reduction of 5% due to data quality cuts [55].The LIGO O2 horizon distance for observing 0.5 M ⊙ mergers (with a signal-to-noise ratio of at least 8 in each of the two detectors) is ∼60 Mpc [56].From Ref. [54], we expect an improvement in the horizon distance of about 50% from O2 to O3.We therefore assume a horizon distance D max ¼ 90 Mpc (z ≈ 0.0192), giving a sensitive time-volume of hVTi ≈ 1.8 × 10 −4 Gpc 3 yr.The resulting posteriors on f PBH are shown in blue in Fig. 1.The gray shaded region shows the current constraint on 0.5 M ⊙ PBHs, coming from the observed merger rate during LIGO O2 [56].We note that the median value of f PBH scales as ffiffiffiffi N p .This is because the merger rate scales predominantly as f 2 PBH , reflecting the number of PBH binaries which are formed in the early Universe.

V. HIGH-REDSHIFT MERGERS
In the case of high-redshift mergers, we consider the sensitivity of the proposed Einstein Telescope (ET) [29] with an observing time of one year.We assume that N BH-BH merger events are observed at redshifts greater than z > 40.We again assume a monochromatic mass function, now at M PBH ¼ 10 M ⊙ , close to the maximum sensitivity of ET.We estimate the mean number of observed merger events by integrating Eq. ( 2) over z > 40, incorporating the full redshift-dependence of the merger rate RðzÞ.The differential time-volume sensitivity dhVTi=dz is calculated as in Eq. ( 2) with the selection function fðzÞ estimated from Ref. [57,58]: this corresponds to ∼5% of mergers detectable at z ¼ 40, decreasing to zero close to z ¼ 100.
The resulting posteriors on f PBH are shown in orange in Fig. 1, where we again cut off the posterior at the current upper limit (gray) from Ref. [35].

VI. DISCOVERING GALACTIC PBHs WITH SKA
Isolated black holes can in principle be detected by observing the broadband spectrum of radiation-spanning from radio waves all the way up to the x-ray bandproduced by the interstellar gas accreting onto them [59].This strategy can in particular be applied to the search for PBHs in the inner Galaxy and results in a bound on their abundance, as shown in Ref. [60,61].A more refined analysis based on state-of-the-art numerical simulation of gas accretion onto BHs recently allowed an even stronger bound to be placed on the abundance of a population of PBHs in the 1-100 M ⊙ mass window [62].
During the next decade, the Square Kilometer Array (SKA) observatory will allow us to perform radio observations with unprecedented sensitivity.In particular, it will become possible to search for a subdominant population of PBHs that amount to a small fraction (even less than a percent) of the dark matter in the Universe [63].
A complex data analysis will be needed to identify a possible population of black holes of primordial origin among the plethora of radio sources that SKA will detect.This will be based on a comparison with other wavelengths and on further considerations on the spatial distribution, mass functions, and other quantities that may allow us to disentangle a population of conventional astrophysical BHs from PBHs [60][61][62][63].Here we assume that, as a result of such data analysis, a certain number of radio sources observed by SKA are clearly associated with a population of PBHs, and we apply a procedure similar to the one described above to infer the cosmic abundance of PBHs from the number of such objects identified in SKA data.We focus in particular on the mid-frequency instrument (SKA1-Mid), which will observe the sky in the 350 MHz-14 GHz domain.We consider 100 M ⊙ PBHs because these massive objects would be more easily detected, given the ∝ M 2 scaling of the accretion rate, and more easily distinguished from a population of astrophysical BHs.The region of interest is the Galactic ridge (−0.9°< l < 0.7°; −0.3°< b < 0.3°), where a large concentration of molecular gas is present.
As a first step, we determine the distribution PðN SKA jf PBH Þ of the number of PBHs detectable by SKA1-MID (1.4 GHz), adopting the Monte Carlo simulations described in Refs.[60,62], and assuming gain G ¼ 15 K=Jy, receiver temperature T rx ¼ 25 K, sky temperature towards the GC T sky ¼ 70 K, bandwidth Δν ¼ 770 MHz, and ≃1000 h exposure, which corresponds to 85 nJy sensitivity.
Then, we apply the Bayesian inversion discussed above to compute the posterior distribution function Pðf PBH jN SKA Þ, where we again assume a Jeffreys prior on the mean number of expected observations.We show in green in Fig. 1 the resulting posteriors on f PBH .The current upper bound for 100 M ⊙ PBHs comes from the radio and x-ray searches mentioned above [62].

VII. CONSTRAINTS ON DM SELF-ANNIHILATION FROM GAMMA RAYS
Primordial black holes seed the growth of very steep dark matter halos that form via approximately radial infall.The theory of secondary infall [64] and simple 1D simulations [65] predict that these ultracompact minihalos (UCMHs) have ρðrÞ ∝ r −9=4 density profiles, which has been confirmed by recent 3D simulations [16].Since f PBH is at or well below the percent level in all but one of our detection scenarios, we can assume that UCMHs form in isolation, so we neglect the effects of PBH-PBH interactions on the UCMH profile.
Due to the steepness of the profile, the WIMP density reaches a maximum value at the "annihilation plateau," where the WIMP annihilation rate becomes equal to the Hubble rate.Due to the large resulting gamma-ray luminosities, UCMHs in the Milky Way would appear as bright point sources with no counterparts in other wavelengths.Previous analyses searching the 3FGL for WIMP subhalos [66][67][68] have identified 19 bright, high-latitude, nonvariable unassociated point sources that are spectrally compatible with annihilating WIMPs.As described in detail in the Appendix, we perform a Monte Carlo simulation to assess the observability of UCMHs by Fermi.We then use this to determine the 95% confidence level (C.L.) upper bound on the WIMP annihilation cross section in the zero-velocity limit ðσv rel Þ 0 .This upper limit depends on the PBHs' spatial distribution which we assume tracks the Milky Way DM distribution.We fix f PBH to the fifth percentile of the posterior Pðf PBH jNÞ, derived in the previous sections for the detection of N PBH candidates.We conservatively assume that all 19 compatible unassociated point sources are UCMHs and set the upper limit on ðσv rel Þ 0 by comparing with the expected number of UCMHs passing cuts on their integrated gamma-ray flux and Galactic latitude (given M PBH , m χ , and N).
Annihilation in UCMHs outside the Milky Way over all redshifts contributes to the diffuse, isotropic extragalactic background (EGB) [69][70][71], which has been measured by Fermi [72].This provides an additional very robust constraint on the DM self-annihilation cross section, since it requires no assumptions about the PBH spatial distribution.To set a conservative bound, we do not assume a particular background model.Instead, we compute the expected gamma-ray flux from UCMHs in each of Fermi's energy bins, and calculate the likelihood of such an excess above the observed flux using the statistical and systematic uncertainties.As for the point-source constraints, we fix f PBH to the fifth percentile for a given detection scenario.
An important difference with regard to standard indirect detection analyses is the scaling of signals with the fractional WIMP abundance f χ ¼ Ω χ =Ω DM for underabundant thermal relics.Typically, the DM annihilation rate depends on the combination f χ 2 ðσv rel Þ 0 , since it factors into terms dependent on the integrated DM density profile squared (J-factor) and the self-annihilation cross section.In the PBH scenario, the DM density profile itself depends on ðσv rel Þ 0 , since this sets the radius of the annihilation plateau.As a result, the DM annihilation rate (and thus the extragalactic diffuse flux from PBHs and expected number of unassociated point sources) depends on the combination f χ 4 ðσv rel Þ 0 ; this is derived in the Appendix.

VIII. RESULTS AND DISCUSSION
For each detection scenario in Table I, we show as a function of WIMP mass the 95% C.L. upper limit on f 4 χ ðσv rel Þ 0 in Fig. 2, where f χ ¼ Ω χ =Ω DM is the fractional contribution of a particle species to the cosmic DM density.This allows us to compare our projections with the theoretical predictions in cases where new particles constitute only a subdominant component of DM.The colored curves show the most stringent constraint arising from gamma-ray observations at a given WIMP mass, assuming annihilation into bb.For our projected limits assuming a small number of PBH detections (solid lines), point-source constraints dominate at low WIMP mass, while diffuse constraints are more relevant at high mass.This can be seen as a "kink" in each of the solid lines, above which diffuse constraints dominate.For larger numbers of PBH detections (dashed lines), diffuse constraints generally dominate (see the Appendix for a more detailed comparison of the limits).We find that a detection of Oð10Þ PBHs with any of the methods described above would rule out large ranges of standard-model extensions with stable relics at the electroweak scale.To illustrate this, we show in dark gray the envelope of the 95% C.L. profile-likelihood contours for the MSSM7 [73] and various GUT-scale SUSY models [74] obtained by the GAMBIT Collaboration.These contours show the zero-velocity limit of the annihilation cross section ðσv rel Þ 0 and so can be directly compared with our projected gamma-ray constraints assuming s-wave annihilation only.In light gray, we also show contours for scalar singlet scenarios [75].Requiring the WIMP to comprise at least 10% of dark matter, there are no models in the right part of the scalar singlet region below the line f χ 4 ðσv rel Þ 0 ¼ 10 −29 cm 3 =s, and in the leftmost sliver of the singlet contour and throughout the dark gray envelope all such models have f χ 4 ðσv rel Þ 0 > 10 −32 cm 3 =s.Importantly, since f χ scales approximately as hσvi −1 , the unitarity limit on the freeze-out annihilation cross section hσvi [76] yields, in the case of s-wave annihilation, a lower bound on the rescaled annihilation cross section today, This is shown as the lower black dotted curve in Fig. 2. In some scenarios, we are able to exclude below this unitarity bound.The unitarity bound can in principle be made lower in the case of p-wave annihilation, resonances, coannihilation, and kinematic thresholds, but for the models shown this is always counterbalanced by values of f χ that far exceed the unitarity lower bound.Note that, although the constraints we have presented do not apply to, e.g., super-WIMP DM candidates like the gravitino [77,78], they do affect asymmetric DM models [79,80], where the symmetric component is reduced via freeze-out.The detection of PBHs therefore places strong constraints on the identity of the particles which make up the dominant fraction of DM.

IX. CONCLUSIONS
Upcoming gravitational wave campaigns and radio surveys have the potential to unequivocally detect primordial black holes and therefore constrain their cosmic abundance.In this case, we have shown that even a small number of detections would place stringent constraints on self-annihilating dark matter.These constraints would be many orders of magnitude stronger than current constraints from gamma-ray telescopes [81][82][83][84], and even a single PBH detection would completely rule out thermal WIMPs, which constitute a substantial fraction of the Universe's dark matter.These scenarios thus rule out large regions of parameter space for proposed models of new physics at the weak scale (even those with subdominant WIMPs), making a compelling case for the search for primordial black holes now and in the near future.

ACKNOWLEDGMENTS
Where necessary, we have used the publicly available WebPlotDigitizer [85] to digitize plots.This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.B. J. K. acknowledges funding from the Netherlands Organization for Scientific Research (NWO) through the VIDI research program "Probing the Genesis of Dark Matter" (680-47-532).D. G. has received financial support through the Postdoctoral Junior Leader Fellowship Programme from la Caixa Banking Foundation (Grant No. LCF/BQ/LI18/ 11630014).We thank Christian Byrnes for sharing preliminary results from Ref. [16].

APPENDIX: GAMMA-RAY CONSTRAINTS
Here we give the details of our procedure for constraining ðσv rel Þ 0 in each PBH detection scenario using their diffuse gamma-ray emission and by treating them as gamma-ray point sources.For each detection scenario, as FIG. 2. Constraints on DM self-annihilation cross section.The solid lines correspond to the 95% C.L. upper limits obtained assuming a small number of PBH detections with LIGO/Virgo O3 (blue), the Einstein Telescope (ET, orange) and SKA (green).The lower dashed lines correspond to constraints which would be obtained if the number of PBH observations is as large as allowed by current limits.The dark gray region is the envelope of 95% C.L. profile likelihood contours for several supersymmetric models, while the light gray region is for singlet scalar scenarios.The horizontal dotted black line indicates the standard thermal relic cross section, 3 × 10 −26 cm 3 =s.The angled dotted black line shows the lower bound from unitarity for s-wave annihilation [20].
point estimates for f PBH , we conservatively use the fifth percentile for fPBH computed using the posteriors Pðf PBH jNÞ found above.These values are collected in Table II.
The key ingredient for both the diffuse and point-source analyses is the WIMP annihilation rate around a PBH.As described in the main text, the ultracompact minihalos (UCMHs) surrounding PBHs have ρ ∝ r −9=4 DM halos.Since the maximum-possible WIMP density at present is ρ max ¼ m χ =ðσv rel Þ 0 t 0 (where t 0 is the age of the Universe), 3the UCMH density profile is constant within the cutoff radius [16] The WIMP annihilation rate in a UCMH is obtained by integrating the UCMH density squared over all space4 We note that the maximum density (and the corresponding cutoff) apply only to the WIMP component of dark matter in the UCMH; we therefore now consider the case of a subdominant WIMP.In standard indirect detection analyses, the gamma-ray flux depends separately on the integral of the WIMP density squared in the observing region (the J-factor) and particle physics factors, giving an overall dependence on f 2 χ ðσv rel Þ 0 for underabundant species making up a fraction f χ of the DM.In the UCMH scenario, the dependence is more complicated, since the WIMP density profile depends on the particle physics through ðσv rel Þ 0 .Since r cut is the radius at which the WIMP density attains its maximum value and ρ χ ðrÞ ¼ f χ ρ DM ðrÞ, we have r cut ∝ f 4=9 χ .In contrast, ρ max does not depend on f χ .This means that the WIMP annihilation rate for an underabundant species scales as Finally, we note that we aim to constrain the zero-velocity limit of the WIMP annihilation cross section ðσv rel Þ 0 , which is typically the most relevant for annihilation signals today.We compare in Fig. 2 with the corresponding quantity calculated for different theories of weak-scale new physics.This comparison (along with our conclusions about the constraining power of a future PBH detection) therefore does not depend strongly on assumptions about the velocity dependence of the cross section.

PBHs as gamma-ray point sources
The natural test statistic for this analysis is N γ , the number of PBHs that pass the Fermi analysis cuts and appear in the 3FGL unassociated point-source list.Since there are 19 observed unassociated point sources, to set an upper bound on ðσv rel Þ 0 at the α ¼ 0.95 level we require the probability of observing at least 19 unassociated point sources to be equal to the corresponding p-value of 0.05.This probability is the cumulative distribution function for a binomial distribution: where N MW ≡ bf PBH ð3 × 10 12 M ⊙ Þ=M PBH c is the number of PBHs in the Milky Way, and p γ ððσv rel Þ 0 ; fPBH ; M PBH ; m χ Þ is the probability that a given PBH appears in the unassociated point-source catalogue.With all other parameters fixed, this defines the upper bound for ðσv rel Þ 0 , which must be determined numerically.We neglect the possibility that astrophysical sources appear in the unassociated point-source catalogue, since that would only strengthen our bounds.We compute p γ using a Monte Carlo procedure.First, a PBH's position is randomly sampled.The PBH spatial distribution is assumed to track the Galactic DM distribution, here taken to be the Einasto profile where r is the PBH's galactocentric distance and the halo parameters are ðρ s ; r s ; αÞ ¼ ð0.033 GeV=cm 3 ; 28.44 kpc; 0.17Þ [86].This parametrization is convenient, since the quantity r α can be efficiently sampled, because its PDF is a gamma distribution with scale θ r ¼ αr α s =2 and shape k r ¼ 3=α; the PBH's Galactic longitude and sine of its latitude are sampled uniformly.
Next, we assume that any PBH with sufficiently large integrated flux above 1 GeV lying far enough outside the Galactic plane appears in the unassociated point-source catalogue [67]: Using the annihilation rate defined above, the differential gamma-ray flux from the PBH is where dN=dE is the energy spectrum of photons per DM annihilation, which is easily integrated over energy.The final MC estimate for p γ is obtained by repeating this procedure and multiplying the fraction of sampled PBHs passing the detectability cuts by N MW .Since there are roughly N MW ¼ 1 × 10 5 -2 × 10 11 PBHs in the Milky Way (depending on N PBH and M PBH ), ðσv rel Þ 0 and thus p γ must be very small to give fewer than N U ¼ 19 point sources.A naive Monte Carlo simulation is extremely inefficient for such small cross sections, but with importance sampling, we can obtain more accurate results for lower computational cost.In detail, Eq. (A8) and the flux threshold from Eq. (A7) define the maximum distance d max at which a PBH is detectable, which means we should only sample PBHs within this distance of Earth.Denoting the distance from the Galactic Center to Earth as d ⊕ , we sample the PBH Galactic longitude uniformly from ½− arcsin d max d ⊕ ; arcsin d max d ⊕ and the sine of its Galactic latitude uniformly from ½−d max =d ⊕ ; d max =d ⊕ .Furthermore, inverse CDF sampling allows us to sample r from the interval ½d ⊕ − d max ; d ⊕ þ d max .The Monte Carlo samples must be reweighted to account for the restricted sampling volume.

Diffuse gamma rays from PBHs
The photons emitted from DM halos around PBHs located across all redshifts contribute to the extragalactic gamma-ray background (EGB).Taking the gamma-ray flux in each of the 26 energy bins Fermi uses for their EGB analysis to be normally distributed, the sum of the squared, standardized fluxes follows a χ 2 distribution, motivating the use of this quantity as the test statistic: The ith bin energies and error bars are denoted by E i and σ i and the differential fluxes in the numerator are the observed flux, flux from UCMHs surrounding PBHs and background flux respectively.
The PBH flux depends on the PBH and dark matter parameters and is well known: [15,86,87] The optical depth τðE; zÞ accounts for the attenuation of gamma rays emitted from redshift z and observed with energy E due primarily to pair production on baryonic matter and scattering off the extragalactic background light, for which we use the tables from Ref. [86].The cosmologically averaged DM density is ρ DM ≈ 30 M ⊙ =kpc 3 .Instead of modeling the astrophysical contributions to the EGB from, e.g., blazars, star-forming galaxies and misaligned active galactic nuclei, we set ϕ ex bg ðE i Þ ¼ ϕ ex obs ðE i Þ.This yields less stringent limits than Ref.
[87] and subsequent mixed DM-PBH analyses utilizing their results [15,16], since their fiducial EGB background model is in tension with the Fermi observations.With this simplification, the observed test statistic reduces to χ 2 obs ¼ P 26 i¼1 ½ϕ ex PBH ðE i Þ=σ i 2 .To set a limit at the level α ¼ 0.95 requires setting the probability of the test statistic exceeding the observed value to the corresponding p-value of 0.05: This implies χ 2 obs ¼ 38.9, the relevant critical value for 26 degrees of freedom.Substituting Eq. (A10) into the definition of χ 2 obs yields the bound on ðσv rel Þ 0 .In Fig. 3, we show separately the diffuse and pointsource constraints arising in each of the detection scenarios we have considered.
The point-source constraints are most important for small N PBH , while the diffuse constraints dominate for the values of N max in Table I.The difference in scaling is because the number of PBHs passing the integrated flux cut is proportional to f PBH times p γ (the probability of a PBH lying in the tail of the integrated flux distribution), which is a very sensitive function of Γ ∝ ðσv rel Þ 1=3 0 , since the distribution is roughly log-normal.This means that holding N γ fixed while increasing f PBH translates into a small decrease in ðσv rel Þ 0 .In contrast, the diffuse extragalactic flux from PBHs scales less strongly with ðσv rel Þ 0 as ϕ ex PBH ∝ f PBH ðσv rel Þ 1=3 0 , so the cross section constraint scales as f 3 PBH .These different scalings for the diffuse and point-source cross section bounds are exhibited in Fig. 3.

Dependence of results on priors
In Fig. 3 detection rate.Recall that to obtain the posterior for f PBH , we must choose a prior for the mean number of PBH candidates Λ, which we expect to observe.Solid lines in Fig. 3 show the projected bounds assuming a Jeffreys prior on the PBH rate.For a Poisson process, the Jeffreys prior (which is invariant under reparametrization) is given by PrðΛÞ ∝ 1=Λ 1=2 .Dot-dashed lines show instead the results assuming a log-flat prior on the mean PBH detection rate, PrðΛÞ ∝ 1=Λ.The Jeffreys prior is that assumed by LIGO/ Virgo [47,48] and corresponds to the results we show in Fig. 2. Instead, the log-flat prior pushes the posterior towards a smaller detection rate, smaller values of f PBH , and therefore weaker constraints on ðσv rel Þ 0 .We see that with only a single PBH detection in gravitational waves (upper-left and middle-left panels), the log-flat prior may weaken the possible constraints by up to an order of magnitude.This is because with a small number of detections, the posterior on f PBH is still rather broad (see Fig. 1), and so a change to a more conservative prior can noticeably lift the constraints.For a larger number of observations (right panels of Fig. 3), the results are insensitive to the choice of prior (the solid and dot-dashed lines are almost indistinguishable).This is to be expected; in this case, the posterior is dominated by the Poissonian likelihood for observing N max PBH candidates.
Even in the scenarios where a change of prior on the rate weakens the projected constraints, our main conclusions are not substantively changed: even a small number of PBH detections can place stringent constraints on thermal WIMP dark matter and, more generally, theories of new physics at the weak scale.

7 FIG. 3 .
FIG.3.Diffuse and point-source 95% C.L. upper limits on ðσv rel Þ 0 .Each row corresponds to a different PBH detection scenario: LIGO O3 (top), Einstein Telescope (middle), and SKA (bottom).Green and red lines show constraints derived from diffuse extragalactic emission and point sources, respectively.Solid lines show constraints assuming a Jeffreys prior on the rate of PBH observations, while dot-dashed lines show results for the more conservative log-flat prior.Note that in the right panels, the solid and dot-dashed lines are largely indistinguishable.In Fig.2, we show results for the Jeffreys prior, taking the stronger of the diffuse and point-source constraints at a given DM mass[20].

TABLE II .
Point estimates for f PBH .The entries in this table are the fifth percentiles of the posteriors Pðf PBH jNÞ for each scenario in TableI.The second column indicates the prior for the merger rate (in the case of gravitational wave scenarios) or the mean number of expected observations (for the SKA scenario).