Where Are LIGO's Big Black Holes?

One of the most fundamental quantities in gravitational-wave (GW) astrophysics is the mass distribution of stellar-mass black holes (BHs). Characterizing this distribution in merging binary systems is crucial to understanding stellar evolution, supernova physics, and the formation of compact binary systems. Prior to the first gravitational-wave detections of binary black holes (BBHs), the sample of ∼20 BHs in X-ray binary systems was used to infer the BH mass distribution (Özel et al. 2010; Farr et al. 2011), providing strong evidence for the existence of a mass gap between the heaviest neutron star (NS; $\sim 2\mbox{--}3\,{M}_{\odot }$) and the lightest BH ($\sim 4\mbox{--}5\,{M}_{\odot }$; but see also Kreidberg et al. 2012). The presence of a mass gap between NSs and BHs has critical implications for supernova explosion theory (Belczynski et al. 2012), and there are several proposed methods to probe this mass gap with GW observations of compact binaries (Littenberg et al. 2015; Kovetz et al. 2017; Mandel et al. 2017). In addition to the low-mass gap, supernova theory suggests that pulsational pair-instability supernovae (PPISN; Heger & Woosley 2002) and pair-instability supernovae (PISN; Fowler & Hoyle 1964; Rakavy et al. 1967; Bond et al. 1984) lead to a second mass gap between ∼50 and $135\,{M}_{\odot }$ for BHs formed from stellar core collapse (Belczynski et al. 2014, 2016a; Marchant et al. 2016; Spera & Mapelli 2017; Woosley 2017). Although studies of the lower-mass gap have to wait for many more binary detections, because LIGO's sensitivity is almost 500 times greater for $50\mbox{--}50\,{M}_{\odot }$ mergers than $3\mbox{--}3\,{M}_{\odot }$ mergers, the existing data already begin to probe the upper mass gap.

For the first four BBH detections (GW150914, LVT151012, GW151226, and GW170104), the LIGO-Virgo Collaboration fit the BBH mass distribution with a power-law parameterization on the primary BBH mass, m₁ (Abbott et al. 2016a, 2017b), inspired by the stellar initial mass function (IMF) (Salpeter 1955; Kroupa 2001). Specifically, Abbott et al. (2016a) use the following one-parameter power law to model the distribution of primary component BH masses:

$$\begin{eqnarray}&&p({m}_{1}| \alpha )\propto {m}_{1}^{-\alpha },\end{eqnarray} \tag{ 1 }$$

where ${M}_{\min }\lt {m}_{1}\lt {M}_{\max }$. The mass ratio between component BHs, $q\equiv {m}_{2}/{m}_{1}\leqslant 1$, is assumed to be uniformly distributed in the allowed range ${M}_{\min }/{m}_{1}\leqslant q\leqslant \min ({M}_{\mathrm{tot},\max }/{m}_{1}-1,1)$. Thus, the marginal distribution of the secondary component mass, m₂, is given by

$$\begin{eqnarray}&&p({m}_{2}| {m}_{1})\propto \displaystyle \frac{1}{\min ({m}_{1},{M}_{\mathrm{tot},\max }-{m}_{1})-{M}_{\min }},\end{eqnarray} \tag{ 2 }$$

and therefore the joint mass distribution is

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}| \alpha )\propto \displaystyle \frac{{m}_{1}^{-\alpha }}{\min ({m}_{1},{M}_{\mathrm{tot},\max }-{m}_{1})-{M}_{\min }}.\end{eqnarray} \tag{ 3 }$$

