Counting Black Holes: The Cosmic Stellar Remnant Population and Implications for LIGO

TL;DR

This paper develops an empirical model linking stellar evolution remnants to gravitational wave black hole mergers, predicting host galaxy populations and informing expectations for future LIGO detections.

Contribution

It introduces a novel empirical approach to interpret black hole merger signals based on galaxy stellar mass and metallicity relationships, connecting stellar remnants to gravitational wave observations.

Findings

Milky Way hosts millions of ~30 M_sun black holes.

Dwarf galaxies like Draco host ~100 such black holes.

Merger rate densities from LIGO are consistent with ~1% black hole binary involvement.

Abstract

We present an empirical approach for interpreting gravitational wave signals of binary black hole mergers under the assumption that the underlying black hole population is sourced by remnants of stellar evolution. Using the observed relationship between galaxy mass and stellar metallicity, we predict the black hole count as a function of galaxy stellar mass. We show, for example, that a galaxy like the Milky Way should host millions of $\sim 30~M_\odot$ black holes and dwarf satellite galaxies like Draco should host $\sim 100$ such remnants, with weak dependence on the assumed IMF and stellar evolution model. Most low-mass black holes ($\sim10 M_\odot$) typically reside within massive galaxies ($M_\star \simeq 10^{11} M_\odot$) while massive black holes ($\sim 50~M_\odot$) typically reside within dwarf galaxies ($M_\odot \simeq 10^9 M_\odot$) today. If roughly $1\%$ of black holes are involved in a binary black hole merger, then the reported merger rate densities from Advanced LIGO can be accommodated for a range of merger timescales, and the detection of mergers with $> 50~M_\odot$ black holes should be expected within the next decade. Identifying the host galaxy population of the mergers provides a way to constrain both the binary neutron star or black hole formation efficiencies and the merger timescale distributions; these events would be primarily localized in dwarf galaxies if the merger timescale is short compared to the age of the universe and in massive galaxies otherwise. As more mergers are detected, the prospect of identifying the host galaxy population, either directly through the detection of electromagnetic counterparts of binary neutron star mergers or indirectly through the anisotropy of the events, will become a realistic possibility.

Equations

$$N_{\mathrm{bh}}(>m_{\mathrm{bh}},M_{\star })$$

$$=N_{\star }(M_{\star })\int \mathcal{P}(Z,M_{\star }){\int }_{{\mathcal{M}}_{\min }(m_{\mathrm{bh}},Z)}^{{\mathcal{M}}_{\mathrm{u}}}\xi ({\mathcal{M}}^{\prime })d{\mathcal{M}}^{\prime }dZ.$$

$$\bar{\mathcal{M}}(M_{\star })={\int }_{0.08M_{\odot }}^{{\mathcal{M}}_l(M_{\star })}{\mathcal{M}}^{\prime }\xi ({\mathcal{M}}^{\prime })d{\mathcal{M}}^{\prime }.$$

$$n_{\mathrm{bh}}(>m_{\mathrm{bh}})={\int }_{M_{\min }}^{\infty }\varphi (M_{\star })N_{\mathrm{bh}}(>m_{\mathrm{bh}},M_{\star })d{M}_{\star }.$$

$$n_{\mathrm{bh}}^{\mathrm{pair}}(>m_{\mathrm{bh}})=\frac{1}{2}ϵn_{\mathrm{bh}}(>m_{\mathrm{bh}}).$$

$${\dot{n}}_{\mathrm{bh}}(>m_{\mathrm{bh}},t)=n_{\mathrm{bh}}(>m_{\mathrm{bh}})\frac{\psi (t)}{{\int }_0^{t_0}\psi (t^{\prime })\mathrm{d}{t}^{\prime }}.$$

$$\mathcal{R}=\frac{1}{2}ϵ{\int }_0^{t_0}{\dot{n}}_{\mathrm{bh}}(t_0-{\tau }^{\prime })\mathrm{P}({\tau }^{\prime })\mathrm{d}{\tau }^{\prime },$$

$$\mathcal{R}=\frac{1}{2}ϵ{\dot{n}}_{\mathrm{bh}}(t_0-\tau ).$$

$${\mathcal{R}}_{\mathrm{prompt}}=\frac{1}{2}ϵ{\dot{n}}_{\mathrm{co},0}.$$

$${\dot{n}}_{\mathrm{ns},0}={\int }_{M_{\min }}^{\infty }\varphi (M_{\star })\int \frac{{\dot{M}}_{\star }(M_{\star },Z_g)}{\bar{\mathcal{M}}}\mathcal{P}(Z_g,M_{\star })$$

$${\int }_{{\mathcal{M}}_{\min }^{\mathrm{ns}}}^{{\mathcal{M}}_{\max }^{\mathrm{ns}}}\xi ({\mathcal{M}}^{\prime })d{\mathcal{M}}^{\prime }\mathrm{d}{Z}_g\mathrm{d}{\mathrm{M}}_{\star }.$$

Full text

††pagerange: Counting Black Holes: The Cosmic Stellar Remnant Population and Implications for LIGO–References††pubyear: 2016

1 Introduction

With the first detection of gravitational waves, the Laser Interferometer Gravitational-Wave Observatory (LIGO) simultaneously confirmed a fundamental prediction of General Relativity and discovered the first known binary black hole (BBH) merger (Abbott et al., 2016d). This first gravitational wave event, GW150914, appears to have been caused by the merger to two fairly massive ($\sim 30M_{\odot}$) black holes. Subsequent detections of a BBH mergers (GW151226 and GW170104, Abbott et al., 2016c, 2017) and a candidate BBH event (LVT151012, Abbott et al., 2016a) have allowed more robust estimates of the local BBH merger rate density and have confirmed the existence of black holes involved in these mergers with masses that range from 7.5 to 36 $M_{\odot}$.

