Probing Primordial-Black-Hole Dark Matter with Gravitational Waves

TL;DR

This paper explores the potential of advanced LIGO to detect or constrain primordial black holes in the 20-100 solar mass range as dark matter candidates by analyzing gravitational wave data and mass spectra.

Contribution

It demonstrates that aLIGO can place strong limits on primordial black hole dark matter in the 20-100 solar mass range using gravitational wave observations.

Findings

aLIGO can detect PBH contribution down to f_PBH<0.5 within 5 years.

Current constraints are less robust in the 20-100 M_sun range.

aLIGO data can test the hypothesis that PBHs constitute all dark matter.

Abstract

Primordial black holes (PBHs) have long been suggested as a candidate for making up some or all of the dark matter in the Universe. Most of the theoretically possible mass range for PBH dark matter has been ruled out with various null observations of expected signatures of their interaction with standard astrophysical objects. However, current constraints are significantly less robust in the 20 M_sun < M_PBH < 100 M_sun mass window, which has received much attention recently, following the detection of merging black holes with estimated masses of ~30 M_sun by LIGO and the suggestion that these could be black holes formed in the early Universe. We consider the potential of advanced LIGO (aLIGO) operating at design sensitivity to probe this mass range by looking for peaks in the mass spectrum of detected events. To quantify the background, which is due to black holes that are formed from dying stars, we model the shape of the stellar-black-hole mass function and calibrate its amplitude to match the O1 results. Adopting very conservative assumptions about the PBH and stellar-black-hole merger rates, we show that ~5 years of aLIGO data can be used to detect a contribution of >20 M_sun PBHs to dark matter down to f_PBH<0.5 at >99.9% confidence level. Combined with other probes that already suggest tension with f_PBH=1, the obtainable independent limits from aLIGO will thus enable a firm test of the scenario that PBHs make up all of dark matter.

Equations

$$R_{\mathrm{PBH}}\approx 2f_{\mathrm{PBH}}^{53/21}{\left(M_{\mathrm{PBH}}/30M_{\odot }\right)}^{-11/21}\mathrm{Gp}{\mathrm{c}}^{-3}\mathrm{y}{\mathrm{r}}^{-1},$$

$$P(M_1)=A_{M_1}{M}_1^{-\alpha }\mathcal{H}(M_1-M_{\mathrm{gap}})e^{-M_1/M_{\mathrm{cap}}},$$

$$P(M_2)=A_{M_2}\mathcal{H}(M_2-M_{\mathrm{gap}})\mathcal{H}(M_1-M_2).$$

$$\frac{d{N}^{\mathrm{BG}}(M_1)}{d{M}_1}$$

$$N_i^{\mathrm{BG}}={\int }_{M_{\text{min},i}}^{M_{\text{max},i}}\frac{dN(M_1)}{d{M}_1}d{M}_1.$$

$$S/N=\sqrt{\sum _i{\left(\frac{N_i^{\mathrm{Tot}}-N_i^{\mathrm{BG}})}{\sqrt{N_i^{\mathrm{Tot}}}}\right)}^2};N^{Tot}\equiv N_i^{\mathrm{PBH}}+N_i^{\mathrm{BG}}$$

$$⟶\sqrt{\sum _i{\left(\frac{N_i^{\mathrm{PBH}}(f_{\mathrm{PBH}},M_{\mathrm{PBH}})}{\sqrt{N_i^{BG}}}\right)}^2}-n_{\sigma }=0,$$

Full text

One of the cornerstones of $\Lambda{\rm CDM}$, the concordance cosmological standard model, is the cold dark matter (DM) component that makes up $\sim\!25\%$ of the energy density in the Universe today. While the evidence for its existence are compelling Ade:2015xua ; Freese:2017idy , the nature of it is still unknown. As the limits on models of particle dark matter (in particular weakly-interacting massive particles, known as WIMPs Jungman:1995df ) are tightening Essig:2012yx ; Aprile:2012nq ; Akerib:2015rjg , it is becoming ever more important to consider alternative models.