The only free parameter in this assumed mass distribution is the power-law slope, α. The minimum BH mass, M_min, is fixed at ${M}_{\min }=5\,{M}_{\odot }$, and the maximum mass, M_max, is fixed at ${M}_{\max }=100\,{M}_{\odot }-{M}_{\min }$. Meanwhile, the total BBH mass, ${M}_{\mathrm{tot}}={m}_{1}+{m}_{2}$, is also restricted: ${M}_{\mathrm{tot}}\leqslant {M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$. (Enforcing ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$ causes a break in the power law at $50\,{M}_{\odot }$.) The choice for M_min is motivated by the empirical lower-mass gap, while M_max and ${M}_{\mathrm{tot},\max }$ are set by the stellar binary matched-filter search, which defines stellar-mass BBHs as those with source-frame total masses ${m}_{1}+{m}_{2}\leqslant 100\,{M}_{\odot }$ (Cannon et al. 2012; Abbott et al. 2016b; Usman et al. 2016). However, LIGO is in principle sensitive to heavier BBHs (Abbott et al. 2017e). BBHs with detector-frame total masses up to $600\,{M}_{\odot }$ can be detected via matched-filtering by the intermediate-mass black hole (IMBH) modeled search (Dal Canton & Harry 2017; Messick et al. 2017; Nitz et al. 2017), and IMBHs of even higher mass can be detected by the unmodeled transient search (Klimenko et al. 2008; Abbott et al. 2017a).

We observe that a key assumption of the distribution in Equation (3) is that BHs in merging binaries follow the same mass distribution from $5\,{M}_{\odot }$ to at least $50\,{M}_{\odot }$ and that there exist BHs as heavy as $95\,{M}_{\odot }$. Meanwhile, LIGO is extremely sensitive to heavy BBHs. The first-order post-Newtonian approximation predicts that for low-mass BBHs and a Euclidean universe, the spacetime volume, VT, to which LIGO can detect a BBH merger of a fixed mass ratio increases with its primary component mass, m₁, roughly as ${VT}\propto {m}_{1}^{5/2}$. In the following section, we find that when accounting for cosmology and taking BBHs over the entire range $10\,{M}_{\odot }\lt {M}_{\mathrm{tot}}\lt 100\,{M}_{\odot }$, it is still a good approximation to take ${VT}\propto {m}_{1}^{k}$, with $k\sim 2.2$. This means that if the BBH mass distribution follows a power law with slope α as in Equation (3), we expect the mass distribution among detected BHs to follow ${m}_{1}^{-\alpha +2.2}$. For a Salpeter IMF ($\alpha =2.35$), this implies an almost flat detected distribution of BBH masses. Thus, the absence of heavy BBHs in the data quickly indicates either that the mass distribution declines steeply toward high masses ($\alpha \gg 2.2$) or that an upper mass gap sharply cuts off the mass distribution.

In this Letter, we show that we can start to distinguish between these two scenarios with the first four LIGO BBH detections (including LVT151012, which has an 87% probability of being astrophysical; Abbott et al. 2016a, 2017b). Using simulated BBH detections, we demonstrate that if an existing mass gap is not accounted for, the non-detection of heavy stellar-mass BHs will quickly bias the power-law fit to distributions that are erroneously too steep. However, by including a maximum BH mass as a free parameter in the analysis, we can simultaneously infer the shape of the mass distribution and the location of a mass gap, if present. We carry out this analysis for the first four BBHs as well as for simulated BBH detections. We find that there is already evidence for an upper mass cutoff at $\sim 40\,{M}_{\odot }$ from the first four detections, a conclusion that is further supported by the two recently reported BBH detections (GW170608 and GW170814; Abbott et al. 2017d, 2017c). We show that with ${ \mathcal O }(10)$ additional detections, the presence and location of the bottom edge of the mass gap will be highly constrained.

As we noted in the previous section, for a given mass ratio, the sensitivity of the LIGO-Virgo search scales with primary component mass roughly as ${{m}_{1}}^{2.2}$. We characterize the sensitivity by the redshifted spacetime volume, VT, for which a given search is sensitive to a BBH system of given masses. If we assume that the rate of BBH coalescences is uniform in comoving volume and source-frame time and neglect the effects of BH spin on the detectability of a source, VT depends only on the power spectral density (PSD) curve characterizing the detectors and the BBH component masses (Abbott et al. 2016d). Under these assumptions, VT is given by

$$\begin{eqnarray}&&{VT}({m}_{1},{m}_{2})=T\int {dz}\displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{1}{1+z}f(z,{m}_{1},{m}_{2}),\end{eqnarray} \tag{ 4 }$$

where T is the search time, V_c is the comoving volume, and $0\lt f(z,{m}_{1},{m}_{2})\lt 1$ is the probability that a BBH system of masses m₁, m₂ at redshift z will be detected. We adopt the cosmological parameters from Ade et al. (2016) throughout the calculation. To calculate the detection probability, $f(z,{m}_{1},{m}_{2})$, we use the semi-analytic approximation from Abbott et al. (2016d). Taking the PSD function corresponding to the early aLIGO high-sensitivity scenario in Abbott et al. (2016c; a good approximation to the PSD during the first and second aLIGO observing runs), we calculate the optimal matched-filter signal-to-noise ratio (S/N) of a BBH with component masses m₁ and m₂ and zero spins located at redshift z. The optimal S/N, ${\rho }_{\mathrm{opt}}$, corresponds to a face-on source that is directly overhead to a single detector. We then generate random angular factors, $0\lt w\lt 1$, from a single-detector antenna power pattern, assuming that binaries are distributed uniformly on the sky with isotropic inclination vectors (Finn & Chernoff 1993; Dominik et al. 2015; Belczynski et al. 2016c). The angular factor, w, characterizes the response of a detector to a source at a given sky location and orientation (so that w = 1 for an overhead, face-on source). From the distribution of angular factors, w, we assign a distribution of S/Ns, $\rho =w{\rho }_{\mathrm{opt}}$, for each source with parameters $({m}_{1},{m}_{2},z)$. Out of this distribution of S/Ns, the fraction that exceeds the single-detector threshold, $\rho \gt 8$, roughly corresponding to a network threshold, $\rho \gt 12$, is taken to be the detection probability $f(z,{m}_{1},{m}_{2})$. This semi-analytic calculation for $f(z,{m}_{1},{m}_{2})$ neglects the effects of non-Gaussian noise, which tends to raise the S/N detection threshold for binaries of very high mass. However, it remains a good approximation for stellar-mass binaries of total masses up to at least ${M}_{\mathrm{tot}}=100\,{M}_{\odot }$ and possibly higher (Abbott et al. 2016d, 2017e).

The expected sensitive redshifted spacetime volume as a function of total mass is shown in Figure 1, calculated for one year of observation (T = 1 year) at the O1–O2 LIGO sensitivity. For example, we note that in O1 and O2 the LIGO detectors probed a volume roughly seven times larger for $75\mbox{--}75\,{M}_{\odot }$ binaries as compared to $25\mbox{--}25\,{M}_{\odot }$ binaries. In particular, since ${m}_{1}={M}_{\mathrm{tot}}/(1+q)$ for a fixed mass ratio q, we can see from Figure 1 that VT scales approximately as m₁^k with $k\sim 2.2$ for ${M}_{\mathrm{tot}}\lesssim 100\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

To calculate the sensitivity to a population of BBHs, the relevant quantity is the population-averaged spacetime volume, $\langle {VT}\rangle$. If we know the distribution of masses across the population of BBHs, ${p}_{\mathrm{pop}}({m}_{1},{m}_{2})$, assuming negligible spins and a constant comoving merger rate, we can calculate the population-averaged sensitive spacetime volume (Equation (15) in Abbott et al. 2016d):

$$\begin{eqnarray}&&\langle {VT}\rangle =\int \int {VT}({m}_{1},{m}_{2}){p}_{\mathrm{pop}}({m}_{1},{m}_{2}){{dm}}_{1}{{dm}}_{2},\end{eqnarray} \tag{ 5 }$$

where the first integral is over ${M}_{\min }\lt {m}_{1}\lt {M}_{\max }$ and the second integral is over ${M}_{\min }\lt {m}_{2}\lt \min ({m}_{1},{M}_{\max ,\mathrm{tot}}-{m}_{1})$. $\langle {VT}\rangle$ relates the specific merger rate, R, to the expected number, Λ, of BBH signals in a given detection period (Abbott et al. 2016d):

$$\begin{eqnarray}&&{\rm{\Lambda }}=R\langle {VT}\rangle .\end{eqnarray} \tag{ 6 }$$

The number of BBH detections, n, follows a Poisson process with mean Λ. To explore the existence of a high-mass gap in Section 4.1, we compare the expected number of low-mass BBH signals, ${{\rm{\Lambda }}}_{\mathrm{low}}=R\langle {VT}{\rangle }_{\mathrm{low}}$, to the expected number of high-mass BBH signals, ${{\rm{\Lambda }}}_{\mathrm{high}}=R\langle {VT}{\rangle }_{\mathrm{high}}$, for different power-law populations, where low (high) mass is defined by the primary component mass ${m}_{1}\leqslant {M}_{\mathrm{cutoff}}$ (${m}_{1}\gt {M}_{\mathrm{cutoff}}$). We define

$$\begin{eqnarray}&&\displaystyle \frac{1}{r}\equiv \displaystyle \frac{{{\rm{\Lambda }}}_{\mathrm{low}}}{{{\rm{\Lambda }}}_{\mathrm{high}}}=\displaystyle \frac{\langle {VT}{\rangle }_{\mathrm{low}}}{\langle {VT}{\rangle }_{\mathrm{high}}},\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}\langle {VT}{\rangle }_{\mathrm{low}} & \equiv & \displaystyle \int {\displaystyle \int }_{{M}_{\min }}^{{M}_{\mathrm{cutoff}}}{VT}({m}_{1},{m}_{2}){p}_{\mathrm{pop}}({m}_{1},{m}_{2}){{dm}}_{1}{{dm}}_{2}\\ \langle {VT}{\rangle }_{\mathrm{high}} & \equiv & \displaystyle \int {\displaystyle \int }_{{M}_{\mathrm{cutoff}}}^{{M}_{\max }}{VT}({m}_{1},{m}_{2}){p}_{\mathrm{pop}}({m}_{1},{m}_{2}){{dm}}_{1}{{dm}}_{2}.\end{eqnarray} \tag{ 8 }$$