As the field now pivots from gravitational wave discovery to gravitational wave astronomy, there are number of questions we hope to explore in more detail. One basic question is the origin of these massive black holes. Heavy ($\gtrsim 30\ M_{\odot}$) black holes are expected to exist as the result of stellar evolution (e.g. Spera et al., 2015; Belczynski et al., 2010, and references therein), and have been predicted to dominate the LIGO signal (e.g. Belczynski et al., 2016b). However, the possibility that the GW150914 event was due to primordial black holes (that would constitute some part of the dark matter) has also been advanced (Bird et al., 2016; Cholis et al., 2016; Carr et al., 2016; Inomata et al., 2016).

Ab initio computation of the BBH merger rate (of stellar remnants) is currently not possible. This calculation requires inputs from multiple fields including galaxy formation and numerical relativity. In this work, we outline a simple way to compute the BBH merger rate for stellar remnant black holes that allows one to assess the uncertainties in the various required ingredients in a transparent manner. The key idea is that we have a good empirical understanding of the overall galaxy number density, stellar ages and metallicities as a function of galaxy mass, and estimates of the initial mass function of stars that formed in these galaxies. Using these ingredients, together with the current generation of stellar evolution codes, we can provide a grounded estimate of the global distribution of black holes as a function of black hole mass and galaxy stellar mass. With this as a starting point, we are able to quantify the astrophysical parameters needed to produce the observed BBH merger signals observed. For example, the LIGO collaboration has reported a global “event-based” merger rate of $\mathcal{R}=55^{+185}_{-46}$ Gpc*-3* yr*-1* for binaries more massive than $5\,M_{\odot}$ each (Abbott et al., 2016b, f), and the detection of GW170104 has reduced the range to $12\leq\mathcal{R}\leq 213$ (Abbott et al., 2017). Below we will demonstrate that such a rate is reasonable with the population of stellar black holes we expect to exist within galaxies in the local universe and we discuss how the host mass of mergers will be a valuable diagnostic for testing scenarios going forward.

Our work is complementary to past work by Belczynski et al. (2016b) and Lamberts et al. (2016), who focus their efforts on understanding the formation of binaries and the details of binary black hole evolution. Both of these papers explicitly focused on the first GW150914 event. Lamberts et al. (2016) concluded that the black holes involved in GW150914 likely formed in a massive galaxy at $z\sim 1$, but that formation in a dwarf galaxy was also likely possible. Belczynski et al. (2016b) suggested that the black holes likely form in low-metallicity systems. Chatterjee et al. (2017) explored the formation of BBH systems specifically in globular clusters and came to qualitatively similar conclusions.

Important to all of these investigations is the realization that massive black hole formation is suppressed in stellar populations with higher metallicites (Spera et al., 2015; Belczynski et al., 2008). For example, according to the calculations of Spera et al. (2015), a star of mass $\mathcal{M}\simeq 90M_{\odot}$ will be required to produce a remnant of mass $m_{\rm bh}=30M_{\odot}$ if its metallicity is $Z/Z_{\odot}=-0.5$. A more metal poor star ($-1.5$) will need to be only $\mathcal{M}\simeq 33M_{\odot}$ to produce a $30M_{\odot}$ BH remnant. These expectations, combined with the long merger times often predicted for isolated BBH systems (see Postnov & Yungelson, 2014, and references therein), have led many authors to conclude that the detected BBH merger signals observed by LIGO will be dominated by the mergers of black holes formed in the early universe. However, processes such as the Kozai-Lidov mechanism (Kozai, 1962; Lidov, 1962) may accelerate the merger timescale both in galaxy centers(VanLandingham et al., 2016) and star clusters (Silsbee & Tremaine, 2017; Kimpson et al., 2016) and low-metallicity star formation is ongoing at low-redshifts, especially in low mass galaxies (Ellison et al., 2008; Mannucci et al., 2010; Lara-López et al., 2010), allowing massive BBH systems to form locally and potentially merge quickly enough to detect. We note that spin constraints limit the merging time of GW150914 to be greater than 100 Myr (Kushnir et al., 2016). Given the large uncertainties, we treat the merging time of BBH systems as a free parameter in our analysis, and explore observational means to constrain the timescale.

Our work is organized as follows: In §2 we describe our approach. In §3 we estimate the local number density of black holes as a function of black hole mass using the local mass-metallicity relation and stellar mass function. We examine how these black holes are distributed among galaxies of a given mass in §3.1 and work out the global number density in §3.2. We move from the total number density to merger rate densities in §4 and discuss how by localizing event rates by host galaxy masses we can begin to constrain more details of BBH merger scenarios in §4.1. We summarize and conclude in §5.

2 Our Approach

Before detailing our predictions in the next section, we here provide an overview of our basic approach.

Our calculation proceeds in two steps. First we determine the cosmic number density of stellar remnants, $n_{\rm bh}$, using empirical measurements of the galaxy population (§3). Next, with this empirically-grounded number density in hand, we calculate binary black hole merger rate densities using two simple parameters that quantify our ignorance of binary formation and binary evolution (§4). Specifically, we introduce a dimensionless parameter $\epsilon\leq 1$ that captures the complicated physics of binary system formation and binary remnant survival. We also introduce a characteristic merger timescale, $\tau$, that parameterizes our uncertainty in binary black hole evolution and the distribution of times it takes for a binary pair to merge after formation. In principle, any merger timescale distribution can be mapped to this characteristic timescale $\tau$.

As we show below, $n_{\rm bh}$ is fairly robust to uncertainties. The larger uncertainties are in the physics of binary formation and evolution encapsulated in $\epsilon$ and $\tau$. We show that the observed LIGO signals can be readily explained with reasonable values for both of these parameters and make predictions for future observations of black hole and neutron star merger rates that provide a means to constrain the physics they encapsulate. We also explore and quantify the degeneracies between $\epsilon$ and $\tau$ and describe how determining the galaxy host mass distribution for mergers could be a way to break these degeneracies.

3 Black hole number density predictions

