Picky Partners: The Pairing of Component Masses in Binary Black Hole Mergers

TL;DR

This paper analyzes binary black hole mergers and finds that black holes tend to pair with similar masses, providing insights into their formation mechanisms based on LIGO/Virgo data.

Contribution

It demonstrates that black hole pairings are non-random and favor similar masses, with a statistical model supporting mass ratio dependence in binary formation.

Findings

Black hole pairs are more likely to have similar masses than random pairing.

The median mass ratio of black hole binaries is approximately 0.91.

Most detected binaries have mass ratios greater than 0.5.

Abstract

We examine the relationship between individual black hole (BH) masses in merging binary black hole (BBH) systems. Analyzing the ten BBH detections from LIGO/Virgo's first two observing runs, we find that the masses of the component BHs comprising each binary are unlikely to be randomly drawn from the same underlying distribution. Instead, the two BHs of a given binary prefer to be of comparable mass. We show that it is $\sim 5$ times more likely that the component BHs in a given binary are always equal (to within 5\%) than that they are randomly paired. If we assume that the probability of a merger between two BHs scales with the mass ratio $q$ as $q^\beta$, so that $\beta=0$ corresponds to random pairings, we find $\beta>0$ is favored at credibility 0.987. By modeling the mass distribution, we find that the median mass ratio is $q_{50\%} = 0.91^{+0.05}_{-0.17}$ at 90\% credibility. While the pairing between BHs depends on their mass ratio, we find no evidence that it depends on the total mass of the system: it is $\sim 6$ times more likely that the pairing depends purely on the mass ratio than on the total mass. We predict that 99\% of BBHs detected by LIGO/Virgo will have mass ratios $q > 0.5$. We conclude that merging black holes do not form random pairings; instead they are selective about their partners, preferring to mate with black holes of a similar mass. The details of these selective pairings provide insight into the underlying formation channels of merging binaries.

Figures

Equations

$$p(m)\propto m^{\gamma },m_{\min }<m<m_{\max },$$