The integration limits on the m₂ integral in Equation (8) are identical to those in Equation (5) so that the total $\langle {VT}\rangle =\langle {VT}{\rangle }_{\mathrm{low}}+\langle {VT}{\rangle }_{\mathrm{high}}$. We can then compute the probability of detecting n_high BBHs with primary component mass ${m}_{1}\gt {M}_{\mathrm{cutoff}}$, given that we have detected n_low BBHs with primary component mass ${m}_{1}\lt {M}_{\mathrm{cutoff}}$. (We ignore mass measurement uncertainties that may prevent us from definitively assigning a BBH to either the low- or high-mass class.) This probability is given by

$$\begin{eqnarray}&&p({n}_{\mathrm{high}}| {n}_{\mathrm{low}})={\displaystyle \int }_{0}^{\infty }{\displaystyle \int }_{0}^{\infty }p({n}_{\mathrm{high}},{{\rm{\Lambda }}}_{\mathrm{high}},{{\rm{\Lambda }}}_{\mathrm{low}}| {n}_{\mathrm{low}})\\ &&\quad d{{\rm{\Lambda }}}_{\mathrm{high}}d{{\rm{\Lambda }}}_{\mathrm{low}}\\ &&\quad =\displaystyle \int \displaystyle \int p({n}_{\mathrm{high}}| {{\rm{\Lambda }}}_{\mathrm{high}})p({{\rm{\Lambda }}}_{\mathrm{high}}| {{\rm{\Lambda }}}_{\mathrm{low}})p\\ &&\qquad \times \,({{\rm{\Lambda }}}_{\mathrm{low}}| {n}_{\mathrm{low}})d{{\rm{\Lambda }}}_{\mathrm{high}}d{{\rm{\Lambda }}}_{\mathrm{low}}\\ &&\quad =\displaystyle \int \displaystyle \int p({n}_{\mathrm{high}}| {{\rm{\Lambda }}}_{\mathrm{high}})\delta ({{\rm{\Lambda }}}_{\mathrm{high}}-r{{\rm{\Lambda }}}_{\mathrm{low}})p\\ &&\qquad \times \,({{\rm{\Lambda }}}_{\mathrm{low}}| {n}_{\mathrm{low}})d{{\rm{\Lambda }}}_{\mathrm{high}}d{{\rm{\Lambda }}}_{\mathrm{low}}\\ &&\quad =\displaystyle \int p({n}_{\mathrm{high}}| r{{\rm{\Lambda }}}_{\mathrm{low}})p({{\rm{\Lambda }}}_{\mathrm{low}}| {n}_{\mathrm{low}})d{{\rm{\Lambda }}}_{\mathrm{low}}\\ &&\qquad \propto \displaystyle \int p({n}_{\mathrm{high}}| r{{\rm{\Lambda }}}_{\mathrm{low}})p({n}_{\mathrm{low}}| {{\rm{\Lambda }}}_{\mathrm{low}}){p}_{0}({{\rm{\Lambda }}}_{\mathrm{low}})d{{\rm{\Lambda }}}_{\mathrm{low}},\end{eqnarray} \tag{ 9 }$$

where in the third line we used the definition of r given by Equation (8) and in the last line we used Bayes's theorem. In Equation (9), terms like $p(n| {\rm{\Lambda }})$ denote the Poisson probability of n with mean Λ. We take the prior ${p}_{0}({{\rm{\Lambda }}}_{\mathrm{low}})$ to be the Jeffrey's prior:

$$\begin{eqnarray}&&{p}_{0}({{\rm{\Lambda }}}_{\mathrm{low}})\propto \displaystyle \frac{1}{\sqrt{{{\rm{\Lambda }}}_{\mathrm{low}}}}.\end{eqnarray} \tag{ 10 }$$

We will return to Equation (9) in Section 4.1.