In this section, we discuss the global number density of black holes at the present time and their distribution as a function of the galaxy stellar mass. In §3.1 we work out the expected specific black hole frequency as a function of black hole mass and galaxy mass. In §3.2 explore the implied cosmic density of black holes. These determinations will ground us as we move forward to estimate merger rates.

3.1 Black Hole Populations Within Galaxies

The number of black holes more massive than $m_{\rm bh}$ that exist within a galaxy will depend on the number of massive stars previously formed in that galaxy with an initial mass larger than some minimum value, $\mathcal{M}_{\rm min}(m_{\rm bh},Z)$. The minimum mass of a star required to produce a remnant black hole of some mass $m_{\rm bh}$ is expected to be a strong function of stellar metallicity $Z$ owing to mass loss from stellar winds. Figure 1 plots $\mathcal{M}_{\rm min}$ as a function of $Z$ for three example black hole remnant masses $m_{\rm bh}>10,30,$and $50M_{\odot}$, as determined by stellar evolution codes PARSEC (Spera et al., 2015, black) and StarTrack (Fryer et al., 2012, blue). Both calculations give similar results, especially at low metallicities111The largest discrepancy between the models is at $m_{\rm bh}>50M_{\odot}$, which is perhaps not unexpected since the fits from Fryer et al. (2012) are extrapolations at this mass range.. We see that for high metallicities ($Z\gtrsim-1.5$), a very large stellar progenitor ($\mathcal{M}\gtrsim 90M_{\odot}$) will be required to produce the massive black holes of the type that have been observed in mergers by LIGO. Lower metallicity populations require less extreme progenitors. We will adopt the PARSEC results as our fiducial choice below.

With $\mathcal{M}_{\rm min}$ in hand, we can determine the total number of black holes more massive than $m_{\rm bh}$ that have ever formed, $N_{\rm bh}(>m_{\rm bh})$, within a galaxy of mass $M_{\star}$ and a total number of stars $N_{\star}(M_{\star})$ by integrating over the stellar IMF $\xi(\mathcal{M})$ and the metallicity distribution function (MDF) of stars expected for a galaxy of that mass $\mathcal{P}(Z,M_{\star})$:

[TABLE]

We set the upper limit on the IMF integral at $\mathcal{M}_{\rm u}=150\ M_{\odot}$, though our results are not strongly sensitive to this choice.222Setting the upper limit to $\infty$ in the subsequent analysis changes our results by $<10\%$. The black hole count is normalized by $N_{\star}(M_{\star})=M_{\star}/\bar{\mathcal{M}}(M_{\star})$, where

[TABLE]

For the upper limit $\mathcal{M}_{l}(M_{\star})$, we chose the stellar mass with main sequence lifetime equal to the average stellar age of galaxies of mass $M_{\star}$ (from Behroozi et al., 2013, see their Fig.  13). For $\mathcal{P}(Z,M_{\star})$ assume that galaxies more massive than $M_{\star}=10^{9}M_{\odot}$ follow a log-normal distribution in $Z$, with mean and standard deviation given by Gallazzi et al. (2005). For smaller galaxies, we use the results of Kirby et al. (2013), who measured resolved-star MDFs for 15 individual local dwarf galaxies with stellar masses $M_{\star}\simeq 10^{3}-10^{8}M_{\odot}$. We assume that these individually observed MDFs are representative for galaxies in the dwarf mass range throughout the universe. Finally, for $\xi(\mathcal{M})$ we adopt a Kroupa (2002) IMF for our fiducial calculations. We have also explore the effects of metallicity-dependent IMF (specifically adopting the IMF of Geha et al., 2013) and find that our results are sensitive at the factor of $\sim 2$ level to this level of variation in the IMF.

Figure 2 shows $N_{\rm bh}(M_{\star})$ as derived from Equation 1 for three choices of black hole mass: $m_{\rm bh}>10,\,30,\,50M_{\odot}$. The $m_{\rm bh}>30M_{\odot}$ results are shown as squares, color coded according the median metallicity of galaxies at each $M_{\star}$. Boxes at $M_{\star}<10^{9}M_{\odot}$ are placed at the stellar masses of the individual galaxies in the Kirby et al. (2013) MDF sample. We see that for galaxies less massive than $M_{\star}\simeq 10^{10}M_{\odot}$, the number of black holes of all masses scales linearly with galaxy mass, $N_{\rm bh}\propto M_{\star}$. For example, we find that there should be roughly one $30M_{\odot}$ black hole per $1000M_{\odot}$ of stars in a galaxy, at least for smaller galaxies. For the $m_{\rm bh}>30$ and $50M_{\odot}$ populations, the linear scaling with $M_{\star}$ breaks down when galaxies become so metal rich that only the low-metallicity tail of the population can be associated with massive black hole formation. But the black hole counts recover and begin increasing monotonically with stellar mass again once galaxies become massive enough that there is no longer a strong trend between $M_{\star}$ and $Z$ (at $M_{\star}\gtrsim 10^{11}M_{\odot}$). Note for the smallest black holes, $N_{\rm bh}\propto M_{\star}$ for all galaxy masses, as there is very little trend between progenitor mass and remnant mass for $m_{\rm bh}\lesssim 10M_{\odot}$ (see Figure 1).

3.2 Cosmic Black Hole Number Density

In order to obtain the global number density of black holes, $n_{\rm bh}(>m_{\rm bh})$, we integrate $N_{\rm bh}(M_{\star})$ over the galactic stellar mass function $\phi(M_{\star})$:

[TABLE]

We adopt the results of Baldry et al. (2012) for $\phi(M_{\star})$ though we have checked that using the stellar mass function from Bernardi et al. (2013) does not change our results significantly. For the minimum mass in the $M_{\star}$ integral we use M${}_{\rm min}=10^{3}M_{\odot}$ and find the number density of black holes is convergent below this galaxy mass.