An especially intriguing candidate to make up the invisible form of matter in the Universe is primordial black holes (PBHs), which are black holes that are formed deep in the radiation era of the infant Universe Carr:1974nx ; Meszaros:1974tb ; Chapline:1975 ; Carr:1975qj ; Clesse:2015wea . Based on various observations, the contribution of PBHs to dark matter has been strongly constrained across more than 30 orders of magnitude of their theoretically possible mass range Carr:2016drx ; Kuhnel:2017pwq ; Carr:2017jsz . Still, in several mass windows existing constraints are less stringent and additional probes are called for. This is especially true for the $20\,M_{\odot}\lesssim M_{\rm PBH}\lesssim 100\,M_{\odot}$ window, which has attracted much interest as a result of the first detection of merging black holes with measured masses of $\sim 30\,M_{\odot}$ by the LIGO observatory Abbott:2016blz , following the demonstration in Ref. Bird:2016dcv that the predicted merger rate for PBHs in this mass range is consistent with the estimated event rate for high mass mergers from the O1 aLIGO data TheLIGOScientific:2016pea .

One could describe the search for PBH dark matter in analogy to the one for particle dark matter, as is illustrated in Fig. 1. The constraints on the former to date have been solely based on “direct” detection searches, involving possible interactions between PBHs and standard astrophysical objects. In this Letter, we consider the prospects of the “indirect” search path for PBH dark matter, namely the production of standard (cosmological) model signals in the form of gravitational waves as a result of PBH self-interaction (or “annihilation”). The key to this approach is to understand and quantify the background as well as possible, and to identify unique features in the dark matter signal that can tell them apart.

Examples of such features are the orbital eccentricity of the coalescing binary and the black hole spins. The former was investigated recently in Ref. Cholis:2016kqi , but unfortunately the prospects for detecting events with a non-zero trace of the initial eccentricity in aLIGO are quite dim. As for the spin, the problem is twofold. The initial spin distribution of PBHs is unclear on the one hand (see Chiba:2017rvs for some estimates), and on the other, the key identifier of PBH mergers in this regard—that the spin of the two black holes should not be aligned—is mimicked by various models of dynamical binary formation Rodriguez:2016vmx ; Talbot:2017yur (whose rate is currently very uncertain), and so it is very hard to distinguish them from the background.

Here we focus on the mass spectrum of merging black hole binaries as a probe of PBH dark matter. The idea is simple: if PBHs in a given mass range account for all (or some fraction) of the dark matter in the Universe, we should see an excess of merger events involving black holes in the corresponding mass bin Kovetz:2016kpi . Provided that this excess is large enough to be differentiated from the expected background from mergers of black holes of stellar origin, a detection, or lack thereof, could either support the existence or place limits on the abundance of PBHs. This abundance is parameterized by the quantity $f_{\rm PBH}$, the fraction of dark matter in PBHs.

Our goal is to provide the most reliable constraints on PBH dark matter—focusing especially on the most motivated question which is whether all of dark matter can be explained by PBHs—and we therefore choose to be very careful and consider the most pessimistic case for PBH dark-matter detection, i.e. the lowest estimated rate of PBH mergers and the highest rate of stellar black hole mergers. We will show that even in this challenging case, gravitational-wave observations by aLIGO at design sensitivity should within a decade either exhibit strong hints for a PBH contribution to dark matter, or rule them out as the single form of dark matter (albeit allowing for a considerable fraction of it to made up of PBHs).

To begin, we review the suggested mechanisms of PBH binary formation. The first model was put forward in Refs. Nakamura:1997sm ; Sasaki:2016jop , where it was demonstrated that a subset of an initial PBH population, assumed to be randomly distributed in space, would be in close enough proximity to overcome the cosmic expansion and form bound pairs. The distribution of the semi-major axis of these binaries will be quite wide. However, as they generically have high initial orbital eccentricity (they avoid a head-on collision due to the influence of the closest neighbors), some will have merger times that allow them to reach the endpoint of their coalescence within the detectable volume of aLIGO. Ref. Sasaki:2016jop predicts a very high merger rate (which is in fact already in tension with existing observations if $f_{\rm PBH}\gg 0.01$). However, a crucial question is whether these early-formed binaries can survive as bound pairs throughout the evolution of the Universe, without being disrupted. This was briefly addressed in Ref. Sasaki:2016jop , where the probability for disruption was calculated to be as small as $\mathcal{O}(10^{-7})$, but only for a binary that resides in a Milky-Way type halo today. Since PBHs binaries in this scenario were created very early on, before the formation of large dark matter halos (in a series of violent cosmic processes), this is at best an underestimate of the actual disruption rate. A thorough reexamination of this model, taking into account the interaction of PBH binaries with other PBHs, the rest of dark matter (if $f_{\rm pbh}\ll 1$) and with baryonic matter, is under way, but is outside the scope of this Letter. In short, Ref. Ali-Haimoud:2017 evaluates the chance for disruption by tidal effects of the smooth halo and encounters with other PBHs (which are more efficient in the first halos, which are denser and have lower velocity dispersions), finding that it is higher than previously estimated. Meanwhile, the effect of a circumbinary accretion disk has been examined in Ref. Hayasaki:2009ug , concluding that it could lead to a decrease in the semi-major axis on a timescale fast enough such that all early-formed binaries with masses in the stellar-mass range will have merged well before redshift $z\sim 1$, and therefore remain outside the reach of aLIGO. Aiming to provide the most robust bounds on PBH dark matter, we shall therefore treat this rate as an optimistic case, and proceed to focus on more conservative scenarios, described below.