Our goal is to jointly infer the shape of the BBH mass distribution along with the lower edge of a potential mass gap, M_max. We therefore follow Abbott et al. (2016a) in fitting a power-law mass distribution to GW BBH mass measurements, but we add the maximum BH mass, M_max, as a free parameter. We leave the minimum BH mass, M_min, fixed at ${M}_{\min }=5\,{M}_{\odot }$. Thus, we consider a two-parameter mass distribution:

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}| \alpha ,{M}_{\max })\\ &&\quad \propto \,\displaystyle \frac{{m}_{1}^{-\alpha }\ { \mathcal H }({M}_{\max }-{m}_{1})}{\min ({m}_{1},{M}_{\mathrm{tot},\max }-{m}_{1})-{M}_{\min }},\end{eqnarray} \tag{ 11 }$$

where ${ \mathcal H }$ is the Heaviside step function that enforces a cutoff in the distribution at ${m}_{1}={M}_{\max }$. For consistency with the LIGO definition of a stellar-mass BH, we restrict ${M}_{\max }\leqslant 100\,{M}_{\odot }$ throughout. Furthermore, as in the LIGO Collaboration's analysis, we enforce ${m}_{1}+{m}_{2}\leqslant {M}_{\mathrm{tot},\max }$, which provides an additional constraint for ${M}_{\max }\gt \tfrac{1}{2}{M}_{\mathrm{tot},\max }$. The LIGO Collaboration fixes ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$ and ${M}_{\max }=100\,{M}_{\odot }-{M}_{\min }$, as this is the definition of a stellar-mass BBH set by the search. This choice corresponds to one of the following assumptions: either (a) BBHs with total source-frame masses ${M}_{\mathrm{tot}}\gt 100\,{M}_{\odot }$ do not exist as part of the population of stellar-mass BBHs or (b) LIGO is not sensitive to BBHs with total source-frame masses ${M}_{\mathrm{tot}}\gt 100\,{M}_{\odot }$, so we cannot constrain their existence. (In the absence of detections with ${M}_{\mathrm{tot},\max }\gt 100\,{M}_{\odot }$, setting ${M}_{\mathrm{tot},\max }\leqslant 100\,{M}_{\odot }$ in the population model, Equation (11), is equivalent to assuming that the sensitivity vanishes for binaries with ${M}_{\mathrm{tot}}\gt 100\,{M}_{\odot }$.) Assumption (a) may not be well motivated, as population-synthesis models that predict stellar BBHs with component masses ${M}_{\max }\sim 50\mbox{--}100\,{M}_{\odot }$ tend to allow ${M}_{\mathrm{tot},\max }\sim 2{M}_{\max }$ (Belczynski et al. 2016b; Eldridge & Stanway 2016). Assumption (b) can also be questioned, as LIGO is sensitive to BBHs with detector-frame total masses up to at least $600\,{M}_{\odot }$ in the IMBH matched-filter search (Abbott et al. 2017e). However, the sensitivity may be lower than expected for very high mass BBHs due to non-stationary instrumental noise (Slutsky et al. 2010) and the absence of precessing and higher-order mode template waveforms, which leads to worse matches between signal and template for very high mass BBHs in the matched-filter search (Capano et al. 2014; Calderón Bustillo et al. 2017a, 2017b). If we had an accurate model of ${VT}({m}_{1},{m}_{2})$ across the mass range $5\lt {m}_{1},{m}_{2}\lt 100\,{M}_{\odot }$ (by performing a large-scale injection campaign), we could set ${M}_{\mathrm{tot},\max }=2{M}_{\max }\leqslant 200\,{M}_{\odot }$ in Equation (11). However, because our calculation of VT may be overestimating the sensitivity to binaries with ${M}_{\mathrm{tot},\max }\gt 100\,{M}_{\odot }$, when fitting Equation (11) in the following sections, we repeat the analysis once under the assumption that ${M}_{\mathrm{tot},\max }=\min (2{M}_{\max },100\,{M}_{\odot })$ and once assuming that ${M}_{\mathrm{tot},\max }=2{M}_{\max }\leqslant 200\,{M}_{\odot }$.

To extract the parameters of our assumed mass distribution (Equation (11)) from data, we use the same hierarchical Bayesian methods as Appendix D of Abbott et al. (2016a), further explained in Mandel et al. (2016). While GW data are noisy and subject to selection effects, both the measurement uncertainties and selection effects are well quantified. The selection effects refer to the mass-dependent detection efficiency. Under the assumptions of negligible BH spins and a uniform comoving merger rate, the detection efficiency is proportional to the sensitive spacetime volume ${VT}({m}_{1},{m}_{2})$ as described in Section 2 and Abbott et al. (2016d).

BBH masses are measured using the LALInference parameter-estimation pipeline, which calculates the posterior probability density function (PDF) of all parameters that govern the waveform given the data, d_i, from a BBH detection (Veitch et al. 2015). For an individual system, measurements of m₁ and m₂ take the form of ${ \mathcal O }$(10,000) posterior samples drawn from the posterior PDF, $p({m}_{1},{m}_{2}| {d}_{i})$. In the following section, we perform our analysis on published mass measurements from the first four BBHs as well as on simulated BBH measurements. We use the fact that the one-dimensional PDFs for the source-frame chirp mass and symmetric mass ratio are well approximated by independent (uncorrelated) Gaussian distributions.

For the first four BBH sources, GW150914, LVT151012, GW151226, and GW170104, we approximate the source-frame chirp mass posterior PDF as a Gaussian with a mean and standard deviation given by the median and 90% credible intervals listed in Table 4 of Abbott et al. (2016a) or Table 1 of Abbott et al. (2017b). In the case that the 90% credible interval is slightly asymmetric about the median, we use the average to estimate the standard deviation. We likewise approximate the posterior PDF of the symmetric mass ratio, $\eta =q/{(1+q)}^{2}$, as a Gaussian truncated to the allowed range $[0,0.25]$, with a mean and standard deviation given by the entry for q in the same tables. Using these approximate chirp mass and symmetric mass ratio distributions, we generate 25,000 posterior samples from the component mass posterior PDFs of each event.

