Primordial Black Holes as Silver Bullets for New Physics at the Weak Scale

TL;DR

This paper discusses how detecting primordial black holes through radio and gravitational wave searches could impose strict constraints on weak-scale extensions of the Standard Model, potentially ruling out many dark matter theories.

Contribution

It demonstrates that upcoming PBH searches can significantly constrain or exclude broad classes of weak-scale new physics models, including popular supersymmetric and scalar dark matter theories.

Findings

Radio and gravitational wave searches can detect PBHs or set constraints.

Detection prospects could rule out large parts of supersymmetric models.

PBH observations can test and constrain weak-scale physics theories.

Abstract

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 waves 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 PBHs 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.

Figures

Tables

Table 1. Table 1: PBH Detection Scenarios. For each of the three detection scenarios we consider, we show the benchmark PBH mass M ⊙ subscript 𝑀 direct-product M_{\odot} which we assume, as well as N min subscript 𝑁 min N_{\mathrm{min}} and N max subscript 𝑁 max N_{\mathrm{max}} , the smallest and largest number of detections we consider. Note that N max subscript 𝑁 max N_{\mathrm{max}} is the largest number of detections which would be compatible with current PBH constraints.

	Sub-solar (O3)	High redshift (ET)	SKA
$M_{\mathrm{PBH}}$ [$M_{\odot }$]	0.5	10	100
$N_{\min }$	1	1	10
$N_{\max }$	80	24000	80

Table 2. Table 2: Point estimates for f PBH subscript 𝑓 PBH f_{\mathrm{PBH}} . The entries in this table are the 5th percentiles of the posteriors P ​ ( f PBH | N ) 𝑃 conditional subscript 𝑓 PBH 𝑁 P(f_{\mathrm{PBH}}|N) for each scenario in Table  1 . The second column indicates the prior for the merger rate (in the case of gravitational wave scenarios) or mean number of expected observations (for the SKA scenario).

	Prior	$0.5M_{\odot }$	$10M_{\odot }$	$100M_{\odot }$
$N_{\min }$	Log-flat	$1.49\times {10}^{-3}$	$7.28\times {10}^{-6}$	$2.56\times {10}^{-5}$
	Jeffreys	$2.81\times {10}^{-3}$	$1.32\times {10}^{-5}$	$2.82\times {10}^{-5}$
$N_{\max }$	Log-flat	$1.10\times {10}^{-1}$	$4.17\times {10}^{-3}$	$3.34\times {10}^{-4}$
	Jeffreys	$1.10\times {10}^{-1}$	$4.17\times {10}^{-3}$	$3.36\times {10}^{-4}$

Equations

$$\mathcal{R}(z)=\frac{1}{2}{n}_{\mathrm{PBH}}P(t_{\mathrm{merge}}=t[z]).$$

$$\Lambda ={\int }_{z_1}^{z_2}\mathcal{R}(z)\frac{\mathrm{d}⟨VT⟩}{\mathrm{d}z}\mathrm{d}z\equiv \overline{\mathcal{R}}⟨VT⟩,$$

$$P(\Lambda |N)\propto P(\Lambda )P(N|\Lambda )\propto \frac{1}{{\Lambda }^{1/2}}{\Lambda }^N\exp \left(-\Lambda \right).$$

$$\begin{array}{rl}P(\overline{\mathcal{R}}|N,⟨VT⟩)& \propto ⟨VT⟩{\left(\overline{\mathcal{R}}⟨VT⟩\right)}^{N-1/2}\exp \left(-\overline{\mathcal{R}}⟨VT⟩\right).\end{array}$$

$$P(f_{\mathrm{PBH}}|N)=P(\overline{\mathcal{R}}|N)\frac{\partial \overline{\mathcal{R}}}{\partial f_{\mathrm{PBH}}}.$$

$$⟨VT⟩=\frac{4\pi }{3}{D}_{\mathrm{avg}}^{3}T.$$

$$f_{\chi }{}^4(\sigma v_{\text{rel}}{)}_0\gtrsim 10^{-38}{\mathrm{cm}}^3{\mathrm{s}}^{-1}{\left(\frac{m_{\chi }}{1\mathrm{TeV}}\right)}^6.$$

$$r_{\mathrm{cut}}$$