The solid black line in Figure 3 shows the results of this calculation of $n_{\rm bh}$ for our fiducial Kroupa IMF assumption. For comparison, the dotted line shows the result for the Geha metallicity-dependent IMF (Geha et al., 2013). Though the Geha IMF yields slightly more black holes, the factor of $\sim 2$ offset is not large given the other uncertainties in this calculation. We will adopt the Kroupa IMF in all the results to follow. In that case, we see, for example, that number density of 30 $M_{\odot}$ black holes is $n_{\rm bh}\sim 10^{14}$ Gpc*-3*. If $\sim 0.1\%$ of these black holes merge over a Hubble time ($t_{H}\sim 10^{10}$ yrs) then we might expect a local rate of $\mathcal{R}\sim 0.001\,n_{\rm bh}/t_{H}\sim 10$ Gpc*-3* yr*-1*, which is comparable to the LIGO estimate for massive black holes based on the $\sim 30~{}M_{\odot}$ pair seen in the GW150914 event ($\mathcal{R}_{30}=3.4^{+8.6}_{-2.8}$ Gpc*-3* yr*-1* Abbott et al., 2016b). In Section 4 we will present a more careful comparison to the inferred LIGO rates.

Figure 3 clearly shows that the overall black hole number density in the universe is fairly high. Whether this provides a consistent and reasonable explanation of the LIGO observations depends largely on the expected fraction of merging BBH and the merger time scale. One question of interest is how is this cosmic abundance of black holes distributed among galaxies? Figure 4 shows the results for various cuts on $m_{\rm bh}$. We see that most low-mass lack holes in the universe reside within massive galaxies, while higher mass black holes tend to reside in dwarfs. This general trend is expected since low-mass black holes tend to track stellar mass, and most of the stellar mass in the local universe is in massive galaxies. Massive black holes tend to reside in $M_{\star}\sim 10^{8-10}M_{\odot}$ galaxies. The most likely host for a single $\sim 30M_{\odot}$ black hole chosen at random in the universe is a galaxy of stellar mass $M_{\star}\sim 10^{10}M_{\odot}$. Of course, just because most black holes live in massive galaxies this does not necessarily imply that most black hole mergers will occur in massive galaxies. We will return to this question in §4.1.

3.3 Comparison to Core Collapse Supernova Rates

A useful test of our methodology is to compare the observed density of core collapse supernova (CCS) remnants to that predicted in our model. If we change the mass limits in Equation 1, we can calculate the global density of CCS remnants using Equation 3. For the minimum stellar mass we use $8~{}M_{\odot}$ and for the upper limit we use 18 $M_{\odot}$, which assumes that most stars above this mass collapse to form black holes with no visible supernovae (Smartt, 2015). Doing this gives a value of roughly one CCS for every 100 solar masses of stars formed, and integrating over the stellar mass function yields a remnant density of $n_{r}=2.5\times 10^{6}$ Mpc*-3*. Current measurements place the local CCS rate density at $0.7\times 10^{-4}$ yr*-1* Mpc*-3* (Strolger et al., 2015). If we assume the CCS rate closely tracks the star formation rate, we can normalize the evolution of the cosmic star formation rate density from Madau & Dickinson (2014) to this value and integrate over the age of the universe to find the total density of CCS remnants. This gives a density of $n_{r}=3.8\times 10^{6}$ Mpc*-3*, which is within a factor of $\sim 1.5$ of our estimated remnant density.

4 Black Hole Merger Rates

In what follows we will assume that black hole mergers occur among binary pairs and that these mergers occur after the birth of the binary pair over a timescale $\tau$. For simplicity, we will focus on merger rates for pairs of black holes each with masses above the same threshold value of $m_{\rm bh}$. The merger timescale $\tau$ is subject to several assumptions and therefore difficult to calculate from first principles (Lipunov et al., 1997; Sipior & Sigurdsson, 2002; Dominik et al., 2013; Belczynski et al., 2016a; Lamberts et al., 2016). Our approach is to treat $\tau$ as a parameter to be constrained.

At any given epoch, the number density of black hole pairs available to merge before $z=0$ can be written in terms of the black hole number density at that time. Specifically for pairs of mass $m_{1},m_{2}>m_{\rm bh}$ we have

[TABLE]

Here we have introduced a new parameter that we refer to as the “binary black hole efficiency”: $\epsilon\equiv f_{{\rm b}\star}\times f_{\mathrm{m_{1}/m_{2}}}\times f_{\rm surv}\times f_{\rm t}<1$. This dimensionless quantity parameterizes our ignorance of merging black holes from massive stars. The value $f_{{\rm b}\star}$ is the massive star binary fraction ($f_{{\rm b}\star}\sim 0.5$; e.g. Sana et al., 2012; Kobulnicky & Fryer, 2007; Pfalzner & Olczak, 2007) and $f_{\mathrm{m_{1}/m_{2}}}$ is the fraction of massive binary systems with mass ratios near unity. Current models predict $f_{\mathrm{m_{1}/m_{2}}}\sim 0.1$ for $\mathrm{m_{1}/m_{2}=0.9}$ (Sana et al., 2012). The fraction of those massive star binaries that survive as black hole pairs after stellar evolution is $f_{\rm surv}\sim 0.1$ (Belczynski et al., 2016a; Lamberts et al., 2016). Finally, $f_{\rm t}$ represents the fraction of binary black holes with orbital configurations that make them available to merge before the present day ($f_{\rm t}<1$). In this work we assume $\epsilon$ is independent of mass $\mathrm{m_{1},m_{2}}$. If it varied significantly in the $10-50M_{\odot}$ mass range, then our predictions for BBH mergers not yet observed by LIGO would be inaccurate. With these assumptions, we find below that the binary efficiency parameter values $\epsilon\simeq 0.01-0.001$ can reproduce the reported black hole merger rate density from LIGO using only stellar remnant black holes.