$$p(m_1,m_2\mid \gamma ,m_{\min },m_{\max })=\{\begin{array}{ll}\frac{2{\left(\gamma +1\right)}^2}{{\left(m_{\max }^{\gamma +1}-m_{\min }^{\gamma +1}\right)}^2}{m}_1^{\gamma }{m}_2^{\gamma }& \text{mbox}if{m}_{\min }<m_2<m_1<m_{\max }\\ 0& \text{mbox}else.\end{array}$$

$$p(m_1\mid \gamma ,m_{\min },m_{\max })=\frac{2\left(\gamma +1\right)m_1^{\gamma }\left(m_1^{\gamma +1}-m_{\min }^{\gamma +1}\right)}{{\left(m_{\max }^{\gamma +1}-m_{\min }^{\gamma +1}\right)}^2},$$

$$p(m_2\mid \gamma ,m_{\min },m_{\max })=\frac{2\left(\gamma +1\right)m_2^{\gamma }\left(m_{\max }^{\gamma +1}-m_2^{\gamma +1}\right)}{{\left(m_{\max }^{\gamma +1}-m_{\min }^{\gamma +1}\right)}^2}.$$

$$p(m_1,m_2\mid \gamma ,m_{\min },m_{\max },\vec{\beta })\propto \{\begin{array}{ll}{m}_1^{\gamma }{m}_2^{\gamma }{f}_p(\frac{m_2}{m_1},m_1+m_2\mid \vec{\beta })& \text{mbox}if{m}_{\min }<m_2<m_1<m_{\max }\\ 0& \text{mbox}else.\end{array}$$

$$f_p(q,M_{\mathrm{tot}}\mid \vec{\beta }=\{{\beta }_q,q_{\min }\})\propto \{\begin{array}{ll}{q}^{{\beta }_q}& \text{mbox}ifq>q_{\min }\\ 0& \text{mbox}else.\end{array}$$

$$f_p(q,M_{\mathrm{tot}}\mid \vec{\beta }=\{{\beta }_M\})\propto M_{\mathrm{tot}}^{{\beta }_M}.$$

$$f_p(q,M_{\mathrm{tot}}\mid \vec{\beta }=\{{\beta }_q,q_{\min },{\beta }_M\})\propto \{\begin{array}{ll}{q}^{{\beta }_q}{M}_{\mathrm{tot}}^{{\beta }_M}& \text{mbox}ifq>q_{\min }\\ 0& \text{mbox}else.\end{array}$$

$$\begin{array}{rl}p(m_1,m_2& \mid \alpha ,\beta ,m_{\min }{m}_{\max })=\\ & (\alpha +1)(\beta +1)\frac{m_1^{\alpha }}{m_{\max }^{\alpha +1}-m_{\min }^{\alpha +1}}\frac{m_2^{\beta }}{m_1^{\beta +1}-m_{\min }^{\beta +1}}.\end{array}$$

$$p(m_1,m_2\mid d)=\int p(m_1,m_2\mid \theta )p(\theta \mid d)\mathrm{d}\theta .$$

$$p(d\mid \theta )\propto e^{-\mu (\theta )}\prod _{i=1}^{N_{\mathrm{obs}}}\int p(d_i|m_1,m_2)\frac{\mathrm{d}\mathcal{R}}{\mathrm{d}{m}_1{m}_2}\left(\theta \right)\mathrm{d}{m}_1{m}_2,$$

$$p(d\mid \theta )\propto \prod _{i=1}^{N_{\mathrm{obs}}}\frac{\int p(d_i\mid m_1,m_2)p(m_1,m_2\mid \theta )\mathrm{d}{m}_1\mathrm{d}{m}_2}{\int VT(m_1,m_2)p(m_1,m_2\mid \theta )\mathrm{d}{m}_1\mathrm{d}{m}_2}.$$

$$p(m_1,m_2,z)\propto d_L(z{)}^2(1+z{)}^2\left(d_C(z)+\frac{(1+z)d_H}{E(z)}\right),$$

$$q_{\min }=m_{\min }/m_{\max }+q_{\mathrm{scale}}(0.95-m_{\min }/m_{\max }),$$

Full text

1 Introduction

Following the first two observing runs (O1 and O2) of advanced LIGO (Aasi et al., 2015) and Virgo (Acernese et al., 2015), the LIGO/Virgo Collaboration (LVC) reported ten detections of merging binary black holes (BBH) (Abbott et al., 2018c), with tens more detections expected from the third observing run (O3), and hundreds of expected detections per year once the gravitational-wave (GW) detector network reaches design sensitivity (Abbott et al., 2018a).

The formation and history of these BBHs remains a fundamental question in GW astrophysics. The proposed formation channels include isolated (Dominik et al., 2015; Belczynski et al., 2016a; Woosley, 2016; Eldridge & Stanway, 2016; Belczynski et al., 2016b; Stevenson et al., 2017; Kruckow et al., 2018; Spera et al., 2019a), dynamical (Mapelli, 2016; Hurley et al., 2016; Rodriguez et al., 2016b; Askar et al., 2017; Chatterjee et al., 2017; Rodriguez et al., 2018; Samsing, 2018; Zevin et al., 2019; Di Carlo et al., 2019), and primordial (Bird et al., 2016; García-Bellido, 2017), with many variants within each model. Different formation channels are expected to leave an imprint on the properties of the BBH population (Barrett et al., 2018; Taylor & Gerosa, 2018; Arca Sedda & Benacquista, 2019), including the mass distribution (Stevenson et al., 2015; Zevin et al., 2017), spin distribution (Rodriguez et al., 2016c; Farr et al., 2017; Vitale et al., 2017b; Farr et al., 2018), and redshift evolution (Fishbach et al., 2018; Vitale & Farr, 2018; Rodriguez & Loeb, 2018). It is therefore possible to learn about the astrophysics of BBH formation by fitting for these population distributions using GW data. In Abbott et al. (2018b), the LVC carried out such an analysis on the first ten BBH detections, fitting the mass, spin and redshift distributions with simple parameterized models. For example, the mass distribution was fit to a model in which the primary mass (the more massive component of a binary) follows a power-law between some minimum and maximum mass, while the secondary mass is distributed with a power-law between the minimum mass and its primary mass partner. Abbott et al. (2018b) additionally considered a slightly more complex model, which replaces the minimum mass cutoff with tapering at the low-mass end and allows for an additional Gaussian component at the high mass end of the primary-mass power-law. In this work, we restrict the population analysis to the ten Abbott et al. (2018c) BBH detections, as the detection efficiency has been previously studied for this sample and is well-understood. Using the wrong detection efficiency leads to selection biases in population inference. In future work we will extend our analysis to include overlapping samples with differing selection effects.

In this work we extend the analysis of Abbott et al. (2018b) by focusing on a particular aspect of the BBH mass distribution: the pairing between the two component BHs in the binary. We ask whether the universe makes merging binary black hole systems by randomly pairing up black holes, or whether the mass of each black hole in a pair influences the mass of its companion. This differs from the analysis of Abbott et al. (2018b), in which the parameterization for the mass distribution does not separate the underlying BH mass distribution and the pairing function. Under the models considered by Abbott et al. (2018b), it is not possible to fit for an underlying mass distribution that is common to both component BHs or quantify the deviation from the random-pairing scenario, as we do in this work.

We expect that the pairing function carries an imprint of the physics by which component BHs find their partners. Despite the different physical processes involved, many formation channels predict a preference for similar component masses (Mandel & Farmer, 2018). Binaries formed via homogeneous chemical evolution are expected to strongly prefer equal mass components . The traditional isolated evolution channel is also expected to favor comparable mass components, . However, the common envelope phase remains poorly understood, and this channel can produce mergers between fairly unequal component masses, especially at lower metallicities (Dominik et al., 2012; Stevenson et al., 2017; Klencki et al., 2018; Giacobbo et al., 2018; Spera et al., 2019b). Some studies have suggested that dynamical evolution also tends to produce more mergers with equal mass components due to the fact that comparable mass binaries have a higher binding energy and form tighter binaries (Rodriguez et al., 2016a; Amaro-Seoane & Chen, 2016). However, other dynamical channels may mildly prefer unequal mass components (Michaely & Perets, 2019). Alternatively, it has been suggested that in dynamical channels, the merger probability depends on the total mass, rather than the mass ratio, as mass segregation and dynamical interactions may favor binaries with larger total masses (O’Leary et al., 2016; Perna et al., 2019). Constraining the BBH pairing function with GW observations allows us to test these different predictions.

The pairing function has been previously studied in the context of the initial mass function for binary stars, where the degree of correlation between component stars (and the dependence on the orbital separation) remains an open question (Pinsonneault & Stanek, 2006; Kouwenhoven et al., 2009; Kroupa et al., 2013; Moe & Di Stefano, 2017; Kroupa & Jerabkova, 2018). It is possible that studying the pairing function for merging BBHs may shed light on the masses of their stellar progenitors, although the relationship between a BH’s mass and its progenitor star’s zero-age main-sequence (ZAMS) mass is complicated by the many stages of evolution undergone by BBHs.

In the stellar context, it has been pointed out that different pairing functions affect the 1-dimensional distribution of mass ratios as well as the 1-dimensional distributions of primary and secondary masses (Kouwenhoven et al., 2009). Because the primary (secondary) mass is defined to be the more (less) massive component in the binary, even randomly drawing two components from the same underlying distribution results in the primary and secondary masses having different distributions. Random draws can also result in very different mass ratio distributions, depending on the shape of the underlying mass distribution. We emphasize that the pairing mechanism cannot be determined by examining any one of these one-dimensional distributions independently. For example, a mass ratio distribution that favors near-unity mass ratios may simply indicate that the underlying BH mass distribution peaks in a narrow mass range, rather than that similar component masses are more likely to partner and merge. It is therefore important to examine the two-dimensional mass distribution in order to analyze whether or not there is a preference for similar-mass components.

This paper explores the BBH pairing function by analyzing the first ten LIGO/Virgo BBHs according to the mass models described in Section 2. The results of the analysis and implications for future detections are found in Section 3. In Section 4 we demonstrate the analysis on mock GW data and forecast the constraints that will be possible with $\sim 50$–$100$ more BBH detections (to be expected at the end of O3 or shortly after the start of O4). We conclude in Section 5. Appendix A describes the details of the hierarchical Bayesian analysis.

2 Mass Distribution Models

In the simplest case, we consider a model in which the component masses in a BBH system are independently drawn from the same underlying power-law distribution:

[TABLE]

where $\gamma$ is the power-law slope, and $m_{\mathrm{min}}$ and $m_{\mathrm{max}}$ are the minimum and maximum mass. We refer to this as the “random pairing” mass distribution (Kouwenhoven et al., 2009). We note that in this case, the marginal distributions of the primary and secondary masses are not identical, because the primary (secondary) is defined as the more massive (less massive) component. Defining $m_{1}$ as the primary mass and $m_{2}$ as the secondary mass, the random pairing power-law distribution takes the form: {widetext}

[TABLE]

This implies that the primary masses follow the distribution:

[TABLE]

while the secondaries follow:

[TABLE]

We reiterate that the distributions in Eqs. 3 and 4 are not the same as the underlying distribution (Eq. 1), even though both masses are separately drawn from this distribution. In particular, the primary mass distribution will tend to favor larger masses compared to the secondary mass distribution. Furthermore, different choices of the underlying power-law parameters ($\gamma$, $m_{\mathrm{min}}$, and $m_{\mathrm{max}}$) will lead to different distributions in the mass ratio $q\equiv m_{2}/m_{1}\leq 1$. If the underlying power-law is steep enough in either direction, mass ratios close to unity will be favored even if the two components are randomly paired.

In order to explore the pairing of two component BHs, we consider mass distributions that contain the random pairing distribution as a sub-model, but allow for deviations parameterized by a pairing function, $f_{p}$.

If each BH mass is drawn from an underlying power-law distribution with slope $\gamma$, and the probability of two masses belonging to a merging binary is given by $f_{p}(q,M_{\mathrm{tot}}\mid\vec{\beta})$, where $\vec{\beta}$ denotes the hyper-parameter(s) of the pairing function, the mass distribution of merging BHs follows: {widetext}

[TABLE]

We first consider pairing functions, $f_{p}$, that depend purely on the mass ratio. As a simple model, we assume that the pairing function follows a power-law in mass ratio with slope $\beta_{q}$, with a minimum mass ratio threshold required for merger, $q_{\mathrm{min}}$. In this model:

[TABLE]

This model reduces to the random pairing model for $\left(\beta_{q},q_{\mathrm{min}}\right)=\left(0,\frac{m_{\mathrm{min}}}{m_{\mathrm{max}}}\right)$.

Second, we consider pairing functions $f_{p}$ that depend on the total mass of the system. A simple model in this case, motivated by the predictions of the dynamical channel of O’Leary et al. (2016) and the primordial BH channel of Kocsis et al. (2018), is a simple power-law in $M_{\mathrm{tot}}$:

[TABLE]

More generally, we may consider a pairing function that depends on both the mass ratio and the total mass:

[TABLE]

We highlight that because we consider models that reduce to random pairing under some choice of parameters, our assumed parameterization differs from the mass distribution models analyzed in Abbott et al. (2018b). The basic power-law model in Abbott et al. (2018b) is defined such that the marginal $p(m_{1})$ distribution follows a power-law, so that the joint mass distribution takes the form (note we have redefined $\alpha$ from Abbott et al., 2018b, as $-\alpha$):

[TABLE]

On the other hand, for the parameterizations we consider in this work, the marginal distribution of primary masses does not follow an exact power-law; instead, these parameterizations allow for the possibility that both masses in a binary are drawn from the same underlying power-law distribution.

Following the methods laid out in Appendix A, we fit the pairing models discussed above to the first ten BBH detections in Section 3.

3 Results

3.1 LVC Model

We begin by recovering the results of Abbott et al. (2018b) under the same mass model given by Eq. 9, equivalent to Model B of Abbott et al. (2018b), in order to demonstrate consistency between our methods. Although we use slightly different assumptions regarding the spin distribution and the selection effects calculation (see Appendix A), we recover nearly-identical posterior distributions on the population hyper-parameters: $\alpha=-1.7^{+1.8}_{-1.5}$, $\beta=6.1^{+5.4}_{-5.7}$, $m_{\mathrm{min}}=7.3^{+1.7}_{-3.7}\ M_{\odot}$, and $m_{\mathrm{max}}=41.8^{+18.6}_{-5.5}\ M_{\odot}$. This is to be compared with $\alpha=-1.6^{+1.7}_{-1.5}$, $\beta=6.7^{+4.8}_{-5.9}$, $m_{\mathrm{min}}=7.9^{+1.2}_{-2.5}\ M_{\odot}$, and $m_{\mathrm{max}}=42.0^{+15.0}_{-5.7}\ M_{\odot}$ found by Abbott et al. (2018b). As described in Appendix A, our prior on $m_{\mathrm{min}}$ starts at 3 $M_{\odot}$ rather than 5 $M_{\odot}$. Furthermore, recall that our convention for the power-law slope $\alpha$ has a sign flip compared to the convention in Abbott et al. (2018b).

With the current set of events, the data cannot distinguish between the mass model of Abbott et al. (2018b) (Eq. 9) and the models we consider in this work, and they all give consistent results for the inferred mass distribution $p(m_{1},m_{2})$. However, the parameters of our models have a different interpretation from the mass model of Abbott et al. (2018b). While the power-law slope $\alpha$ of Eq. 9 refers to the power-law of the primary mass distribution, the power-law slope $\gamma$ of this work refers to the underlying mass distribution power-law from which both primary and secondary BHs are drawn. The additional parameters $\beta_{q}$, $q_{\mathrm{min}}$ and $\beta_{M}$ in our models allow us to explore whether the pairing between the two component masses is random or whether (and how) it depends on the mass ratio , according to Eq. 8).

3.2 Random Pairing

Figure 1 shows the results of fitting the random-pairing model (Eq. 2 with free parameters $\gamma$, $m_{\mathrm{min}}$, $m_{\mathrm{max}}$), the mass ratio power-law model (Eq. 6 with free parameters $\gamma$, $m_{\mathrm{min}}$, $m_{\mathrm{max}}$, $\beta_{q}$) and the total mass power-law (Eq. 7 with free parameters $\gamma$, $m_{\mathrm{min}}$, $m_{\mathrm{max}}$, $\beta_{M}$) to the ten BBHs from the first two observing runs. We use flat priors on all free parameters, with uninformative prior bounds listed in Appendix A. In each case, the parameters of Eq. 8 that are not left free are fixed to the default values $\beta_{q}=\beta_{M}=0$, $q_{\mathrm{min}}=m_{\mathrm{min}}/m_{\mathrm{max}}$.

If we fix the pairing to be random, we find $\gamma=-1.1^{+1.0}_{-0.9}$. However, as shown in the left-hand, middle panel in Figure 1, and explained in Section 2, this does not imply that the one-dimensional marginal distributions of the primary and secondary masses follow this common power-law; the primary masses follow a flatter distribution while the secondary masses follow a steeper distribution. Note that the inferred mass ratio distribution (bottom row, left-hand panel) in this case is inferred to be nearly flat across the range $\sim 0.15$–$1$. This is a consequence of this particular fit to the random-pairing model; in general the marginal mass ratio distribution can slope significantly upwards depending on the value of $\gamma$. The lower-limit on the mass ratio in the random-pairing model is given by the ratio $m_{\mathrm{min}}/m_{\mathrm{max}}$, and is constrained to $\sim 0.15$ in this case due to the measurements $m_{\mathrm{min}}\sim 7\ M_{\odot}$ and $m_{\mathrm{max}}\sim 40\ M_{\odot}$.

3.3 Mass Ratio Dependent Pairing

The effect of introducing $\beta_{q}$ as a free parameter, while still fixing $q_{\mathrm{min}}=m_{\mathrm{min}}/m_{\mathrm{max}}$ and $\beta_{M}=0$, is shown the middle column of Figure 1. Under this model extension, the data display a clear preference for mass ratios close to unity (bottom row, middle column), which implies more overlap between the primary and secondary mass distributions (middle row, central panel). We infer $\beta_{q}=7.0^{+4.5}_{-5.5}$, and find that $\beta_{q}\leq 0$ is ruled out with probability $0.987$. This suggests that the random pairing model ($\beta_{q}=0$) is strongly disfavored by the data. Meanwhile, the underlying mass distribution power-law slope is inferred to be a bit steeper than in the random-pairing case, with $\gamma=-1.3^{+0.9}_{-0.8}$ compared to $\gamma=-1.1^{+1.0}_{-0.9}$ for the random-pairing model. All models essentially agree on $m_{\mathrm{min}}\sim 6.7^{+1.9}_{-3.2}\ M_{\odot}$ and $m_{\mathrm{max}}\sim 41.9^{+18.2}_{-5.7}\ M_{\odot}$.

We fit the model with both $\beta_{q}$ and $q_{\mathrm{min}}$ left free in Figure 2, which displays the posterior distributions on the five hyper-parameters (three parameters to characterize the 1-dimensional mass distribution, and two – $\beta_{q}$ and $q_{\mathrm{min}}$ – to characterize the pairing function) as a corner plot (Foreman-Mackey, 2016). We find a strong preference for near-equal mass ratios, inferring that 99% of merging BBHs have mass ratios between $q_{1\%}=0.66^{+0.25}_{-0.28}$ and unity. For reference, under Model B from Abbott et al. (2018b), we find $q_{1\%}=0.59^{+0.13}_{-0.33}$. Meanwhile, we find that the median mass ratio is $q_{50\%}=0.91^{+0.05}_{-0.17}$. This agrees very closely with the findings of Roulet & Zaldarriaga (2019), who find a preference for population distributions with an average mass ratio $\bar{q}=0.89^{+0.08}_{-0.18}$. In fact, the current set of detections is consistent with all binaries consisting of equal component masses (to within 5%, with $q_{\mathrm{min}}=0.95$; the maximum value permitted by our prior). We find that $q_{\mathrm{min}}=0.95$ is five times more likely than $(\beta_{q},q_{\mathrm{min}})=(0,m_{\mathrm{min}}/m_{\mathrm{max}})$, meaning that it is five times more likely that all binaries consist of equal component masses than that they are randomly paired. Although we do not include the events of Venumadhav et al. (2019a, b) or Nitz et al. (2019) in our analysis (), we note that all of their detections are also consistent with mass ratios of unity. This suggests an even stronger preference for near-equal component masses in the underlying population.

3.4 Posterior predictive distributions

In the following, because there is no evidence that the data prefers nonzero $\beta_{M}$, we fix $\beta_{M}=0$ in the population model of Eq. 8 (reducing to the model of Eq. 6). Using the recovered posteriors on the population hyper-parameters, shown as a corner plot in Figure 2, we calculate the posterior population distribution $p(m_{1},m_{2}\mid\rm data)$, shown in Figure 3. We define the posterior population distribution as in Abbott et al. (2018b); this refers to the distribution of true mass values marginalized over the hyper-parameter posteriors for a given population model:

[TABLE]

Here, $p(\theta\mid d)$ refers to the posterior distribution on the population model’s hyper-parameters inferred from the ten BBH events. Figure 3 shows draws from the posterior population distribution on the true masses of the underlying BBH population (blue) as well as the true masses of detected systems (orange), found by applying selection effects to the underlying mass distribution. For comparison, the mass measurements of the O1 and O2 detections are shown in black.

Using the posterior population distribution on the masses $m_{1}$ and $m_{2}$, we calculate the corresponding distribution on the mass ratio $q$ in Figure 4.

We expect 90% of detected events to have their true masses fall within the orange credible region of Figure 3. In terms of the mass ratio distribution, we expect that 90% of detected events will have mass ratios $q>0.73$, and 99% of detected events will have mass ratios $q>0.51$. We can take these predictions one step further by simulating the measured mass ratio values for detected events, which accounts for measurement uncertainty in addition to the selection effects. Given the true masses of a detected event, we generate a mock posterior to represent how those masses would be measured in LIGO/Virgo data. The mock posteriors are generated according to the prescription described in Appendix B. We summarize the expected mass ratio posteriors from anticipated detections as the green dashed line (median) and shaded band (symmetric 90% interval) in Figure 4. We refer to this green band as the “posterior predictive process.” Based on the first ten detections, and assuming that all detections are described by the same population model assumed here, we expect that 90% of future detections will have recovered mass ratio posteriors that lie within the shaded band. We see that measurement uncertainty plays a significant role in shifting the observed mass ratio posteriors (with the default flat-in-component-mass priors) away from 1 relative to $q_{\mathrm{true}}$.

4 Simulations

We expect to have tens more BBH detections by the end of LIGO/Virgo’s third observing run in mid-2020, and hundreds of detections within a few more years (Abbott et al., 2018a). In this section, we explore the expected mass distribution constraints from tens to hundreds of detections under the models considered here. We perform our analysis on mock GW detections that we generate from known underlying distributions. We follow a simplified yet realistic method for generating mock measurements from the underlying population and ensure that the mock primary and secondary masses are measured with uncertainties typical to second-generation GW detectors (Vitale et al., 2017a). The method for generating mock detections is described in Appendix B.

The expected constraints from 60 detections (similar to what we expect by the end of O3; Abbott et al., 2018a) are shown in Figure 5 for a simulated population described by Eq. 6 with $m_{\mathrm{min}}=7\ M_{\odot}$, $m_{\mathrm{max}}=40\ M_{\odot}$, $\gamma=-1$, $\beta_{q}=6$, and q_{\mathrm{min}}=m_{\mathrm{min}}/m_{\mathrm{max}}{\color[rgb]{0,0,0}=0.175}.

We find with 60 events from this simulated population, we can typically rule out random-pairing with Bayes factors $\gtrsim 1000$. These projections are conservative because the deviations from random-pairing in the chosen mock population are not very large compared to the values of $q_{\mathrm{min}}$ and $\beta_{q}$ that are favored by the first ten events. The parameters that govern the pairing mechanism will become increasingly well-constrained, although with large correlations between them. The simulated 60 events shown here yield $\beta_{q}=3.4^{+4.6}_{-5.9}$ and $q_{\mathrm{scale}}=0.5^{+0.4}_{-0.4}$ ($q_{\mathrm{min}}=0.5^{+0.3}_{-0.3}$). Meanwhile, the parameters of the underlying power-law mass distribution will also become well-constrained. With 50 more events, we expect to constrain $m_{\mathrm{max}}$ to a couple of solar masses (this particular realization yields $m_{\mathrm{max}}=40.6^{+2.4}_{-1.7}$) and the power-law slope $\gamma$ to a 90% credible interval of $<1$ ($\gamma=-1.1^{+0.4}_{-0.3}$). Note that if $\beta_{M}$ is left free in the pairing function, the measurement of $\gamma$ will become less informative, as we constrain a linear combination of $\beta_{M}$ and $\gamma$ (see Section 3). With 100 events, these constraints will improve roughly as $1/N$ and $1/\sqrt{N}$ for $m_{\mathrm{max}}$ and $\gamma$, respectively111Because $m_{\mathrm{max}}$ is a sharp feature, its measurement converges faster than the typical $1/\sqrt{N}$ (see e.g. Chakrabarty et al., 2003; Johnson, 2007). This makes it particularly useful as a feature to constrain cosmology (Farr et al., 2019).: 100 simulated events gives $m_{\mathrm{max}}=38.9^{+1.4}_{-0.9}\ M_{\odot}$ and $\gamma=-0.8^{+0.3}_{-0.3}$. It may take more events for the constraints on $m_{\mathrm{min}}$ to become interesting, because the detector sensitivity is a steep function of BH mass, and most detections are at the high end of the mass function (Fishbach & Holz, 2017), or beyond (Fishbach et al., 2019). For a flat prior starting at $3\,M_{\odot}$, 60 mock events give $m_{\mathrm{min}}=6.5^{+1.3}_{-2.4}$. However, the constraints on $m_{\mathrm{min}}$, like the constraints on $\gamma$, are less informative if $\beta_{M}$ is allowed to vary.

5 Conclusion

We have fit the mass distribution of merging BBHs with a simple model that parameterizes the pairing function between the two components in a binary. We highlight the importance of comparing the full two-dimensional mass distribution of BBHs, because it is impossible to disentangle the overall BH mass distribution from the pairing function when considering only one-dimensional distributions of the primary/ secondary mass or the mass ratio.

Based on the first ten LIGO/Virgo BBH detections, we conclude that component BHs are not randomly paired in a binary; rather, the pairing likely favors components of comparable masses. We find that it is five times more likely that mergers only take place between equal (to within 5%) mass BHs than that component BHs are randomly drawn from the same underlying distribution. Our fits imply that 99% of mass ratios among the population of merging BBHs are greater (closer to equal mass) than $q_{\mathrm{1\%}}=0.66^{+0.25}_{-0.28}$. This is to be compared with an expected value of $q_{\mathrm{1\%}}=0.15^{+0.07}_{-0.06}$ for the random pairing scenario. We predict that among detected BBHs, 90% will have mass ratios $q>0.73$, and 99% will have mass ratios $q>0.51$. Meanwhile, we find no evidence that the pairing function depends on the total mass of the system, contrary to the predictions of some dynamical and primordial BH formation channels (O’Leary et al., 2006; Kocsis et al., 2018; Perna et al., 2019).

By the end of O3, the details of the pairing function will be better constrained (compare the joint posterior on $q_{\mathrm{min}}$ and $\beta_{q}$ in Figure 2—the current constraints—with Figure 5—the constraints we expect by the end of O3). We hope that these results will enable detailed comparisons with the predictions of the full 2-dimensional merger rate $\mathcal{R}(m_{1},m_{2})$ from population synthesis simulations.

As usual, our results rely on the assumption that there is a single population of BBHs that is adequately described by our simple parameterized model (see e.g. Doctor et al., 2019). One way to test the validity of this assumption with future detections is to compare them against the posterior predictive distribution (for example, Figures 3 and 4) inferred from the model. We conclude that the universe does not assemble its black-hole binaries at random, and future constraints of the pairing function we have introduced above will yield important insights into these formation processes.

We are grateful to the LIGO and Virgo collaborations for publicly releasing posterior samples on the Gravitational Wave Open Science Center. We thank Eve Chase and Reed Essick for their helpful comments on the manuscript, and Reed Essick for suggesting the concept of a posterior predictive process for Figure 4. MF was supported by the NSF Graduate Research Fellowship Program under grant DGE-1746045. MF and DEH were supported by 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. DEH also gratefully acknowledges support from the Marion and Stuart Rice Award.

Appendix A Methods

We carry out a hierarchical Bayesian analysis to fit the hyper-parameters for each of the mass models discussed in Section 2. We fit only for the distribution of primary and secondary masses, and fix the distributions of all other BBH intrinsic and extrinsic source parameters. We fix the underlying redshift distribution to follow a merger rate that is uniform in comoving volume and source-frame time. We assume that the underlying population is isotropic on the sky, with isotropic inclination angles. For definiteness we fix the spin distribution of both binary components to be uniform in spin magnitude and isotropic in spin tilt. Although this distribution is not necessarily favored by the data, the correlation between the inferred spin distribution and the inferred mass distribution is negligible, as shown in Abbott et al. (2018b), which fit simultaneously for the mass and spin distribution. In particular, despite using a different spin model, we recover the results of Abbott et al. (2018b) under the same mass model.

The likelihood is given by the inhomogeneous Poisson process likelihood (Loredo, 2004; Mandel et al., 2019; Abbott et al., 2018b). For $N_{\mathrm{obs}}$ independent events, the likelihood of the data $d$ given hyper-parameters $\theta$ is:

[TABLE]

where $p(d_{i}\mid m_{1},m_{2})$ denotes the likelihood of an individual event’s data given its component masses, $\frac{\mathrm{d}\mathcal{R}}{\mathrm{d}m_{1}m_{2}}\left(\theta\right)$ is the differential merger rate density, which integrates to the total merger rate density $\mathcal{R}$ and is given by $\mathcal{R}p(m_{1},m_{2}\mid\theta)$, and $\mu(\theta)=\mathcal{R}\langle VT\rangle_{\theta}$ denotes the expected number of detections given $\mathcal{R}$ and the sensitive spacetime volume $\langle VT\rangle_{\theta}$ of the detector network to a given population of BBHs with hyper-parameters $\theta$.

We assume that the merger rate density is uniform in comoving volume and source-frame time, and calculate $\langle VT\rangle_{\theta}$ according to a semi-analytic prescription (Finn & Chernoff, 1993; Finn, 1996). Following Abbott et al. (2018b), we assume that a single-detector signal-to-noise (SNR) threshold of 8 is necessary and sufficient for detection, and that the detector’s noise curve is described by the Early High Sensitivity power spectral density (PSD) for advanced LIGO (Abbott et al., 2018a). The validity of these assumptions is discussed in Abbott et al. (2018b). Unlike in Abbott et al. (2018b), in this work we do not calibrate the $VT(m_{1},m_{2})$ to the results of injection campaigns into the detection pipelines. As we demonstrate by explicitly comparing our results to those of Abbott et al. (2018b) in Section 3, using the uncalibrated $VT$ calculation leads to a slight bias in our inference of the overall-merger rate, with the median shifting by a factor of $\sim 1.7$, as expected from Fig. 9 in Abbott et al. (2018b). However, this does not affect the inferred shape of the mass distribution, which is our primary interest in this work. We also neglect the effect of non-zero spins in the estimation of $VT$, as spins have a sub-dominant effect on the sensitivity (Abbott et al., 2018b), especially given that existing detections disfavor a significant population of highly spinning systems.

Note that if we marginalize over the rate density $\mathcal{R}$ with a flat-in-log prior, the likelihood takes the form (Fishbach et al., 2018; Mandel et al., 2019):

[TABLE]

To get the individual-event likelihood term $p(d_{i}\mid m_{1},m_{2})$ that appears in Eq. A1, we use the publicly available IMRPhenomPv2 posterior samples for the ten BBH detections in O1 and O2 (Gravitational Wave Open Science Center, 2018). There is a negligible difference between the mass posteriors derived with the IMRPhenomPv2 waveform (Husa et al., 2016; Khan et al., 2016) and the SEOBNRv3 waveform (Pan et al., 2014). The individual-event posteriors were calculated under priors that are flat in detector-frame masses, and “volumetric” in luminosity distance, $d_{L}$. In terms of source-frame masses and cosmological redshift $z$, the default event-level prior is therefore (Abbott et al., 2018b):

[TABLE]

where $d_{C}$ is the comoving distance and $d_{H}=c/H_{0}$ is the Hubble distance, and $E(z)=H(z)/H_{0}$ (Hogg, 1999). We divide out by these priors in our analysis to get a term that is proportional to the likelihood rather than the posterior. We fix the cosmological parameters to the best-fit Planck 2015 values (Ade et al., 2016) throughout for consistency with Abbott et al. (2018b) and Abbott et al. (2018c).

We sample from the overall likelihood of Eq. A1 using PyMC3 (Salvatier et al., 2016). In all models considered, we choose priors that are flat over $m_{\mathrm{min}}$, $m_{\mathrm{max}}$ and the power-law slope $\gamma$ within the ranges $3\ M_{\odot}<m_{\mathrm{min}}<10\ M_{\odot}$, $35\ M_{\odot}<m_{\mathrm{max}}<100\ M_{\odot}$ and $-4<\gamma<2$. We take a flat-in-log prior on the rate $p(\mathcal{R})\propto 1/\mathcal{R}$. Unless they are fixed to some value, we take a flat prior on $\beta_{q}$ in the range $-4<\beta_{q}<12$ and a flat prior on $\beta_{M}$ in the range $0<\beta_{M}<12$. Because the prior range of the minimum mass ratio $q_{\mathrm{min}}$ depends on two other free parameters, $m_{\mathrm{min}}$ and $m_{\mathrm{max}}$, we introduce another parameter $q_{\mathrm{scale}}$, defined so that:

[TABLE]

and, unless it is fixed, we sample over $q_{\mathrm{scale}}$ with a flat prior from 0 to 1, so that $m_{\mathrm{min}}/m_{\mathrm{max}}<q_{\mathrm{min}}<0.95$. We restrict the upper limit of $q_{\mathrm{min}}$ to slightly below 1 in order to avoid sampling issues, as the mass ratio of any individual GW event is measured with a finite resolution, and this prevents $q_{\mathrm{min}}$ from being resolved arbitrarily close to $q_{\mathrm{min}}=1$. When using Model B from Abbott et al. (2018b), we use their same priors, with the exception of lower prior boundary for $m_{\mathrm{min}}$, which we take to be $3\ M_{\odot}$ rather than $5\ M_{\odot}$.

For those models which contain random-pairing as a subset, we quantify the evidence for the random-pairing hypothesis versus the full model by calculating the Savage-Dickey density ratio (SDDR), which is defined as the ratio of the posterior probability to the prior probability at the given point in parameter space (Dickey, 1971).

Appendix B Simulated Detections

In generating mock detections, we assume that the underlying population follows a uniform in comoving volume and source-frame time merger rate, with isotropic sky positions and inclinations, and zero spins. The true component masses are drawn from the given population distribution. We note that the assumptions of fixed redshift and spin distributions are unlikely to affect the inference of the pairing function (mass ratios are measured independently of redshift, and excluding spins did not make a difference in the O1 and O2 analysis); however, these distributions can be fit jointly with the mass distribution and marginalized over (Abbott et al., 2018b).

Given the true parameters of the binary, we calculate the SNR of the signal in a single detector, assuming that the noise is described by the Mid-High Sensitivity PSD as expected for O3 for the LIGO detectors (Abbott et al., 2018a). We assume that the binary is then detected if it passes a single-detector SNR threshold of 8. In order to assign measured component masses to each detected binary, we assume that the fractional uncertainty on the source-frame chirp mass follows $\frac{\sigma_{\mathcal{M}}}{\mathcal{M}}=\frac{8}{\rho}\left(0.01+\left(\frac{0.2z}{1+z}\right)^{2}\right)^{1/2}$, where $z$ is the true redshift, while the uncertainty on the symmetric mass ratio $\eta\equiv\frac{m_{1}m_{2}}{(m_{1}+m_{2})^{2}}$ follows $\sigma_{\eta}=0.03\frac{8}{\rho}$, where $\rho$ is the single-detector SNR of the source. Given a true value of $\mathcal{M}$ and $\eta$ for each binary, we randomly draw $\mathcal{M}_{\mathrm{obs}}$ from a log-normal distribution centered on $\mathcal{M}$ with standard deviation $\sigma_{\mathcal{M}}$, and $\eta_{\mathrm{obs}}$ from a normal distribution centered on $\eta$ with standard deviation $\sigma_{\eta}$. With these values of $\mathcal{M}_{\mathrm{obs}}$ and $\eta_{\mathrm{obs}}$ and their assumed known distributions about the true chirp mass $\mathcal{M}$ and symmetric mass ratio $\eta$, we generate mock posterior samples for the component masses $m_{1}$ and $m_{2}$ under flat priors using the Monte-Carlo sampler PyStan (Carpenter et al., 2017). These uncertainties are typical of the O2 detections, and result in typical 90% measurement uncertainties on the source-frame component masses of $\approx 50\%$, with a distribution of uncertainties that matches the one in Vitale et al. (2017a).