$$\times {\left(\frac{m_{\chi }}{100\mathrm{GeV}}\right)}^{-4/9}{\left(\frac{M_{\mathrm{PBH}}}{M_{\odot }}\right)}^{1/3}.$$

$$\Gamma$$

$$\Gamma \propto {\left[f_{\chi }^4(\sigma v_{\text{rel}}{)}_0\right]}^{1/3}.$$

$$\mathrm{Pr}(N_{\gamma }\ge 19)=\sum _{N_{\gamma }=19}^{N_{\mathrm{MW}}}B(N_{\gamma }|N_{\mathrm{MW}},p_{\gamma })=0.05,$$

$${\rho }_{\mathrm{Ein}}(r)$$

$$|b|>20^{\circ },{\Phi }_{>1\mathrm{GeV}}>7\times 10^{-10}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}.$$

$${\varphi }_{\mathrm{PBH}}^{\mathrm{point}}(E)$$

$${\chi }^2=\sum _{i=1}^{26}{\left(\frac{{\varphi }_{\mathrm{obs}}^{\mathrm{ex}}(E_i)-[{\varphi }_{\mathrm{PBH}}^{\mathrm{ex}}(E_i)+{\varphi }_{\mathrm{bg}}^{\mathrm{ex}}(E_i)]}{{\sigma }_i}\right)}^2.$$

$${\varphi }_{\mathrm{PBH}}^{\mathrm{ex}}(E)$$

$$\mathrm{Pr}({\chi }^2\ge {\chi }_{\mathrm{obs}}^2)={\int }_{{\chi }_{\mathrm{obs}}^2}^{\infty }P({\chi }^2)\mathrm{d}{\chi }^2$$

Full text

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 center of galaxies Gondolo and Silk (1999); Gondolo (2000); Ullio et al. (2001); Bertone et al. (2001); Merritt et al. (2002); Bertone and Merritt (2005); Merritt et al. (2007) and intermediate-mass BHs Bertone et al. (2005); Zhao and Silk (2005); Bringmann et al. (2009). The argument has been more recently extended to the case of primordial black holes (PBHs), which can form before Big Bang nucleosynthesis Hawking (1971); Carr and Hawking (1974) 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_{\mathrm{PBH}}=\Omega_{\mathrm{PBH}}/\Omega_{\mathrm{DM}}$ Lacki and Beacom (2010); Eroshenko (2016); Boucenna et al. (2018); Adamek et al. (2019).

In this letter, 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\lesssim 1\,M_{\odot}$ with LIGO/Virgo; ii) The detection of GWs produced by the merger of $\mathcal{O}(10\,M_{\odot})$ 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 1.

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 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 weak-scale 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 \faGithub, and a link below each figure (\faFileCodeO) provides the python code with which it was generated. An archived version is available on Zenodo \faTags.

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_{\odot}$ Slany and Stuchlik (2008) or perhaps larger Shaposhnikov and Titarchuk (2009). In standard theories of stellar evolution, compact objects do not form below the Chandrasekhar mass of $\sim 1.4\,M_{\odot}$ Chandrasekhar (1931a, b, 1935). As a result, the detection of sub-solar mass BHs would strongly suggest a primordial origin (though see Refs. Shandera et al. (2018); Kouvaris et al. (2018) for alternative exotic formation mechanisms). The second case is the observation of BH-BH mergers at high redshift with proposed detectors such as Einstein Telescope (ET) Hild et al. (2011) and Cosmic Explorer (CE) Sathyaprakash et al. (2012). 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 Mapelli et al. (2019). Above around $z\sim 40$, the merger rate of astrophysical BHs should be negligible Koushiappas and Loeb (2017). High-redshift mergers therefore point towards a primordial BH population which formed much earlier Chen and Huang (2019).

PBH Merger Rates: Setting the stage.

To calculate the PBH merger rate as function of the PBH abundance, we follow Refs. Ali-Haïmoud et al. (2017); Kavanagh et al. (2018) (though see also Refs. Hayasaki et al. (2016); Raidal et al. (2017); Chen and Huang (2018); Ballesteros et al. (2018); Belotsky et al. (2019); Bringmann et al. (2019); Raidal et al. (2019)). In this scenario, pairs of PBHs decouple from the Hubble flow and form binaries before matter-radiation equality Nakamura et al. (1997); Ioka et al. (1998); Sasaki et al. (2016). For a given value of $(f_{\mathrm{PBH}},M_{\mathrm{PBH}})$, we calculate the comoving PBH merger rate (in the source frame) as a function of redshift as:

[TABLE]

Here, $n_{\mathrm{PBH}}=f_{\mathrm{PBH}}\Omega_{\mathrm{DM}}\rho_{c,0}/M_{\mathrm{PBH}}$ is the comoving number density of PBHs in terms of the critical density today $\rho_{c,0}$. The merger time distribution is $P(t_{\mathrm{merge}})$, which we evaluate at $t_{\mathrm{merge}}=t[z]$, the age of the Universe at redshift $z$. The merger time depends on the semi-major axis $a$ and eccentricity $e$ of the binary orbit Peters (1964). The distribution $P(t_{\mathrm{merge}})$ is thus obtained from the distribution of $P(a,e|M_{\mathrm{PBH}})$, calculated assuming a uniform initial distribution of PBHs and accounting for torques due to all other PBHs and smooth density perturbations Ali-Haïmoud et al. (2017). We also account for the impact of local DM halos on the binary properties Kavanagh et al. (2018), for which we assume a DM density profile with slope $\gamma=9/4$ Adamek et al. (2019), leading to a roughly 30% enhancement in the merger rate.

To obtain projected constraints on $\mathcal{R}$ (or equivalently $f_{\mathrm{PBH}}$) from future GW observation campaigns, we follow LIGO/Virgo Abbott et al. (2016a, b) and perform a Bayesian analysis111We will assume that the threshold for detection of mergers is high enough that there is no contamination from terrestrial ‘events’ (i.e. $p_{\mathrm{astro}}\simeq 1$ for all observations).. We assume that the number of observed mergers is Poisson distributed, $P(N|\Lambda)=\mathrm{Poiss}(\Lambda)$, with mean:

[TABLE]

where $\mathrm{d}{\langle{VT}\rangle}/\mathrm{d}z=T\frac{1}{1+z}\frac{\mathrm{d}V_{c}}{\mathrm{d}z}f(z)$ is the differential time-volume sensitivity for an observation time $T$ Abbott et al. (2016b, a). We calculate the comoving volume as a function of redshift $\mathrm{d}V_{c}/\mathrm{d}z$ Hogg (1999) assuming the Planck 2018 best fit cosmological parameters Aghanim et al. (2018). The selection function $f(z)$ encodes the detection probability for a given source at redshift $z$.

As in Refs. Abbott et al. (2016a, b), we assume a Jeffreys’ prior222The Jeffreys prior is an uninformative prior which is invariant under reparametrization of the parameter of interest Jeffreys (1946). on $\Lambda$ (though see Appendix A.3 for a more detailed discussion of the prior dependence of our results). Using Bayes’ Theorem Bayes (1764), the resulting posterior on $\Lambda$, given $N$ observations, is:

[TABLE]

The corresponding posterior on the redshift-averaged merger rate $\overline{\mathcal{R}}=\Lambda/{\langle{VT}\rangle}$ is then:

[TABLE]

We also incorporate a log-normal prior on the total time-volume sensitivity ${\langle{VT}\rangle}$, with a 30% relative uncertainty in order to account for possible calibration uncertainties Abbott et al. (2016b). We marginalize over ${\langle{VT}\rangle}$ to obtain the marginal posterior for $\overline{\mathcal{R}}$, which can be transformed into a posterior on $f_{\mathrm{PBH}}$:

[TABLE]

We now apply this formalism to the specific cases of sub-solar mass PBHs observed with LIGO and high-redshift mergers observed with ET.

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_{\mathrm{PBH}}=0.5\,M_{\odot}$.

LIGO O3 is sensitive to such mergers only in a small range of redshifts, $z\lesssim 0.02$. We therefore assume that the PBH merger rate is constant over this range, replacing it with the redshift-averaged value $\overline{\mathcal{R}}$. We estimate the sensitive time-volume using the approximation of Ref. Magee et al. (2018):

[TABLE]