For our set of simulated BBH detections, we generate a set of component masses from an underlying mass distribution. To each BBH system, we assign a redshift from a redshift distribution that is uniform in the merger-frame comoving volume. Given the simulated masses and redshift of each BBH, we randomly generate its single-detector S/N from the antenna power pattern, using the PSD corresponding to the early aLIGO high-sensitivity scenario (as described in Section 2). Out of this population, the set of detections are those simulated BBHs with a single-detector S/N satisfying $\rho \gt 8$. Given the true component masses and the S/N of each mock BBH detection, we produce realistic mass measurements by generating 5000–10,000 posterior samples for the component masses following the prescription in Equation (1) of Mandel et al. (2017). Given true values for the chirp mass, ${{ \mathcal M }}^{T}$, symmetric mass ratio, ${\eta }^{T}$, and S/N, ${\rho }^{T}$, we draw chirp mass posterior samples from a Gaussian distribution centered at $\bar{{ \mathcal M }}$ with standard deviation ${\sigma }_{{ \mathcal M }}$ and symmetric mass ratio posterior samples from a Gaussian distribution centered at $\bar{\eta }$ with standard deviation ${\sigma }_{\eta }$. We only keep posterior samples with $0.01\lesssim \eta \lesssim 0.25$. The variables $\bar{{ \mathcal M }}$ and $\bar{\eta }$ are drawn from Gaussian distributions:

$$\begin{eqnarray}\bar{{ \mathcal M }} & \sim & N({{ \mathcal M }}^{T},{\sigma }_{{ \mathcal M }}),\\ \bar{\eta } & \sim & N({\eta }^{T},{\sigma }_{\eta }),\end{eqnarray} \tag{ 12 }$$

where ${\sigma }_{{ \mathcal M }}$, ${\sigma }_{\eta }$ scale inversely with the S/N and are given in Mandel et al. (2017).

Once we have samples from the posterior PDF, $p({m}_{1},{m}_{2}| {d}_{i})$, for each event (both real and simulated) and we have calculated the detection efficiency, ${P}_{\det }({m}_{1},{m}_{2})\propto {VT}({m}_{1},{m}_{2})$, we follow Appendix D in Abbott et al. (2016a) to fit Equation (11). The likelihood for a single BBH detection given the parameters of the mass distribution, α and M_max, is given by

$$\begin{eqnarray}p({d}_{i}| \alpha ,{M}_{\max }) & \propto & \displaystyle \frac{\displaystyle \int p({d}_{i}| {m}_{1},{m}_{2})p({m}_{1},{m}_{2}| \alpha ,{M}_{\max }){{dm}}_{1}{{dm}}_{2}}{\beta (\alpha ,{M}_{\max })}\\ & \propto & \displaystyle \frac{\langle p({m}_{1},{m}_{2}| \alpha ,{M}_{\max })\rangle }{\beta (\alpha ,{M}_{\max })},\end{eqnarray} \tag{ 13 }$$

where $\langle \ldots \rangle$ denotes an average over the $({m}_{1},{m}_{2})$ posterior samples. This is valid because for each event, $p({d}_{i}| {m}_{1},{m}_{2})\,\propto p({m}_{1},{m}_{2}| {d}_{i})$, as the prior on ${m}_{1},{m}_{2}$ is taken to be flat. Therefore, we can calculate the integral in the first line of Equation (13) by taking the average of $p({m}_{1},{m}_{2}| \alpha ,{M}_{\max })$ over the mass posterior samples. Meanwhile, $\beta (\alpha ,{M}_{\max })$ is defined as

$$\begin{eqnarray}&&\beta (\alpha ,{M}_{\max })\equiv \int p({m}_{1},{m}_{2}| \alpha ,{M}_{\max }){VT}({m}_{1},{m}_{2}){{dm}}_{1}{{dm}}_{2}.\end{eqnarray} \tag{ 14 }$$

The likelihood for the data across all events ${\boldsymbol{d}}=\{{d}_{i}\}$ is the product of the individual event likelihoods given by Equation (13).

Furthermore, if we fix M_max and assume a prior ${p}_{0}(\alpha | {M}_{\max })$, we can calculate the Bayesian evidence in favor of a given M_max:

$$\begin{eqnarray}&&p({d}_{i}| {M}_{\max })=\int p({d}_{i}| \alpha ,{M}_{\max }){p}_{0}(\alpha | {M}_{\max })d\alpha .\end{eqnarray} \tag{ 15 }$$

We can then calculate the Bayes factor between two power-law models that differ in their choice of M_max. Recall that the default LIGO analysis fixes ${M}_{\max }=100\,{M}_{\odot }-{M}_{\min }$. For a sample of N-detected BBHs (assumed to be independent), the cumulative Bayes factor $K({M}_{\max },100\,{M}_{\odot })$ between a power-law model that fixes ${M}_{\max }=M$ and one that fixes ${M}_{\max }=100\,{M}_{\odot }$ is a product of the single-event evidence ratios:

$$\begin{eqnarray}&&K(M,100\,{M}_{\odot })=\displaystyle \prod _{i=1}^{N}\displaystyle \frac{p({d}_{i}| {M}_{\max }=M)}{p({d}_{i}| {M}_{\max }=100\,{M}_{\odot })}.\end{eqnarray} \tag{ 16 }$$

We calculate the cumulative Bayes factor $K(M,100\,{M}_{\odot })$ in Section 4.2.

4.1. Non-detection of Heavy BBHs