A second model of PBH binary formation was proposed in Ref. Bird:2016dcv . In this model, PBHs form binaries in close two-body encounters, as a result of energy loss from gravitational-wave emission as they pass each other by. The rate for this process to occur was calculated for dark matter halos of different masses (and densities and velocity dispersions), and when integrated over a full mass function (cutting off at the low mass end where halos would be too small not to evaporate by dynamical relaxation), the total merger rate was found to be

[TABLE]

where we have used the fact that the evaporation time is primarily determined by the number of black holes in the halo, as their density in this limit—and thus the dynamical time—is roughly constant with mass. This rate for PBHs with mass $M_{\rm PBH}=30\,M_{\odot}$ is consistent with the recent LIGO $90\%$-interval estimate for black holes as massive as the ones in the first detection, $0.5-12\,{\rm Gpc^{-3}yr^{-1}}$ TheLIGOScientific:2016pea . Of the many assumptions in this calculation, by far the most daring is the extrapolation of the halo mass function, velocity dispersion distribution, density profile and mass-concentration relation to very low ($\sim 10^{3}\,M_{\odot}$) dark matter halo masses, orders of magnitude below what can be observed or even simulated. The good news however, is that as these binaries have a very small separation when they form, they merge very quickly and are not susceptible to the effect of interfering processes such as mentioned above for early-formed binaries. While a bias against a more-optimistic higher rate is not too worrisome for the purposes of this work, it is more crucial to have a clear idea for the lowest reasonable bound on the merger rate. In work to appear Nishikawa , a merger rate is calculated for the same two-body formation mechanism of PBH binaries, but focusing on the contribution to the merger rate from encounters that occur inside dark matter density spikes around supermassive black holes (SMBHs) at the centers of galaxies Gondolo:1999ef . Integrating over the mass function of SMBHs from $\sim 10^{5}\,M_{\odot}$ to $\sim 10^{9}\,M_{\odot}$, this is found to yield roughly $>10\%$ of the total rate from all halos, providing a highly-conservative lower floor for the total PBH merger rate.

We will therefore proceed by adopting $R_{\rm PBH}$ in Eq. (1) as a conservative estimate for the rate of PBH mergers in the local Universe, addressing the pessimistic case of a rate ten times lower, and bearing in mind that if early-Universe binaries are somehow found to be stable to disruption, the rate could be order(s) of magnitude higher.

The next step is to make a prediction for the background, which consists of mergers of stellar black holes that contribute to the mass spectrum of detected events. To do this, we require an assumption for the astrophysical merger rate (which depends on factors such as metallicity and merger time-delay distributions) and for the mass function of stellar black holes. For the merger rate, we follow Ref. TheLIGOScientific:2016htt . We use a simple approximation, $R(z)\simeq 97(1+z)^{2}\,{\rm Gpc^{-3}yr^{-1}}$, which provides a good fit to their results at redshifts $z<1$ TheLIGOScientific:2016pea . As for the black hole mass function, we follow Ref. Kovetz:2016kpi and consider a simple ansatz whereby the probability distribution function (PDF) of the black hole mass is described by a simple power law (motivated by the slope of the stellar initial mass function, which has been corroborated by numerous observations in the $1-100\,M_{\odot}$ mass range Salpeter (1955); Kroupa (2001)). We impose sharp and exponential cutoffs at the lower and upper end of the mass spectrum, respectively, to take into account the neutron-star—BH transition threshold and the increasing wind-driven mass loss of high mass stars in the Wolf-Rayet phase. Denoting the black hole mass by $M_{1}$, the PDF is given by

[TABLE]