Here, $T$ is the analyzable live-time of the detectors and $D_{\mathrm{avg}}$ is the average detector range. This is obtained from $D_{\mathrm{max}}$, the horizon distance of the detector, by averaging the detector response over both orientation and location of the binary, which gives $D_{\mathrm{avg}}\approx D_{\max}/2.26$. We assume a duration of one year for O3 with a 70% duty cycle Abbott et al. (2018a). We also assume a further reduction of 5% due to data quality cuts Abbott et al. (2018b). The LIGO O2 horizon distance for observing 0.5 $M_{\odot}$ mergers (with a signal-to-noise of at least 8 in each of the two detectors) is $\sim 60\,\mathrm{Mpc}$ Abbott et al. (2018c). From Ref. Abbott et al. (2018a), we expect an improvement in the horizon distance of about 50% from O2 to O3. We therefore assume a horizon distance $D_{\mathrm{max}}=90\,\mathrm{Mpc}$ ($z\approx 0.0192$), giving a sensitive time-volume of $\langle VT\rangle\approx 1.8\times 10^{-4}\,\mathrm{Gpc}^{3}\,\mathrm{yr}$.

The resulting posteriors on $f_{\mathrm{PBH}}$ are shown in blue in Fig. 1. The grey shaded region shows the current constraint on $0.5\,M_{\odot}$ PBHs, coming from the observed merger rate during LIGO O2 Abbott et al. (2018c). We note that the median value of $f_{\mathrm{PBH}}$ scales as $\sqrt{N}$. This is because the merger rate scales predominantly as $f_{\mathrm{PBH}}^{2}$, reflecting the number of PBH binaries which are formed in the early Universe.

High-redshift mergers.

In the case of high-redshift mergers, we consider the sensitivity of the proposed Einstein Telescope (ET) Hild et al. (2011) 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 mono-chromatic mass function, now at $M_{\mathrm{PBH}}=10\,M_{\odot}$, 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 $\mathcal{R}(z)$. The differential time-volume sensitivity $\mathrm{d}{\langle{VT}\rangle}/\mathrm{d}z$ is calculated as in Eq. (2) with the selection function $f(z)$ estimated from Ref. Hall (2019a, b): this corresponds to $\sim 5\%$ of mergers detectable at $z=40$, decreasing to zero close to $z=100$.

The resulting posteriors on $f_{\mathrm{PBH}}$ are shown in orange in Fig. 1, where we again cut off the posterior at the current upper limit (grey) from Ref. Kavanagh et al. (2018).

Discovering Galactic PBHs with SKA.

Isolated black holes can in principle be detected by observing the broad-band spectrum of radiation — spanning from radio waves all the way up to the X-ray band — produced by the interstellar gas accreting onto them Fender et al. (2013). 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 Gaggero et al. (2017); Ivanov et al. (2019). 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 Manshanden et al. (2018).

During the next decade, the Square Kilometre 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 percent) of the Dark Matter in the Universe Bull et al. (2018).

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  Gaggero et al. (2017); Manshanden et al. (2018); Bull et al. (2018); Ivanov et al. (2019). Here we assume that, as a result of such data analysis, a certain number of radio sources observed by SKA is clearly associated to 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\,{\rm MHz}$–$14\,{\rm GHz}$ domain. We consider $100$ M⊙ PBHs because these massive objects would be more easily detected, given the $\propto 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^{\circ}<l<0.7^{\circ};-0.3^{\circ}<b<0.3^{\circ}$), where a large concentration of molecular gas is present.