The formation rate density of black hole pairs that can merge will depend on the birthrate density of black holes: $\dot{n}_{\rm bh}^{\rm pair}=0.5\,\epsilon\,\dot{n}_{\rm bh}$. Here, the over-dot implies differentiation with respect to time. We will assume that the black hole formation rate density tracks the observed shape of the global star formation rate (SFR) density $\psi(t)$ (with $t=t_{0}=13.7$ Gyr corresponding to the present day) such that

[TABLE]

For $\psi(t)$ we used the parameterization of Madau & Dickinson (2014). The SFR density peaks at $z\sim 2$, corresponding to $t\simeq 3.4$ Gyr after the Big Bang and a lookback time of $10.3$ Gyr.

Now let us assume that for every binary black hole pair that is born that there is an distribution of times $\mathrm{P(\tau^{\prime})}$ for them to merge. In this case, the cosmic black hole merger rate density today ($t=t_{0}$) can be written as an integral over the black hole birth rate density:

[TABLE]

where $\dot{n}_{\rm bh}$ is evaluated at the black hole mass of relevance for the merger rate. We note that $\mathrm{P(\tau^{\prime})}$ is the average distribution of merger times; the full distribution depends on many other underlying factors such as the orbit of the binary system and the environment it is in.

For simplicity, we treat $\mathrm{P}(\tau^{\prime})$ as a delta function centered on a characteristic timescale: $\mathrm{P}(\tau^{\prime})=\delta(\tau^{\prime}-\tau)$. This allows for our results to be cast in terms of two effective parameters: the merging efficiency $\epsilon$ and the characteristic timescale $\tau$, and results in a present day BBH merger rate density given by:

[TABLE]

Note that Equations 7 and 5 imply that for a fixed value of $\epsilon$, a merger timescale that matches the lookback time to the peak in cosmic star formation ($\tau\sim 10$ Gyr) will produce the largest local merger rate. Thus, in order to match the observed local merger rate, a case with $\tau\sim 10$ Gyr will require the smallest values of $\epsilon$.

The rationale for this simple single-timescale approach is that it allows us to readily explore the relationships between merger timescales, the unknown binary merger efficiency, and the host galaxies of merging events. Though the assumption is clearly a major simplification, any physically-motivated $\tau^{\prime}$ distribution can in principle be mapped to a delay time $\tau$. For example, one well-motivated assumption for the merger time distribution is $\mathrm{P(\tau^{\prime})\propto 1/\tau^{\prime}}$ (Dominik et al., 2013). Using $\tau=1$ Gyr yields roughly the same density, with nearly the same constraints on $\epsilon$, as assuming $\mathrm{P}(\tau^{\prime})\propto 1/\tau^{\prime}$ in Equation 6. Even if underlying merger timescale distribution is multi-modal, a combination of delta-function models can be used. For example, Belczynski et al. (2016a) and Lamberts et al. (2016), predict a bimodal distribution in birth times of massive BBH merger progenitors, with one peak at redshifts $z\sim 2$ and the other at much lower redshift, $z\sim 0.2$. In this case, a combination of delta functions, one with $\tau\sim 1$ and another with $\tau\sim 10$ Gyr reproduces such a model.

Figure 5 shows the predicted local merger rate of $m_{\rm bh}>5M_{\odot}$ black holes as a function of merger timescale $\tau$ for two choices of our binary efficiency parameter $\epsilon=0.01$ and $0.001$. The shaded band shows the total observed range from Abbott et al. (2017): $12\leq\mathcal{R}\leq 213$. We see that for shorter timescales ($\tau\lesssim 2$ Gyr), $\epsilon=0.01$ matches the data better. For longer timescales (close to the peak of cosmic star formation, $\tau\simeq 10$ Gyr) the lower efficiency of $0.1\%$ is more consistent with the measurement. Note that Abbott et al. (2016b) also quote a event-based rate for binary $m_{\rm bh}\simeq 30M_{\odot}$ mergers like GW150914 of $\mathcal{R}_{30}=3.4^{+8.6}_{-2.8}$ Gpc*-3* yr*-1*. Our predicted binary merger rates for $m_{\rm bh}>30M_{\odot}$ black holes also agree well with their $\mathcal{R}_{30}$ range, producing curves like those in Figure 5 shifted down by approximately an order of magnitude, with $\epsilon=0.01$ working better for $\tau\lesssim 2$ Gyr and $\epsilon=0.001$ consistent for $\tau\simeq 10$ Gyr (just as in the $m_{\rm bh}>5M_{\odot}$ rate case).

The degeneracy between $\tau$ and $\epsilon$ values is clearer in Figure 6. The band shows the range of parameter values that are consistent with the reported LIGO rates for merging pairs with $m_{\rm bh}>5M_{\odot}$. For $\tau\lesssim 2$ Gyr, efficiencies of $\epsilon\simeq 0.002-0.03$ are required. The efficiencies need to be smaller if the typical merger timescale approaches the lookback time of peak star formation $\tau\simeq 10$ Gyr, $\epsilon\simeq 0.0002-0.004$. The sharp uptick in required efficiency as $\tau\rightarrow t_{0}=13.7$ Gyr is driven by the fact that the star formation rate drops to zero as we approach the big bang. As the merger timescale approaches age of the Universe, reproducing the observed rates requires virtually every black hole that is present in the early universe to end up merging today.

Figure 7 displays our predicted merger rates for black holes of various masses ($\mathrm{m_{bh}}\geq 5$, $30$, and $50~{}M_{\odot}$) as a function of $\tau$ for $\epsilon=0.01$. As previously discussed, for $\tau<4$ Gyr, this choice of $\epsilon$ is consistent with the reported merger rate for $>5M_{\odot}$ BBH mergers, though the overall amplitude of the lines is linearly proportional to $\epsilon$. The aim of this figure is to illustrate how the rates vary with compact object mass. The $m_{\rm bh}>30M_{\odot}$ BBH merger rate density, for example, is lower by a factor of $\sim 8$ at fixed $\tau$. The the $2-\sigma$ limit from Abbott et al. (2016b) for massive black holes of this kind is $\mathcal{R}_{30}=0.6-12.4~{}\mathrm{Gpc^{-3}yr^{-1}}$, which matches our predictions for this choice of $\epsilon$ as long as $\tau<3$ Gyr (with larger $\tau$ requiring smaller $\epsilon$ as in Figure 6).