where $\mathcal{H}$ is the heaviside function and $A_{M_{1}}$ is an overall normalization. We shall use as fiducial values $\alpha=2.35$ (following Kroupa (2001)), $M_{\rm gap}=5\,M_{\odot}$ (motivated by current observations Bailyn et al. (1998); Ozel et al. (2010); Farr et al. (2011) and by some theoretical works Belczynski et al. (2012); Fryer et al. (2012); Kochanek (2014)) and $M_{\rm cap}=60\,M_{\odot}$. The latter choice naturally affects our constraints strongly at masses beyond the cutoff, as with a vanishing background, PBHs of that mass would be easier to detect. However, very massive PBH binaries will have a merger frequency too small to be detected by aLIGO, which offsets this effect. We deliberately avoid choosing a lower value, to reflect our ignorance and so as not to artificially strengthen the bound on $f_{\rm PBH}$. With future observations, it might be possible to improve on these somewhat arbitrary choices, perhaps leading to more tightened constraints.

Our observable is the total number of detected events as a function of the black hole mass. In order to determine the background contribution, we need to model the mass ratio between the black hole binaries as well. Setting $M_{1}$ from now on to be the mass of the heavier black hole in each binary, we follow Ref. TheLIGOScientific:2016pea and assume that the mass of the lighter black hole, $M_{2}$, has a uniform distribution in the range $[M_{\rm gap},M_{1}]$, given by

[TABLE]

The observable redshift volume of gravitational waves from BH mergers depends on the instrumental properties. The signal-to-noise ratio ($S/N$) for a single interferometer detector is given by $({S/N})^{2}=\int_{f_{\rm min}}^{f_{\rm max}}df{4h_{c}^{2}(f)}/5{S_{n}(f)(2f)^{2}}$, where $h_{c}(f)$ is the observed strain amplitude and $S_{n}(f)=h_{n}^{2}(f)$ is the strain noise amplitude (for more details, see Ref. Cholis:2016kqi and references therein). We follow convention and set the detection threshold at $S/N>8.0$ Abadie et al. (2010). We use the approximated analytical model for aLIGO noise of Ref. Ajith (2011), setting $f_{\rm min}=10\,{\rm Hz}$, above which the curve matches the official LIGO curve very well Shoemaker (2010).

Lastly, inferring the mass of the merging black holes from the gravitational-wave signal involves an associated uncertainty. We model this by convolving the mass function with a log-normal distribution reflecting a $5\%$ relative mass error for aLIGO observations (see Ref. Kovetz:2016kpi ). This choice suggests a minimal width for our binning of the mass function when calculating the signal-to-noise.

We are now equipped to make a theoretical prediction for the total number of detected background merger events over a time $T_{\rm obs}$ with a given mass $M_{1}$ (the mass of the heavier BH in the one-dimensional (1D) case), or two masses $M_{1},M_{2}$ (in the 2D case), which is given by

[TABLE]

Here the observable redshift volume is defined by $z_{\rm max}$, the maximum redshift up to which a BH merger with masses $M_{1},M_{2}$ can be detected; $H(z)$ is the Hubble parameter and $\chi(z)$ is the radial comoving distance. In the 2D case, we use $dN(M_{1},M_{2})/dM_{1}dM_{2}$, dropping the first integration in Eq. (5). Integrating within each mass bin $i$ with edges $[M_{\textrm{min},i},M_{\textrm{max},i}]$, we finally get

[TABLE]

We then divide $N(M_{1})$ into $50$ logarithmic bins from $4\,M_{\odot}$ to $120\,M_{\odot}$ and likewise $N(M_{1},M_{2})$ into 100 bins (10 along each axis). The bin width is chosen such that the mass measurement error is subdominant, and the variance in each bin is set by the Poisson error $\sigma_{i}^{2}=N_{i}$.

The most powerful way to constrain the abundance of a PBH population with a narrow mass distribution is to look at the 2D mass distribution of detected events, and check for a peak in the mass bin(s) surrounding the central value of that distribution. This approach makes use of all the available data, but since the background needs to be calculated for each mass bin, making accurate forecasts requires a good understanding of the black-hole mass function as well as the binary mass ratio. To somewhat relax this model dependence, we can also choose to limit ourselves to using only the number counts for the heavier BH in each binary (the mass-ratio will still enter the signal-to-noise calculation through the fiducial choice of the $M_{2}$ PDF, but since each column in the one-dimensional $M_{1}$ distribution is an integral quantity over the full $M_{2}$ mass range, the dependence on this choice is weaker in the 1D case). The widely used convention when forecasting limits on the fraction of dark matter in PBHs is to assume a delta-function PBH mass function. In practice, we assume Gaussian PDFs for $M_{\rm PBH}=M^{\rm PBH}_{1}=M^{\rm PBH}_{2}$ with a $5\%$ width (the precise choice is unimportant as long as this width is smaller than the measurement error). Together with the rate in Eq. (1), we get $N^{\rm PBH}_{i}(f_{\rm PBH},M_{\rm PBH})$ for each value of $f_{\rm PBH}$ and $M_{\rm PBH}$ using the prescription in Eqs. (5),(6).