We first give an example of how the detection of only four BBHs with primary masses ${m}_{1}\leqslant 40\mbox{--}50\,{M}_{\odot }$ is inconsistent with certain (possibly correct) power-law mass distributions unless a mass gap is imposed. For a given mass distribution, we can use Equation (9) to calculate the probability of not detecting any BBHs above a certain mass, ${n}_{\mathrm{high}}=0$, given that we have detected n_low BBHs below the cutoff mass. To do this, we must first compute the ratio $\langle {VT}{\rangle }_{\mathrm{low}}/\langle {VT}{\rangle }_{\mathrm{high}}$ as defined in Equation (7) and then apply Equation (9). The results of this calculation for $p({n}_{\mathrm{high}}=0| {n}_{\mathrm{low}}=4)$ and $p({n}_{\mathrm{high}}=0| {n}_{\mathrm{low}}=10)$ are displayed in Figure 2. We show the results for two choices of cutoff mass: ${M}_{\mathrm{cutoff}}=41\,{M}_{\odot }$ (green curve) is motivated by the 95% credible upper limit on the primary mass of GW150914, the heaviest BBH detected, and ${M}_{\mathrm{cutoff}}=50\,{M}_{\odot }$ (blue and orange curves) is motivated by PPISN and PISN supernova models, which predict a mass gap starting at $40\mbox{--}50\,{M}_{\odot }$ (depending also on details of binary evolution; Belczynski et al. 2016a; Woosley 2017). We also vary the maximum total mass, ${M}_{\mathrm{tot},\max }$, of the "high-mass" population between $100\,{M}_{\odot }$ (blue and green curves), which is currently the maximum total mass that the aLIGO search includes in the definition of a stellar-mass BBH, to $200\,{M}_{\odot }$ (orange curve). We note that a power law with slope $\alpha =1$ is the "flat in log" population that the LIGO-Virgo Collaboration uses to compute the lower limits on the BBH merger rate, and we calculate that unless a mass gap is imposed, detecting four BBHs with primary masses ${m}_{1}\lt 41\,{M}_{\odot }$ is inconsistent with this population at the 96% level if we restrict ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$, or at the $\gt 99.9 \%$ level if we assume that the BBH population and the detectors' sensitivity extends up to ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$. Furthermore, unless a mass gap is imposed, there is already some tension (inconsistency at the 93% level) with the $\alpha =2.35$ population that LIGO-Virgo uses to compute the upper rate limits if we assume ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$. If the BBH mass distribution has an upper cutoff at $40\mbox{--}50\,{M}_{\odot }$, the inferred merger rates calculated without assuming the cutoff would be 1.4–2.1 times higher for the "flat in log" population and 1.1–1.4 times higher for the $\alpha =2.35$ population.

Zoom In Zoom Out Reset image size

4.2. Bayesian Evidence in Favor of Mass Gap