Figure 7 also includes neutron star-neutron star (NS-NS) merger rates, which were computed in a similar manner as our BBH merger rates. Specifically, we calculate the neutron star density assuming a minimum stellar mass for producing a NS of $8~{}M_{\odot}$ (as we did in the CCS estimate in §3.3), and a maximum stellar mass equal to the minimum needed to form a black hole. The upper limit on the NS-NS binary merger rate density reported in Abbott et al. (2016g) is $\mathcal{R}_{NS}<12,600~{}\mathrm{Gpc^{-3}yr^{-1}}$. The $\epsilon=0.01$ case plotted is clearly well below this observational limit, which provides a weak constraint $\epsilon\lesssim 0.1$ for large $\tau$ and $\epsilon\lesssim 1$ for small $\tau$. We may further check our model using the Milky Way’s binary NS population and the short gamma ray burst (GRB) density. Kim et al. (2006) and Enrico Petrillo et al. (2013) estimate the binary NS merger rate should be $\mathcal{R}_{NS}\simeq 10^{2}-10^{3}~{}\mathrm{Gpc^{-3}yr^{-1}}$, which is consistent with the predictions shown in Figure 7 for $\tau\lesssim 5$ Gyr.

Having confirmed the consistency of our model with previous theoretical explorations and observational constraints, we now predict the merger rate density for even more massive compact objects – a regime that has not yet been probed observationally. Our expected rate density for black hole binary mergers each with $m_{\rm bh}>50M_{\odot}$ is $\mathcal{R}_{50}\gtrsim 1\,(\epsilon/0.01)~{}\mathrm{Gpc^{-3}yr^{-1}}$. With a rate density this high, we expect that a massive merger of this kind should be detected within the next decade. Mergers involving at least one black hole of this high mass should be more common.

4.1 Breaking degeneracies with host galaxy masses

One of the goals of gravitational wave astronomy is to constrain the astrophysics that underlies black hole merger detections, including 1) the physics of black hole binary formation and 2) the processes that drive subsequent mergers. We have parameterized these two global uncertainties using two simplifying parameters: the merger timescale $\tau$ and the binary black hole efficiency $\epsilon$. As demonstrated in Figure 6, current constraints on the merger rate provide only degenerate constraint on these parameters, and in particular allow a vast range of characteristic merger timescales, from fairly prompt mergers, $\tau\simeq 100$ Myr, to mergers that have taken a Hubble time to occur.

One possible way to break this degeneracy is to identify the host galaxy mass distribution for observed merger events. Small galaxies today have ongoing star formation, while larger galaxies tend to be quenched (e.g. Mannucci et al., 2010). Thus, binary mergers that occur soon after formation will more likely be seen in small galaxies. Mergers over timescales comparable to the age of the universe, however, will more closely track the overall stellar mass distribution. Most stars are in massive galaxies today (Baldry et al., 2012; Bernardi et al., 2013). Thus mergers detected locally that have take a long time to occur will be biased to reside within large galaxies.

An expanded network of gravitational wave detectors, including Advanced Virgo and the planned LIGO-India project, should be able to localize gravitational wave sources within a few square degrees (Nissanke et al., 2013; Abbott et al., 2016e). With enough detections, cross-correlating merger locations with galaxy counts on degree scales could enable constraints on the host mass for BBH mergers, as massive galaxies cluster more strongly with other galaxies than do lower mass sytems (e.g. Zehavi et al., 2012; Raccanelli et al., 2016). More precise determinations of host mass distributions would be enabled if there are electromagnetic counterparts to mergers. Unfortunately, BBH mergers are not expected to produce significant EM radiation except in extreme cases (Loeb, 2016, e.g.), though see Perna et al. (2016) for a more plausible scenario. On the other hand, NS-NS mergers are expected to produce short gamma-ray bursts (e.g. Narayan et al., 1992; Rosswog et al., 2003; Nakar & Piran, 2011). The Advanced LIGO/Virgo detector network should detect tens of NS-NS mergers per year (Abadie et al., 2010), which could enable a promising avenue for mapping out the host distributions for these mergers with some precision.

In order to provide some quantitative insight into how well host masses can help discriminate models with different timescales, we consider two extreme examples: one prompt ($\tau\sim 100$ Myr) and another comparable to the Hubble time ($\tau\sim 10$ Gyr). The rate of prompt mergers will be governed by the compact object birthrate in the low-z universe, $\dot{n}_{\rm co,0}\equiv\dot{n}_{\rm co}(t\simeq t_{0})$. Specifically, the merger rate for prompt mergers will be

[TABLE]

While Equation 5 can be used to provide a good estimate for $\dot{n}_{\rm co,0}$, it cannot be used to determine the host mass distribution for newly formed black holes or neutron stars. This is because Equation 5 is normalized by the total remnant density and contains no information on when and where those remnants were born – only where they are today.

In order to accurately determine the host mass distribution for compact objects being born recently, we can perform an analysis similar to the one we used in calculating $n_{\rm bh}$ (Equations 1 and 3). Starting with neutron star birthrates, we assume that they track massive star formation at a given gas-phase metallicity $Z_{g}$. Specifically, we need the minimum ($\mathcal{M}_{\rm min}^{\rm ns}=8M_{\odot}$) and maximum stellar progenitor mass that will produce a neutron star: $\mathcal{M}_{\rm max}^{\rm ns}(Z_{g})=\mathcal{M}_{\rm min}(m_{\rm co},Z)$ with $m_{\rm co}=5M_{\odot}$ is the same progenitor mass limit discussed in § 3.1. The observed mass-metallicity-SFR relation $\dot{M}_{\star}(M_{\star},Z_{g})$ (Mannucci et al., 2010; Lara-López et al., 2010) then provides a means to estimate the birth rate density by integrating over metallicity distribution function and stellar mass function:

[TABLE]

Note that the observed gas-phase metallicity relation provides results for the oxygen abundance and we are interested in the overall metallicity. We account for this following Peeples & Somerville (2013) and use the $[\alpha/\mathrm{Fe}]$-mass relation from Thomas et al. (2005) to determine $Z_{g}$ from the gas-phase oxygen abundance.

The birthrate calculation for black holes is the same as the above calculation for neutron stars except for the limits of the integral over the IMF. For black holes of mass $>m_{\rm bh}$ the lower limit is $\mathcal{M}_{\rm min}(m_{\rm bh},Z)$ and the upper limit is $\mathcal{M}_{\rm u}$ (as defined for Equation 1).

We emphasize that both Equations 5 and 9 give almost identical answers for the global compact object birthrate rate today when Equation 5 is evaluated at $t=t_{0}$. This must be the case if $\psi(t)$ is normalized self-consistently. However, Equation 9 now allows us to determine the host mass distribution for newly formed objects. That is, we can differentiate Equation 8 with respect to galaxy stellar mass to derive the black hole merger rate density per host mass in the case that black hole mergers occur promptly after black hole binary formation.

Figure 8 shows the results of these calculations for binary NS-NS mergers in cyan, all BBH mergers in black, and BBH mergers with $\mathrm{m_{bh}>30~{}M_{\odot}}$ in yellow. Dashed lines indicate prompt mergers of compact objects formed recently ($\tau\sim 100$ Myr), derived by differentiating Equation 8. Solid lines indicate long-timescale mergers of objects formed at redshift $z=2$ ($\tau\sim 10$ Gyr), derived by differentiating Equation 5 (and thus Equation 3) with respect to stellar mass. For this figure, we have used merging efficiencies of $\epsilon=0.007$ for prompt mergers and $\epsilon=0.001$ for delayed mergers in order that they both produce the same global rates, specifically the mean reported LIGO rate for the overall and massive black hole populations.

The median host galaxy mass for prompt NS-NS mergers is $M_{\star}=1.4\times 10^{9}~{}M_{\odot}$ while for delayed mergers it is $M_{\star}=3.9\times 10^{10}~{}M_{\odot}$. The main difference between the distributions is in the low-mass tail, where prompt BBH merger rate density varies slowly with host mass as $M_{\star}^{1/7}$ while long timescale BBH merger rate density increases more sharply with host mass as $M_{\star}^{1/2}$. This is caused by the shape of the mass-metallicity-SFR relation; while the total number of black holes formed in a galaxy is largely independent of its mean stellar metallicity (2) below a certain threshold, the number forming locally depends on the specific star formation rates of the host galaxies, which changes only moderately with stellar mass at these scales (Leitner, 2012; Tomczak et al., 2016, e.g.). In the prompt scenario, this implies 47% of NS-NS mergers occur in dwarf galaxies with $M_{\star}<10^{9}~{}M_{\odot}$ while in the delayed-merging scenario only $6\%$ are hosted by dwarfs. Given the expected binary NS-NS detection rate for the Advanced LIGO/VIRGO network (Abadie et al., 2010) and the likelihood for electromagnetic counterpart signals, it should be possible to map out the host galaxy population in coming decade and to explore the question of whether these mergers have occurred long after formation (thus tracing stellar mass and the most massive galaxies) or promptly after formation (tracing star formation and dwarf galaxies). The LIGO India detector will only enhance the ability to perform this experiment.

The dependence of host mass distribution on merger timescale also exists for BBH mergers. For all $>5M_{\odot}$ BBH mergers, the median host galaxy masses are $9.5\times 10^{8}~{}M_{\odot}$ and $2.6\times 10^{10}~{}M_{\odot}$ for the short and long merger timescales, respectively. The fraction of mergers hosted by $M_{\star}<10^{9}~{}M_{\odot}$ dwarf galaxies is 51% and $9\%$ in the two cases. For the more massive $>30M_{\odot}$ “GW150914-like” mergers, 95% occur in dwarfs for the prompt case while just 24% are hosted by dwarfs in the long-timescale case.

To assess the feasibility of discriminating between the prompt and the long-timescale scenarios for BBH mergers, we can draw an analogy with the ultra-high energy cosmic rays. The origin of the ultra-high energy cosmic rays (UHECRs) with measured energies in excess of $10^{19}$ GeV is a mystery. The strong energy losses at these energies due to interactions with the Cosmic Microwave Background imply that the sources must lie within a Gpc. At energies above $5\times 10^{19}$ GeV, the abrupt (Greisen-Zatsepin-Kuzmin or GZK) cutoff due to the $\Delta^{+}$ resonance leads to a dramatic decrease in the implied distance to the sources. The angular resolution of the measurements and the smearing due to bending in magnetic fields imply that the actual direction of the UHECR events can be reconstructed to only about 5 degrees, a situation similar to that for gravitational wave detections.

The BBH mergers are detectable to distances of a Gpc and unlike the case with the UHECRs, we expect some constraint on the source distance from the gravitational wave detections. For the UHECR case, it seems like the 69 events from the Pierre Auger Observatory is a large enough set of events to test for correlation with putative local sources (Abreu et al., 2010). No clear consensus has been reached on the origin of these UHECR events. In particular, only a fraction of the UHECRs can be attributed to the local catalogs of AGNs or 2MASS galaxies that have been cross correlated (Khanin & Mortlock, 2016). While the origin of UHECR has not been clarified by these analyses, they have demonstrated the feasibility of testing various hypotheses for the origin of UHECRs with upwards of about 50 events.