We now derive our forecast for the limits that aLIGO at design sensitivity can impose with the planned $6$ years of observation by solving the following equation for $f_{\rm PBH}$

[TABLE]

[TABLE]

where we set a desired signal-to-noise ratio of $n_{\sigma}=3$ or $5$ standard deviations and have assumed under the null hypothesis ($N_{\rm PBH}=0$) that the (Poisson) error in each bin is $\sqrt{N^{\rm BG}_{i}}$ (using the Gaussian approximation for the Poisson distribution is valid given the number of events in each bin). The result—in the form of $3-$ and $5\!-\!\sigma$ limits on $f_{\rm PBH}$ for each PBH mass—is shown in Fig. 2. We see that based on the rate in Eq. (1), the scenario in which all the dark matter is in the form of PBHs can be strongly tested (ruled out at $\gg 5\sigma$) when using either the full 2D or 1D mass spectra of observed BH mergers.

Note that we have assumed the stellar-BH mass-function parameters can be held fixed, rather than fitting for them in tandem with the amplitude of the PBH contribution and marginalizing over them. Under the assumption that it is smooth, consistent with TheLIGOScientific:2016pea , the effect of this approximation on our results is small. We emphasize that if it is found to deviate strongly from Eq. (2), future data from planned next generation GW experiments such as Cosmic Explorer Evans:2016mbw , Einstein Telescope ET , DECIGO Seto:2001qf and LISA Amaro-Seoane et al. (2017)—which will be sensitive enough at low frequencies to be able to detect BH mergers at high redshifts, well beyond the peak of the star-formation rate— will allow a straightforward discrimination between the stellar and primordial merging-BH populations, based on their very different redshift distributions Koushiappas:2017kqm . Another source of uncertainty is the overall background amplitude, whose current $90\%$-confidence range still spans roughly an order of magnitude TheLIGOScientific:2016pea .

To incorporate the various modeling uncertainties, we use the dependence of Eq. (1) on $f_{\rm pbh}$ and of Eq. (7) on the background and PBH rates, to show in Fig. 2 a band encompassing an underestimate (overestimate) of up to a factor of $400\%$ ($200\%$) in the background (signal) rates. As can be seen, our conclusion that the scenario of PBH DM can be tested convincingly ($\gg 5\sigma$) in the range $20\,M_{\odot}\lesssim M_{\rm PBH}\lesssim 100\,M_{\odot}$ is quite robust. even with a PBH merger rate ten times lower, which as described above should be taken as an ultra-conservative bound, $f_{\rm PBH}=1$ can still be rejected at $5\sigma$ confidence across the $20-100\,M_{\odot}$ range using the 2D information. Naturally, if the rate is much higher Sasaki:2016jop ; Ali-Haimoud:2017 , a null detection will yield even stronger limits on $f_{\rm PBH}$, motivating efforts to better understand the mechanisms of PBH binary formation and disruption, using both analytic methods and simulations. Eventual bounds will depend on additional real-world factors, such as the operating sensitivity of aLIGO during its six year run. In particular, we stress that the constraints can be significantly improved if the signal-to-noise threshold for a detected event in each interferometer is reduced, as should be done when performing a statistical analysis of an ensemble of events Kovetz:2016kpi .

While the constraints we forecast are weaker than some that have been already claimed or projected by other methods, they constitute a unique and independent method of testing the scenario, by focusing on the self-interaction of PBHs rather than their interaction with other astrophysical objects. They are subject to different systematics and modeling assumptions, and since the analysis presented here heavily errs on the side of caution and still finds promising results, they represent a truly robust test, achievable within a decade, of the important cosmological scenario that dark matter is made of PBHs.

We thank Yacine Ali-Haïmoud, Simeon Bird, Marc Kamionkowski and especially Ilias Cholis for discussions. This work was supported by NSF Grant No. 0244990, NASA NNX15AB18G and the Simons Foundation.