We have seen that assuming a single power-law mass distribution over the entire mass range $5\mbox{--}100\,{M}_{\odot }$ can rule out shallow power-law slopes in the absence of detections with component masses ${m}_{1}\gt 40\mbox{--}50\,{M}_{\odot }$. The absence of high-mass detections will continue to push the inferred power-law slope to steeper values unless we allow for an upper mass gap in the analysis. To study this point further, we simulate mock BBH mass measurements from a power-law population with slope $\alpha =2.35$ and an upper mass cutoff at ${M}_{\max }=41\,{M}_{\odot }$ (Equation (11), but we follow the canonical analysis used by the LVC (see Equation (3)) and fix ${M}_{\max }=95\,{M}_{\odot }$ and ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$ when inferring the power-law slope. While the bias on the inferred slope α may be small with ${ \mathcal O }(10)$ detections, with ${ \mathcal O }(100)$ detections, the canonical analysis will rule out the correct power-law slope (see Figure 3). (Although for hundreds of detections, a non-parametric fit to the mass distribution should be considered.) If we follow the canonical LIGO-Virgo analysis but set ${M}_{\max }=100\,{M}_{\odot }$ and ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$ rather than ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$, the presence of a mass gap will bias the power-law inference even more significantly, as the true population has ${M}_{\mathrm{tot},\max }=82\,{M}_{\odot }$. These results show that failing to account for an upper mass gap may lead to incorrect conclusions about the low-mass distribution. While we demonstrated this for an assumed power-law model, this caveat applies to any parameterized fit to the BBH mass distribution.

Zoom In Zoom Out Reset image size

With the first four LIGO BBH detections, varying M_max does not drastically bias the inference on α when fitting Equation (3) (see the solid lines in Figure 3). However, we can distinguish the model favored by the data by calculating the cumulative Bayes factor. Following Equations (15)–(16), we calculate this factor between two single-parameter power-law models with different fixed values of M_max. We choose to compare two cutoff values, $41\,{M}_{\odot }$ (the 95% upper limit on the heaviest component BH detected) and $100\,{M}_{\odot }$, and take the prior ${p}_{0}(\alpha | {M}_{\max })$ in Equation (15) to be a top hat over the wide range $-2\lt \alpha \lt 7$.

For the first four BBH detections, $K(41\,{M}_{\odot },100\,{M}_{\odot })=13$ if we assume the detectable BBH population only extends to ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$ (so that ${M}_{\max }=100\,{M}_{\odot }$ is really ${M}_{\max }\,=95\,{M}_{\odot }$). If we instead assume that the ${M}_{\max }=100\,{M}_{\odot }$ population is fully detectable up to total binary masses of ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$, the Bayes factor increases to $K(41\,{M}_{\odot },100\,{M}_{\odot })=90$, suggesting that there is already strong support for an upper mass cutoff at ${M}_{\max }\sim 40\,{M}_{\odot }$ over a cutoff at ${M}_{\max }\sim 100\,{M}_{\odot }$ within the assumed power-law model (Kass & Raftery 1995). The Bayes factor also depends on the choice of prior on α. We choose to be relatively uninformative in our prior, excluding only very steeply declining mass distributions ($\alpha \gt 7$) and allowing for moderately upward sloping mass distributions ($-2\lt \alpha \lt 0$), but it is clear from Figure 2 that a prior that favors large positive values of α (steeply declining power-law slopes) will lower the evidence in favor of a mass cutoff ${M}_{\max }\lt 100\,{M}_{\odot }$, while placing greater prior support on low values of α (shallow or downward sloping power laws) will raise the evidence in favor of a mass cutoff. If we enforce $\alpha \gt 0$ in the prior in order to agree with other astrophysical mass distributions, the Bayes factors change to $K(41\,{M}_{\odot },100\,{M}_{\odot })=5$ if we restrict ${M}_{\mathrm{tot},\max }\leqslant 100\,{M}_{\odot }$ or $K(41\,{M}_{\odot },100\,{M}_{\odot })=21$ if we allow ${M}_{\mathrm{tot},\max }\leqslant 200\,{M}_{\odot }$.

We anticipate that a set of 10 BBH detections with primary component masses ${m}_{1}\leqslant 41\,{M}_{\odot }$ will yield a Bayes factor $K(41\,{M}_{\odot },100\,{M}_{\odot })\gt 150$, providing very strong evidence for an upper mass gap. We assume that the underlying BBH population (and aLIGO's sensitivity) extends to total masses of $2{M}_{\max }$ in either case, so ${M}_{\mathrm{tot},\max }\leqslant 200\,{M}_{\odot }$. We take a flat prior on α in the range $[-2,7]$. With 191 events from a simulated $\alpha =2.35$, ${M}_{\max }=41\,{M}_{\odot }$ population, the single-event evidence ratios range from $K(41\,{M}_{\odot },100\,{M}_{\odot })=1.0$ to $K(41\,{M}_{\odot },100\,{M}_{\odot })=16.9$, with a median of $K(41\,{M}_{\odot },100\,{M}_{\odot })=2.6$. With a subset of 10 BBH detections from this population, we get $K(41\,{M}_{\odot },100\,{M}_{\odot })\gt 150$ in more than 99% of cases. If we detect 10 BBHs with primary component masses ${m}_{1}\leqslant 50\,{M}_{\odot }$, we likewise expect very strong evidence for a mass cutoff, with $K(50\,{M}_{\odot },100\,{M}_{\odot })\gt 150$ more than 95% of the time. The Bayes factor only compares two values of the mass cutoff; we fit for the value of M_max favored by a given set of detections in Section 4.3.

4.3. Joint Power Law–Maximum Mass Fit

In this section, we fit the two-parameter mass distribution of Equation (11). We calculate the likelihood $p({\boldsymbol{d}}| \alpha ,{M}_{\max })$ as the product of individual event likelihoods in Equation (13). We take flat priors on α and M_max, with $-2\leqslant \alpha \leqslant 7$ as before and ${M}_{\max }\leqslant 100\,{M}_{\odot }$. The minimum allowed value of M_max for a given set of detections is set by the lower-mass bound of the heaviest detected component BH. For simplicity, we take the lower-mass bound to be the minimum posterior sample. The upper bound ${M}_{\max }\leqslant 100\,{M}_{\odot }$ is motivated by the LIGO stellar-mass BBH search as well as by population-synthesis studies, which usually predict that the BH mass distribution would extend to $80\mbox{--}130\,{M}_{\odot }$ were it not for a pair-instability mass gap (Belczynski et al. 2016a; Eldridge & Stanway 2016; Spera et al. 2016). We calculate the likelihood function on a 500 × 100 grid of $(\alpha ,{M}_{\max })$ values in the allowed prior range and verify that increasing the resolution of the $(\alpha ,{M}_{\max })$ grid does not change our results. In fact, the resolution in the M_max dimension is limited by the finite number of posterior samples that are used to represent the component mass posterior PDFs for each event. To reduce these artificial discontinuities in the M_max dimension of the likelihood evaluation, we apply a two-dimensional smoothing spline before displaying the results.

The results of the joint power law–maximum mass analysis for the set of four detected BBHs is shown in Figure 4. We compute the joint likelihood twice: once fixing ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$, so that the population of stellar-mass BBHs is restricted to total masses ${M}_{\mathrm{tot}}\leqslant 100\,{M}_{\odot }$ regardless of M_max (top right panel) and once fixing ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$, so that the maximum total mass of the population is allowed to extend to ${M}_{\mathrm{tot},\max }=2{M}_{\max }$. We calculate the marginal posterior PDFs of α and M_max (top left and bottom right panels) under the assumption of a uniform prior on α in the range $[-2,7]$ and a uniform prior on M_max in the range $[29\,{M}_{\odot },100\,{M}_{\odot }]$ ($29\,{M}_{\odot }$ is the minimum posterior sample we generated for the primary component mass of GW150914).

Zoom In Zoom Out Reset image size

It is clear that properly accounting for our uncertainty on M_max when fitting the power-law mass distribution increases the support for shallow power-law slopes that would otherwise be ruled out under the assumption that the mass distribution extends continuously to $\sim 100\,{M}_{\odot }$. Allowing for freedom in M_max shifts the preferred values of α to shallower slopes, even allowing for negative α, as compared to the canonical analysis that fixes ${M}_{\max }=95\,{M}_{\odot }$ (orange dotted–dashed curve in Figure 4). Furthermore, the first four BBH detections already start to constrain M_max. The marginal posterior PDF $p({M}_{\max })$ peaks strongly at ${M}_{\max }\sim 40$, and the 95% upper limits on the inferred $p({M}_{\max })$ are $76.6\,{M}_{\odot }$ if assuming ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$ (or $90.7\,{M}_{\odot }$ if we conservatively assume ${M}_{\mathrm{tot},\max }=100\,{M}_{\odot }$). Taking ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$ rather than $100\,{M}_{\odot }$ allows the detectable BBH population to extend to $2{M}_{\max }$, thereby increasing the expected sensitivity to BBHs with primary component masses ${M}_{\max }\gt 50\,{M}_{\odot }$. Thus, the non-detection of heavy BBHs yields tighter constraints on the inferred M_max when we assume ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$, but the peak of the M_max distribution remains unchanged.

To explore the impact of future detections on the inferred mass distribution, we repeat this analysis for three simulated BBH populations, two with a power-law slope $\alpha =2.35$ and one with a power-law slope $\alpha =1$. One of the $\alpha =2.35$ populations, as well as the $\alpha =1$ population, has a mass gap starting at ${M}_{\max }=50\,{M}_{\odot }$, while the other $\alpha =2.35$ population has a mass gap starting at ${M}_{\max }=40\,{M}_{\odot }$.

The results of this calculation, assuming sensitivity up to ${M}_{\mathrm{tot},\max }=200\,{M}_{\odot }$, are shown in Figure 5, where each column corresponds to 40 detections from one of the three simulated populations. We see that 40 detections yield strong constraints on both the slope and maximum mass of the population. If the true population has a cutoff at ${M}_{\max }=40\,{M}_{\odot }$ (right column in Figure 5) rather than $50\,{M}_{\odot }$, we get tighter constraints on M_max. Similarly, we expect better constraints on the maximum mass for shallower mass distributions, as the non-detection of heavy BHs is more striking for shallow mass distributions (see Figure 2). The maximum mass is indeed better constrained for the population with power-law slope $\alpha =1$ (middle column) than for the population with the same mass cutoff ${M}_{\max }=50\,{M}_{\odot }$ but steeper power-law slope $\alpha =2.35$ (left column). As the first four LIGO detections currently favor ${M}_{\max }\lt 50\,{M}_{\odot }$ and moderately shallow power-law slopes $\alpha \lt 3$ (Figure 4), we expect that fewer than 40 detections will strongly constrain M_max.

Zoom In Zoom Out Reset image size

5.1. Effect of Redshift Evolution

5.2. Distribution of Mass Ratios

In fitting the mass distribution of BBHs (Equation (11)), we assumed that the distribution of mass ratios, q, is uniform in the allowed range ${M}_{\min }/{m}_{1}\lt q\lt \min ({M}_{\mathrm{tot},\max }/{m}_{1}-1,1)$. In particular, we assumed that for a given primary component mass, m₁, the marginal distribution of ${m}_{2}={{qm}}_{1}$ is given by Equation (2). However, many BBH formation models predict a preference for equal-mass mergers (Dominik et al. 2012; Rodriguez et al. 2016). To explore the effects of our assumed mass ratio distribution, we generalize Equation (2):

$$\begin{eqnarray}&&p({m}_{2}| {m}_{1})\propto \displaystyle \frac{{m}_{2}^{k}}{\min {({m}_{1},{M}_{\mathrm{tot},\max }-{m}_{1})}^{k+1}-{M}_{\min }^{k+1}},\end{eqnarray} \tag{ 17 }$$

so that k = 0 reduces to Equation (2) while $k\gt 0$ favors more equal-mass ratios. We find that the choice of $k\geqslant 0$ does not noticeably impact our results, and we recover consistent posteriors on $(\alpha ,{M}_{\max })$ if we fix k = 6 rather than k = 0. However, we note that there is currently no evidence that the distribution of mass ratios, $p(q| {m}_{1})$, deviates from the assumed uniform distribution. Although all of the events so far are consistent with mass ratios close to unity, this is not surprising given the selection effects that favor more equal-mass systems. For a fixed primary mass, m₁, assuming full matched-filter sensitivity, we would expect five detections with $q\gt 0.5$ for every detection with $q\lt 0.5$, and two detections with $q\gt 0.7$ for every detection with $q\lt 0.7$, even if we take q to be uniform in the range $[0,1]$ rather than $[{M}_{\min }/{m}_{1},1]$. We can explicitly check if the data favor $k\gt 0$ if we incorporate Equation (17) into the power-law mass distribution, so that Equation (11) becomes

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}| \alpha ,{M}_{\max })\propto \displaystyle \frac{{m}_{1}^{-\alpha }{m}_{2}^{k}\ { \mathcal H }({M}_{\max }-{m}_{1})}{\min {({m}_{1},{M}_{\mathrm{tot},\max }-{m}_{1})}^{k+1}-{M}_{\min }^{k+1}}.\end{eqnarray} \tag{ 18 }$$