If the origin of the BBH mergers is related to stellar mass remnants, then our analogy with UHECR correlation studies suggests that looking for correlation large, clustered galaxies could be a fruitful way forward. Because the galaxy clustering bias begins to increase around masses of $M_{\star}\gtrsim 10^{10.7}$ at $z<0.3$ (e.g. Zu & Mandelbaum, 2015) we use this as our mass threshold. For the long timescale scenario, our prediction is that 30% of BBH mergers with masses $>5M_{\odot}$ and 42% of NS mergers occur in $M_{\star}\gtrsim 10^{10.7}$ hosts, while the corresponding numbers for the prompt scenario are 5.4% and 4.6%. Thus, we would expect a cross correlation with these highly biased galaxies only in the long timescale scenario (the situation for the most massive $>30M_{\odot}$ BBH mergers is even more stark, though less observable: 6.4% and $1.0\times 10^{-3}\%$, respectively). Testing a model that includes templates for both the red and dwarf galaxies, it may be possible to infer the relative contributions of each to the observed mergers. Given our correlated predictions for the NS-NS mergers and their possibile EM counterparts, nailing down the timescale for the mergers seems a likely possibility.

Note that in calculating our host mass distributions for long-timescale mergers we have assumed that black hole pairs that formed near the peak in cosmic star formation rate density ($z=2$) have distributed themselves like the bulk of the black holes in the Universe today. This assumption is conservative in the sense that it biases mergers to occur in lower mass hosts than they otherwise would. In reality, black holes that formed at $z=2$ will reside in slightly higher mass galaxies than the bulk of the black hole population (since half the stars formed after this time, and later star formation occurs in smaller galaxies). Given this, it is possible that the dichotomy in host mass populations for prompt and long merger timescale populations is sharper than that in Figure 8.

5 Conclusions

In this paper we have worked through empirically-derived expectations for the stellar remnant black hole population in the Universe and used these as a basis for interpreting gravitational wave signals such as those being detected by Advanced LIGO and eventually Advanced Virgo (LIGO Scientific Collaboration et al., 2015; Acernese et al., 2015). We have quantified our uncertainties using two parameters: the binary black hole efficiency $\epsilon$ (Equation 4), and the typical merger timescale $\tau$ (Equations 6 and 7).

Stellar-remnant black holes should be abundant in the local universe. For example, $m_{\rm bh}>30M_{\odot}$ black holes should have a local number density of $0.9-2\times 10^{14}\ \mathrm{Gpc^{-3}}$ with the range reflecting variations between our fiducial Kroupa (2002) and metallicity dependent (Geha et al., 2013) IMFs. This corresponds to an occupation rate of $\sim 1$ per 1000 $M_{\odot}$ of stars formed in galaxies with $M_{\star}\lesssim 10^{10}M_{\odot}$ (see Figures 2 and 3). Such an abundant black hole population provides an ample source for binary systems that eventually merge for reasonable choices of parameters that characterize the merger process.

If $\epsilon\simeq 1\%$ of stellar remnant black holes end up in a binary configuration that eventually merges, then the current LIGO merger rate constraints can be accommodated as long as the typical merger timescale is $\tau\lesssim 5$ Gyr (See Figure 6). If mergers tend to occur over a timescale that coincides with the peak in cosmic star formation rate density ($\tau\simeq 10$ Gyr) then the efficiency of binary mergers would need to be smaller ($\epsilon\simeq 0.1\%$) in order to be consistent with the observed rates.

Though our approach is not well suited for ab initio calculations, it does provide fairly robust scalings because the uncertain/unknown parameters are reasonably constant for all compact objects in our calculations. For example, for any $\epsilon$ or $\tau$, $50~{}M_{\odot}$ black holes should have merger rate densities that are a factor of $7\pm 1$ smaller than merger rates of binary $30~{}M_{\odot}$ black holes (see Fig. 7). This range accounts for uncertainties in the faint end of the stellar mass function (taken from Geller et al., 2012; Baldry et al., 2012; Lan et al., 2016). Scaling from the event-based rate derived for GW150914, we would therefore predict the rate for $50~{}M_{\odot}$ black holes binary mergers to be $\mathcal{R}_{50}=8^{+27}_{-6}~{}\mathrm{Gpc^{−3}yr^{-1}}$. This places $50~{}M_{\odot}$ black hole binary mergers at the limit where detection by LIGO within the next decade should be expected. In principle, the mass specturm of detected compact objects will provide information on the galaxy SMF, as few massive detections would imply a shallow faint-end slope of $\alpha\sim-1.3$, while a large number would support slopes closer to$\alpha\sim-1.7$.

Given the degeneracy between merger timescale and binary efficiency in producing the observed range of merger rates, we have explored one possible avenue for breaking this degeneracy. In Figure 8, we showed that for very short timescale “prompt” mergers, which occur soon after black hole formation, the host galaxy population is expected to track the local star formation, and therefore be skewed towards smaller galaxies. For example, about half of the BBH mergers with $m_{\rm bh}>5M_{\odot}$ should occur in hosts with stellar masses $<10^{9}M_{\odot}$ in the prompt merger scenario, while only 10% of such events should be hosted by these dwarf galaxies in the long timescale scenario.

As we move towards an era where a global network of gravitational wave detectors is likely, we can expect source localization to provide a means towards discriminating scenarios of this kind. For BBH mergers with no electromagnetic counterparts, the host distribution may in principle be inferred by searching for correlation of these events with background galaxy population, in much the same way as has been attemped with ultra-high energy cosmic rays measured by the Pierre Auger Observatory (§4.1). For NS-NS mergers, we expect electromagnetic counterparts, which would be more direct way to determine the host masses. If we are able to map out the source distribution for NS-NS and BBH mergers directly through electromagnetic counterparts or indirectly through the anisotropy of source distribution, then we will be able to constrain formation and evolution scenarios for binary black hole and neutron star merger events.

Acknowledgments

Support for this work was provided by NASA through Hubble Space Telescope grants HST-GO-12966.003-A and HST-GO-13343.009-A. We thank A. Lamberts, S. Garrison-Kimmel and our referee for useful discussions.