As a first step, we determine the distribution $P(N_{\mathrm{SKA}}|f_{\mathrm{PBH}})$ of the number of PBHs detectable by SKA1-MID (1.4 GHz), adopting the Monte Carlo simulations described in Gaggero et al. (2017); Manshanden et al. (2018), and assuming gain G = 15 K/Jy, receiver temperature $T_{\rm rx}$ = 25 K, sky temperature towards the GC $T_{\rm sky}$ = 70 K, bandwidth $\Delta\nu$ = 770 MHz, and $\simeq 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_{\mathrm{PBH}}|N_{\mathrm{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_{\mathrm{PBH}}$. The current upper bound for $100\,M_{\odot}$ PBHs comes from the radio and X-ray searches mentioned above Manshanden et al. (2018).

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 Bertschinger (1985) and simple 1D simulations Mack et al. (2007) predict that these ultra-compact minihalos (UCMHs) have $\rho(r)\propto r^{-9/4}$ density profiles, which has been confirmed by recent 3D simulations Adamek et al. (2019). Since $f_{\mathrm{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 Bertoni et al. (2015); Hooper and Witte (2017); Schoonenberg et al. (2016) have identified 19 bright, high-latitude, non-variable unassociated point sources that are spectrally compatible with annihilating WIMPs. As described in detail in Appendix A, we perform a Monte Carlo simulation to assess the observability of UCMHs by Fermi. We then use this to determine the 95% confidence level (CL) upper bound on the WIMP annihilation cross-section in the zero-velocity limit $(\sigma v_{\text{rel}})_{0}$. This upper limit depends on the PBHs’ spatial distribution which we assume tracks the Milky Way DM distribution. We fix $f_{\mathrm{PBH}}$ to the 5th percentile of the posterior $P(f_{\mathrm{PBH}}|N)$, 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 $(\sigma v_{\text{rel}})_{0}$ by comparing with the expected number of UCMHs passing cuts on their integrated gamma-ray flux and galactic latitude (given $M_{\mathrm{PBH}}$, $m_{\chi}$ and $N$).

Annihilation in UCMHs outside the Milky Way over all redshifts contributes to the diffuse, isotropic extragalactic background (EGB) Taylor and Silk (2003); Ullio et al. (2002); Ando (2005), which has been measured by Fermi Ackermann et al. (2015a). 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_{\mathrm{PBH}}$ to the 5th 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_{\chi}=\Omega_{\chi}/\Omega_{\mathrm{DM}}$ for under-abundant thermal relics. Typically, the DM annihilation rate depends on the combination $f_{\chi}{}^{2}(\sigma v_{\text{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 $(\sigma v_{\text{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_{\chi}{}^{4}(\sigma v_{\text{rel}})_{0}$; this is derived in Appendix A.

Results and discussion.

For each detection scenario in Table 1 we show as function of WIMP mass the 95% CL upper limit on $f_{\chi}^{4}(\sigma v_{\text{rel}})_{0}$ in Fig. 2, where $f_{\chi}=\Omega_{\chi}/\Omega_{\mathrm{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 $\bar{b}b$. 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 Appendix A for a more detailed comparison of the limits).

We find that a detection of $\mathcal{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 grey the envelope of the 95% CL profile-likelihood contours for the MSSM7 Athron et al. (2017a) and various GUT-scale SUSY models Athron et al. (2017b) obtained by the GAMBIT collaboration. These contours show the zero-velocity limit of the annihilation cross section $(\sigma v_{\text{rel}})_{0}$ and so can be directly compared with our projected gamma-ray constraints assuming $s$-wave annihilation only. In light grey, we also show contours for scalar singlet scenarios Athron et al. (2017c). 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_{\chi}{}^{4}(\sigma v_{\text{rel}})_{0}=10^{-29}{\mathrm{~{}cm}}^{3}/{\mathrm{s}}$, and in the leftmost sliver of the singlet contour and throughout the dark grey envelope all such models have $f_{\chi}{}^{4}(\sigma v_{\text{rel}})_{0}>10^{-32}{\mathrm{~{}cm}}^{3}/{\mathrm{s}}$.

Importantly, since $f_{\chi}$ scales approximately as $\langle\sigma v\rangle^{-1}$, the unitarity limit on the freeze-out annihilation cross-section $\langle\sigma v\rangle$ Griest and Kamionkowski (1990) yields, in the case of $s$-wave annihilation, a lower bound on the rescaled annihilation cross-section today,

[TABLE]

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, co-annihilation and kinematic thresholds but for the models shown this is always counterbalanced by values of $f_{\chi}$ 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 Gómez-Vargas et al. (2017); Steffen (2006), they do affect asymmetric DM models Petraki and Volkas (2013); Zurek (2014), 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.

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 currents constraints from gamma-ray telescopes Abramowski et al. (2011); Ackermann et al. (2012, 2015b); Abdallah et al. (2018) 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 sub-dominant WIMPs), making a compelling case for the search for primordial black holes now and in the near future.

Acknowledgements.

Where necessary, we have used the publicly available WebPlotDigitizer Rohatgi (2018) to digitise plots. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. BJK acknowledges funding from the Netherlands Organization for Scientific Research (NWO) through the VIDI research program “Probing the Genesis of Dark Matter” (680-47-532). DG has received financial support through the Postdoctoral Junior Leader Fellowship Programme from la Caixa Banking Foundation (grant n. LCF/BQ/LI18/11630014). We thank Christian Byrnes for sharing preliminary results from Adamek et al. (2019).

Appendix A Gamma ray constraints

Here we give the details of our procedure for constraining $(\sigma v_{\text{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 point estimates for $f_{\mathrm{PBH}}$ we conservatively use the 5th percentile for $\hat{f}_{\mathrm{PBH}}$ computed using the posteriors $P(f_{\mathrm{PBH}}|N)$ found above. These values are collected in Table 2.

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 $\rho\propto r^{-9/4}$ DM halos. Since the maximum-possible WIMP density at present is $\rho_{\mathrm{max}}=m_{\chi}/(\sigma v_{\text{rel}})_{0}\,t_{0}$ (where $t_{0}$ is the age of the universe),333This is conservative since $\rho_{\mathrm{max}}$ was not derived for particles moving along radial trajectories Lacki and Beacom (2010). the UCMH density profile is constant within the cutoff radius Adamek et al. (2019)

[TABLE]

The WIMP annihilation rate in a UCMH is obtained by integrating the UCMH density squared over all space 444Changing the bound of the radial integral to the minihalo’s tidal radius (see Adamek et al. (2019), Sec. 2-E) does not appreciably change the annihilation rate.

[TABLE]

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 sub-dominant 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_{\chi}^{2}(\sigma v_{\text{rel}})_{0}$ for under-abundant species making up a fraction $f_{\chi}$ of the DM. In the UCMH scenario, the dependence is more complicated since the WIMP density profile depends on the particle physics through $(\sigma v_{\text{rel}})_{0}$. Since $r_{\mathrm{cut}}$ is the radius at which the WIMP density attains its maximum value and $\rho_{\chi}(r)=f_{\chi}\rho_{\mathrm{DM}}(r)$, we have $r_{\mathrm{cut}}\propto f_{\chi}^{4/9}$. In contrast, $\rho_{\mathrm{max}}$ does not depend on $f_{\chi}$. This means that the WIMP annihilation rate for an under-abundant species scales as

[TABLE]

Finally, we note that we aim to constrain the zero-velocity limit of the WIMP annihilation cross section $(\sigma v_{\text{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.

A.1 PBHs as gamma-ray point sources

The natural test statistic for this analysis is $N_{\gamma}$, 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 $(\sigma v_{\text{rel}})_{0}$ at the $\alpha=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:

[TABLE]

where $N_{\mathrm{MW}}\equiv\lfloor f_{\mathrm{PBH}}(3\times 10^{{12}}~{}M_{\odot})/M_{\mathrm{PBH}}\rfloor$ is the number of PBHs in the Milky Way and $p_{\gamma}((\sigma v_{\text{rel}})_{0},\hat{f}_{\mathrm{PBH}},M_{\mathrm{PBH}},m_{\chi})$ 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 $(\sigma v_{\text{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_{\gamma}$ 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

[TABLE]

where $r$ is the PBH’s galactocentric distance and the halo parameters are $(\rho_{s},r_{s},\alpha)=(0.033{\mathrm{~{}GeV/cm}}^{3},~{}28.44{\mathrm{~{}kpc}},~{}0.17)$ Cirelli et al. (2011). This parameterization is convenient since the quantity $r^{\alpha}$ can be efficiently sampled because its PDF is a gamma distribution with scale $\theta_{r}=\alpha r_{s}^{\alpha}/2$ and shape $k_{r}=3/\alpha$; 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 Hooper and Witte (2017):

[TABLE]

Using the annihilation rate defined above, the differential gamma-ray flux from the PBH is

[TABLE]

where $\mathrm{d}N/\mathrm{d}E$ is the energy spectrum of photons per DM annihilation, which is easily integrated over energy. The final MC estimate for $p_{\gamma}$ is obtained by repeating this procedure and multiplying the fraction of sampled PBHs passing the detectability cuts by $N_{\mathrm{MW}}$.

Since there are roughly $N_{\mathrm{MW}}=1\times 10^{{5}}-2\times 10^{{11}}$ PBHs in the Milky Way (depending on $N_{\mathrm{PBH}}$ and $M_{\mathrm{PBH}}$), $(\sigma v_{\text{rel}})_{0}$ and thus $p_{\gamma}$ must be very small to give fewer than $N_{\mathcal{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. 15 and the flux threshold from Eq. 14 define the maximum distance $d_{\mathrm{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_{\oplus}$, we sample the PBH galactic longitude uniformly from $[-\arcsin\frac{d_{\mathrm{max}}}{d_{\oplus}},\arcsin\frac{d_{\mathrm{max}}}{d_{\oplus}}]$ and the sine of its galactic latitude uniformly from $[-d_{\mathrm{max}}/d_{\oplus},d_{\mathrm{max}}/d_{\oplus}]$. Furthermore, inverse CDF sampling allows us to sample $r$ from the interval $[d_{\oplus}-d_{\mathrm{max}},d_{\oplus}+d_{\mathrm{max}}]$. The Monte Carlo samples must be reweighted to account for the restricted sampling volume.

A.2 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 $\chi^{2}$ distribution, motivating the use of this quantity as the test statistic:

[TABLE]

The $i$th bin energies and error bars are denoted by $E_{i}$ and $\sigma_{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: Boucenna et al. (2018); Ando and Ishiwata (2015); Cirelli et al. (2011)

[TABLE]

The optical depth $\tau(E,z)$ accounts for 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 Cirelli et al. (2011). The cosmologically averaged DM density is $\rho_{\mathrm{DM}}\approx 30~{}M_{\odot}/{\mathrm{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 $\phi_{\mathrm{bg}}^{\mathrm{ex}}(E_{i})=\phi_{\mathrm{obs}}^{\mathrm{ex}}(E_{i})$. This yields less stringent limits than Ref. Ando and Ishiwata (2015) and subsequent mixed DM-PBH analyses utilizing their results Boucenna et al. (2018); Adamek et al. (2019) since their fiducial EGB background model is in tension with the Fermi observations. With this simplification the observed test statistic reduces to $\chi_{\mathrm{obs}}^{2}=\sum_{i=1}^{26}\left[\phi_{\mathrm{PBH}}^{\mathrm{ex}}(E_{i})/\sigma_{i}\right]^{2}$.

To set a limit at the level $\alpha=0.95$ requires setting the probability of the test statistic exceeding the observed value to the corresponding $p$-value of 0.05:

[TABLE]

This implies $\chi_{\mathrm{obs}}^{2}=38.9$, the relevant critical value for 26 degrees of freedom. Substituting Eq. 17 into the definition of $\chi^{2}_{\mathrm{obs}}$ yields the bound on $(\sigma v_{\text{rel}})_{0}$.

In Fig. 3, we show separately the diffuse and point source constraints arising in each of the detection scenarios we have considered.

The point source constraints are most important for small $N_{\mathrm{PBH}}$ while the diffuse constraints dominate for the values of $N_{\mathrm{max}}$ in Table 1. The difference in scaling is because the number of PBHs passing the integrated flux cut is proportional to $f_{\mathrm{PBH}}$ times $p_{\gamma}$ (the probability of a PBH lying in the tail of the integrated flux distribution), which is a very sensitive function of $\Gamma\propto(\sigma v_{\text{rel}})_{0}^{1/3}$ since the distribution is roughly log-normal. This means that holding $N_{\gamma}$ fixed while increasing $f_{\mathrm{PBH}}$ translates into a small decrease in $(\sigma v_{\text{rel}})_{0}$. In contrast, the diffuse extragalactic flux from PBHs scales less strongly with $(\sigma v_{\text{rel}})_{0}$ as $\phi_{\mathrm{PBH}}^{\mathrm{ex}}\propto f_{\mathrm{PBH}}(\sigma v_{\text{rel}})_{0}^{1/3}$, so the cross section constraint scales as $f_{\mathrm{PBH}}^{3}$. These different scalings for the diffuse and point-source cross section bounds are exhibited in Fig. 3.

A.3 Dependence of results on priors

In Fig. 3, we also show the projected constraints on $(\sigma v_{\text{rel}})_{0}$ for two different choices of prior on the mean PBH detection rate. Recall that to obtain the posterior for $f_{\mathrm{PBH}}$, we must choose a prior for the mean number of PBH candidates $\Lambda$ 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 $\mathrm{Pr}(\Lambda)\propto 1/\Lambda^{1/2}$. Dot-dashed lines show instead the results assuming a log-flat prior on the mean PBH detection rate, $\mathrm{Pr}(\Lambda)\propto 1/\Lambda$. The Jeffreys’ prior is that assumed by LIGO/Virgo Abbott et al. (2016a, b) 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_{\mathrm{PBH}}$ and therefore weaker constraints on $(\sigma v_{\text{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_{\mathrm{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_{\mathrm{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.