We fit the above Equation (18) for k, marginalizing over α and M_max, and find that, for the first four LIGO detections, the likelihood $p({\boldsymbol{d}}| k)$ peaks mildly at k = 0, but is very broad. Thus, the first four BBHs mildly favor a uniform distribution of mass ratios. Future detections will continue to test this assumption.

5.3. Extending to Non-power-law Mass Distributions

5.4. Are there BBHs Beyond the Gap?

We have shown that given LIGO's extremely high sensitivity to BBHs with component masses $40\leqslant {m}_{1}\leqslant 100\,{M}_{\odot }$, it is statistically significant that the first four detections have been less massive than $40\,{M}_{\odot }$. We present a two-parameter model for the BBH mass distribution, consisting of a power law with slope α and a cutoff at M_max, and find that the first four detections already provide evidence for a cutoff to the mass distribution at ${M}_{\max }\sim 40\,{M}_{\odot }$. This cutoff may be the lower edge of a PPISN/ PISN upper mass gap. Furthermore, LIGO-Virgo have recently announced two more BBH detections (GW170608 and GW170814), both of which are less massive than $40\,{M}_{\odot }$ and only strengthen our conclusions (Abbott et al. 2017c, 2017d). We find that within ${ \mathcal O }(10)$ BBH detections, the location of the bottom edge of the upper mass gap will be significantly constrained. Our model assumes that all BBHs belong to a single population described by the same power law, so that the detection of a binary with mass in the mass "gap" would reset the lower edge of the gap beyond the mass of the newly detected binary. However, we expect to quickly converge on the true maximum mass of the population within ≲40 detections. At this point, the detection of a binary in the mass gap will be statistically inconsistent with this single population, and may indicate a subpopulation of BBHs that did not form directly from stellar-collapse (e.g., primordial BHs or BHs formed through previous mergers). The BBH spin distribution will provide further constraints on the existence of these subpopulations and will allow us to measure the fraction of BBHs that formed through previous mergers (Fishbach et al. 2017; Gerosa & Berti 2017).

We thank Will Farr for valuable discussions. M.F. was supported by the NSF Graduate Research Fellowship Program under grant No. 1746045. M.F. and D.E.H. were partially supported by NSF CAREER grant PHY-1151836 and NSF grant PHY-1708081. They were also supported by the Kavli Institute for Cosmological Physics at the University of Chicago through NSF grant PHY-1125897 and an endowment from the Kavli Foundation. D.E.H. thanks the Niels Bohr Institute for its hospitality while part of this work was completed and acknowledges the Kavli Foundation and the DNRF for supporting the 2017 Kavli Summer Program.