Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo

I Introduction

The second LIGO/Virgo observing run (O2) spanned nine months between November 2016 through August 2017, building upon the first, four-month run (O1) in 2015. The LIGO/Virgo gravitational-wave (GW) interferometer network is comprised of two instruments in the United States (LIGO) (LIGO Scientific Collaboration et al. 2015; Abbott et al. 2016a) and a third in Europe (Virgo) (Acernese et al. 2015), the latter joining the run in the summer of 2017. In total, ten binary black hole (BBH) mergers have been detected to date (Abbott et al. 2018). The BBHs detected possess a wide range of physical properties. The lightest so far is GW170608 (Abbott et al. 2017a) with an inferred total mass of $18.7_{-0.7}^{+3.3}$$\mathit{M_{\odot}}$. GW170729 (Abbott et al. 2018)—exceptional in several ways—is likely to be the heaviest BBH to date, having total mass $85.2_{-11.2}^{+15.4}$$\mathit{M_{\odot}}$, as well as the most distant, at redshift $0.48_{-0.20}^{+0.19}$. Both GW151226 and GW170729 show evidence for at least one black hole with a spin greater than zero (Abbott et al. 2016b; Abbott et al. 2018).

By measuring the distributions of mass, spin, and merger redshift in the BBH population, we may make inferences about the physics of binary mergers and better understand the origin of these systems. We employ Bayesian inference and modelling (Gelman et al. 2004; Mandel 2010; Foreman-Mackey et al. 2014; Hilbe et al. 2017; ase 2018) which, when applied to parameterized models of the population, is able to infer population-level parameters — sometimes called hyperparameters to distinguish them from the event-level parameters — while properly accounting for the uncertainty in the measurements of each event’s parameters (Mandel 2010; Hogg et al. 2010).

The structure and parameterization of BBH populations models are guided by the physical processes and evolutionary environments in which BBH are expected to form and merge. Several BBH formation channels have been proposed in the literature, each of them involving a specific environment and a number of physical processes. For example, BBHs might form from isolated massive binaries in the galactic field through common-envelope evolution (Bethe & Brown 1998; Portegies Zwart & Yungelson 1998; Belczynski et al. 2002; Voss & Tauris 2003; Dewi et al. 2006; Belczynski et al. 2007; Belczynski et al. 2008; Dominik et al. 2013; Belczynski et al. 2014; Mennekens & Vanbeveren 2014; Spera et al. 2015; Tauris et al. 2017; Eldridge & Stanway 2016; Stevenson et al. 2017b; Chruslinska et al. 2018; Mapelli et al. 2017; Giacobbo et al. 2018; Mapelli & Giacobbo 2018; Kruckow et al. 2018; Giacobbo & Mapelli 2018) or via chemically homogeneous evolution (Marchant et al. 2016; de Mink & Mandel 2016; Mandel & de Mink 2016). Alternatively, BBHs might form via dynamical processes in stellar clusters (Portegies Zwart & McMillan 2000; Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Grindlay et al. 2006; O’Leary et al. 2006; Sadowski et al. 2008; Ivanova et al. 2008; Downing et al. 2010; Downing et al. 2011; Clausen et al. 2013; Ziosi et al. 2014; Rodriguez et al. 2015; Rodriguez et al. 2016a; Mapelli 2016; Askar et al. 2017; Banerjee 2017; Chatterjee et al. 2017) and galactic nuclei (Antonini & Perets 2012; Antonini & Rasio 2016; Petrovich & Antonini 2017), evolution of hierarchical triple systems (Antonini et al. 2014; Kimpson et al. 2016; Antonini et al. 2017; Liu & Lai 2018), gas drag and stellar scattering in accretion disks surrounding super-massive black holes (McKernan et al. 2012; Bartos et al. 2017; Stone et al. 2017). Finally, BBHs might originate as part of a primordial black hole population in the early Universe (Carr & Hawking 1974; Carr et al. 2016; Sasaki et al. 2016; Inomata et al. 2017; Inayoshi et al. 2016; Bird et al. 2016; Ali-Haïmoud et al. 2017; Clesse & García-Bellido 2017; Chen & Huang 2018; Ando et al. 2018), where their mass spectrum is typically proposed as having power law behavior, but spanning a much wider range of masses than stellar mass BH. Each channel contributes differently to the distributions of the mass, spin, distance, and orbital characteristics of BBHs.

There are several processes common to most pathways through stellar evolution which affect the properties of the resultant BBH system. Examples include mass loss (Vink et al. 2001; Vink & de Koter 2005; Gräfener & Hamann 2008) and supernovae (O’Connor & Ott 2011; Fryer et al. 2012; Janka 2012; Ugliano et al. 2012; Ertl et al. 2016; Sukhbold et al. 2016). The mass of the compact object left after the supernova is directly related to its pre-supernova mass and the supernova mechanism itself. Metallicity has been shown (Kudritzki & Puls 2000; Vink et al. 2001; Brott et al. 2011) to have important effects on stellar mass loss through winds — line-driven winds are quenched in metal-poor progenitors, enabling large black holes to form through direct collapse or post-supernova mass fallback (Heger et al. 2003; Mapelli et al. 2009; Belczynski et al. 2010; Spera et al. 2015). This also, in turn, might suppress supernova kicks (Fryer et al. 2012) and hence enhance the number of binaries which are not disrupted.

Theoretical and phenomenological models of BBH formation are explored by population synthesis. This requires modelling not only of stellar evolution but also the influence of their evolutionary environments. For instance, isolated evolution in galactic fields requires prescriptions for binary interactions, such as common envelope physics, as well as mass transfer episodes (see reviews in Kalogera et al. 2007; Vanbeveren 2009; Postnov & Yungelson 2014), and more recently, the effects of rapid rotation de Mink et al. 2009; Mandel & de Mink 2016; Marchant et al. 2016. Meanwhile, BBH formation in dense stellar clusters (Ziosi et al. 2014; Rodriguez et al. 2015; Rodriguez et al. 2016a; Mapelli 2016; Askar et al. 2017; Banerjee 2017) is impacted primarily by dynamical interactions within the cluster (Fregeau 2004; Morscher et al. 2013), but also by cluster size and initial mass functions (Scheepmaker et al. 2007; Portegies Zwart et al. 2010; Kremer et al. 2019). GW observations provide an alternative to sharpen our understanding of those processes.

Electromagnetic observations and modeling of systems containing black holes have led to speculation about the existence of potential gaps in the black hole mass spectrum. Both gaps may be probed using data from current ground-based gravitational-wave interferometers, and as such, have been the target of parametric studies. At low masses, observations of X-ray binaries (XRB) combined via Bayesian population modeling (Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011b) suggest a minimum black hole mass well above the largest neutron star masses. While the existence and nature of this gap is still uncertain (Kreidberg et al. 2012), it is proposed to exist between the most massive neutron stars (Özel & Freire 2016; Freire et al. 2008; Margalit & Metzger 2017) ($2.1-2.5\mathit{M_{\odot}}$) and the lightest black holes $\sim 5\mathit{M_{\odot}}$. It is possible to constrain the existence of this lower mass gap with GW observations (Littenberg et al. 2015; Mandel et al. 2015; Kovetz et al. 2017; Mandel et al. 2017). In Section III, we find our current GW observations do not inform the upper edge of this gap, inferring a minimum mass on the primary black hole at $m_{\textrm{min}}{}\lesssim 9$ $\mathit{M_{\odot}}$. Our volumetric sensitivity to BBH systems with masses less than 5 $\mathit{M_{\odot}}$ is small enough that we expect (and observe) no events in the lower gap region. Thus, our ability to place constraints in this region is severely limited.

Recently, there have been claims of an upper cutoff in the BBH mass spectrum based on the first few LIGO detections (Fishbach & Holz 2017; Talbot & Thrane 2018; Wysocki et al. 2018; Bai et al. 2018; Roulet & Zaldarriaga 2019). This might be expected as a consequence of a different supernova type, called the (pulsational) pair-instability supernova (Heger & Woosley 2002; Belczynski et al. 2016b; Woosley 2017; Spera & Mapelli 2017; Marchant et al. 2018). Evolved stars with a Helium core mass $\gtrsim{}30\mathit{M_{\odot}}$ are expected to become unstable, because efficient pair production softens their equation of state. For Helium core mass $\sim{}30$ – $64\mathit{M_{\odot}}$, the star undergoes a sequence of pulsations, losing mass until stability is reestablished (Woosley et al. 2007). The enhanced mass loss during pulsational pair instability is expected to affect the final collapse of the star, leading to smaller black hole masses. The fate of a star with He core mass $\sim{}64$ – $135\mathit{M_{\odot}}$ is more dramatic: the entire star is disrupted by a pair instability supernova, leaving no remnant (Fowler & Hoyle 1964; Barkat et al. 1967; Rakavy & Shaviv 1967). From the combination of pair instability and pulsational pair instability, it is expected that pair-instability supernovae should leave no black hole remnants between $\sim 50$ – $150\mathit{M_{\odot}}$ because the progenitor star is partially or entirely disrupted by the explosion. It is also possible that contributions from the merger of previous merger products — second generation mergers (O’Leary et al. 2016; Gerosa & Berti 2017; Fishbach et al. 2017; Rodriguez et al. 2018b) — could occupy this gap. Primordial BHs could also span numerous decades of the mass spectrum (Georg & Watson 2017), but their number density in either mass gap is dependent on the behavior of fluctuations in the early Universe (Byrnes et al. 2018). Nonetheless, consistent with prior work, we find that all our mass models have almost no merging black holes above $\sim 45$ $\mathit{M_{\odot}}$.

Observational constraints on the BBH merger rate (Abbott et al. 2016c; Abbott et al. 2018) generally assume a rate density which is uniform in the comoving volume. As first shown in Fishbach et al. 2018, it is also possible to search for redshift evolution in the rate density using current data. Different redshift-dependent evolutionary behavior is possible (Dominik et al. 2013; Mandel & de Mink 2016; Rodriguez et al. 2016a; Mapelli et al. 2017; Rodriguez & Loeb 2018) with different environments and stellar evolution scenarios (O’Shaughnessy et al. 2010; Belczynski et al. 2016a). For instance, theoretical models of isolated evolution through common envelope lead to a distribution of times to merger $p(t_{{\rm GW}})\propto t_{{\rm GW}}^{-1}$ (Dominik et al. 2012; Belczynski et al. 2016a). This would imply that many isolated binaries will coalesce near their formation redshift and produce a BBH merger rate that approximately tracks the star formation rate, peaking near $z\sim 2$. We find in Section IV that the current sample of BBH mergers does not provide enough information to confidently constrain any but the most extreme models. While we place more posterior mass on merger rates that increase with increasing redshift than those that decrease, the scenario of a uniform rate in comoving volume is comfortably within our constraints.

Black hole spin measurements also provide a powerful tool to discriminate between different channels of BBH formation (Mandel & O’Shaughnessy 2010; Abbott et al. 2016d; Rodriguez et al. 2016c; Vitale et al. 2017; Gerosa & Berti 2017; Farr et al. 2017; Farr et al. 2018; Gerosa et al. 2018). For example, BBHs formed in a dynamic environment will have no preferred direction for alignment, producing isotropically oriented spins (Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Mandel & O’Shaughnessy 2010; Rodriguez et al. 2015; Rodriguez et al. 2016c; Stone et al. 2017). However, some evidence has been presented for correlation in spin direction due to the natal environment of the progenitor stars within the cluster (Corsaro et al. 2017). In contrast, isolated binaries are expected to preferentially produce mergers with alignment between the spins of the constituent black holes and the orbital angular momentum of the system (Tutukov & Yungelson 1993; Kalogera 2000; Grandclément et al. 2004; Belczynski et al. 2016a; Rodriguez et al. 2016c; Mandel & de Mink 2016; Marchant et al. 2016; Stevenson et al. 2017b; O’Shaughnessy et al. 2017; Gerosa et al. 2018). Other effects occurring in stellar systems like hierarchical triples could also produce a weak preference for certain spin-orbit misalignments (Rodriguez & Antonini 2018). All of our parameterized models point to preferences against high spin magnitudes when the spin tilts are aligned with the orbital angular momentum. In Section V, we find that the dimensionless spin magnitude inference prefers distributions which decline as the spin magnitude increases from zero, but our ability to distinguish between assumed distributions of spin orientation is very limited.

GW170817, the first binary neutron star merger observed through GW emission (Abbott et al. 2017b), was detected by GW observatories and associated with a short GRB (Abbott et al. 2017c)) in August of 2017. A subsequent post-merger transient (AT 2017gfo) was observed across the electromagnetic spectrum, from radio (Alexander et al. 2017), NIR/optical (Coulter et al. 2017; Soares-Santos et al. 2017; Chornock et al. 2017; Cowperthwaite et al. 2017; Nicholl et al. 2017; Pian et al. 2017), to X-ray (Troja et al. 2017; Margutti et al. 2017) and $\gamma$-ray (Abbott et al. 2017c; Goldstein et al. 2017; Savchenko et al. 2017). Unfortunately, with only one confident detection, it is not yet possible to infer details of binary neutron star populations more than to note that the gravitational-wave measurement is mostly compatible with the observed Galactic population (Özel et al. 2012). However, if GW170817 did form a black hole, it would also occupy the lower mass gap described previously.

We structure the paper as follows. First, notation and models are established in Section II. Section III describes our modeling of the black hole mass distribution, followed by rate distributions and evolution in Section IV. The black hole spin magnitude and orientation distributions are discussed in Section V. We conclude in Section VI. Studies of various systematics are presented in Appendix A. In Appendix B we present additional studies of spin distributions with model selection for a number of zero-parameter spin models and mixtures of spin orientations. To motivate and enable more detailed studies, we have established a repository of our samples and other derived products The data release for this work can be found at https://dcc.ligo.org/LIGO-P1800324/public..

II Data, Notation, and Models

In this work, we analyze the population of 10 BBH merger events confidently identified in the first and second observing run (O1 and O2) (Abbott et al. 2018). We do not include marginal detections, but these likely have a minimal impact our conclusions here (Gaebel et al. 2019). Ordered roughly from smallest to most massive by source-frame chirp mass, the mergers considered in this paper are GW170608, GW151226, GW151012, GW170104, GW170814, GW170809, GW170818, GW150914, GW170823, and GW170729.

The individual properties of those 10 sources were inferred using a Bayesian framework, with results summarized in Abbott et al. 2018. For BBH systems, two waveform models have been used, both calibrated to numerical relativity simulations and incorporating spin effects, albeit differently: IMRPhenomPv2 (Hannam et al. 2014; Husa et al. 2016; Khan et al. 2016), which includes an effective representation (Schmidt et al. 2015) of precession effects, and SEOBNRv3 (Pan et al. 2014; Taracchini et al. 2014; Babak et al. 2017), which incorporates all spin degrees of freedom. The results presented in this work use IMRPhenomPv2; a discussion of potential systematic biases in our inference are discussed in Appendix A. We also refer to Appendix B in (Abbott et al. 2018) for more details on comparisons between those two waveform families.

To assess the stability of our results to statistical effects and systematic error we focus on one modestly exceptional event. Both GW151226 and GW170729 exhibit evidence for measurable black hole spin, but GW170729 in particular is an outlier by several other metrics as well. In addition to spins, it is also more massive and more distant than any of the other events in the catalog. All events used in the population analysis have confident probabilities of astrophysical origin, but GW170729 is the least significant, having the smallest odds ratio of astrophysical versus noise origin (Abbott et al. 2018). As we describe in Sections III and IV, this event has an impact on our inferred merger rate versus both mass and redshift. To demonstrate the robustness of our result, we present these analyses twice: once using every event, and again omitting GW170729.

A coalescing compact binary in a quasi-circular orbit can be completely characterized by its eight intrinsic parameters, namely its component masses $m_{i}$ and spins $\bm{S}_{i}$, and its seven extrinsic parameters: right ascension, declination, luminosity distance, coalescence time, and three Euler angles characterizing its orientation (e.g., inclination, orbital phase, and polarization). Binary eccentricity is also a potentially observable quantity in BBH mergers, with several channels having imprints on eccentricity distributions, e.g. (Quinlan & Shapiro 1987; Kocsis & Levin 2012; Samsing et al. 2014; Fragione et al. 2018; Rodriguez et al. 2018a). However, our ability to parameterize (Huerta et al. 2014; Huerta et al. 2017; Klein et al. 2018; Hinder et al. 2018) and measure (Coughlin et al. 2015; Abbott et al. 2016d; Abbott et al. 2017d; Lower et al. 2018) eccentricity is an area of active development. For low to moderate eccentricity at formation, binaries are expected to circularize (Peters 1964; Hinder et al. 2008) before entering the bandwidth of ground-based GW interferometers. We therefore assume zero eccentricity in our models.

In this work, we define the mass ratio as $q=m_{2}/m_{1}$ where $m_{1}\geq m_{2}$. The frequency of gravitational wave emission is directly related to the component masses. However, due to the expansion of spacetime as the gravitational wave is propagating, the frequencies measured by the instrument are redshifted relative to those emitted at the source (Thorne 1983). We capture these effects by distinguishing between masses as they would be measured in the source frame, denoted as above, and the redshifted masses, $(1+z)m_{i}$, which are measured in the detector frame. Meanwhile, the amplitude of the wave scales inversely with the luminosity distance (Misner et al. 1973). We use the GW measurement of the luminosity distance to obtain the cosmological redshift and therefore convert between detector-frame and source-frame masses. We assume a fixed Planck 2015 (Planck Collaboration et al. 2016) cosmology throughout to convert between a source’s luminosity distance and its redshift (Hogg 1999).

We characterize black hole spins using the dimensionless spin parameter $\bm{\chi}_{i}=\bm{S}_{i}/m_{i}^{2}$. Of particular interest are the magnitude of the dimensionless spin, $a_{i}=|{\bm{\chi}}_{i}|$, and the tilt angle with respect to the orbital angular momentum, $\hat{\bm{L}}$, given by $\cos t_{i}=\hat{\bm{L}}\cdot\hat{{\bm{\chi}}}_{i}$. We also define an overall effective spin, $\chi_{\textrm{eff}}{}$ (Damour 2001; Racine 2008; Ajith et al. 2011), which is a combination of the individual spin components along to orbital angular momentum:

$\chi_{\textrm{eff}}$ is approximately proportional to the lowest order contribution to the GW waveform phase that contains spin for systems with similar masses. Additionally, $\chi_{\textrm{eff}}$ is conserved throughout the binary evolution to high accuracy (Racine 2008; Gerosa et al. 2015).

II.2 Model Features

The current sample is not sufficient to allow for a high-fidelity comparison with models (e.g., population synthesis) which include more detailed descriptions of stellar evolution and environmental influences. As such, we adopt the union of the parameterizations presented in Talbot & Thrane 2017; Fishbach & Holz 2017; Wysocki et al. 2018; Talbot & Thrane 2018; Fishbach et al. 2018. This allows for better facilitation of comparison between models, and the ability to vary the subsets of parameters influencing the mass and spin distributions while leaving others fixed.

The general model family has 8 parameters to characterize the mass model; 3 to characterize each black hole’s spin distribution; one parameter describing the local merger rate, $\mathcal{R}_{0}$; and one parameter characterizing redshift dependence. We refer to the set of these population parameters as $\theta$. All of the population parameters introduced in this section are summarised in Table 1.

II.3 Parameterized Mass Models

The power-law distribution considered previously (Abbott et al. 2016c; Abbott et al. 2017e) modeled the BBH primary mass distribution as a one-parameter power-law, with fixed limits on the minimum and maximum allowed black hole mass. With our sample of ten binaries, we extend this analysis by considering three increasingly complex models for the distribution of black hole masses. The first extension, Model A (derived from Fishbach & Holz 2017; Wysocki et al. 2018), allows the maximum black hole mass $m_{\textrm{max}}{}$ and the power-law index $\alpha$ to vary. In Model B (derived from Kovetz et al. 2017; Fishbach & Holz 2017; Talbot & Thrane 2018) the minimum black hole mass $m_{\textrm{min}}{}$ and the mass ratio power-law index $\beta_{q}$ are also free parameters. However, the priors on Model B and C enforce a minimum of 5 $\mathit{M_{\odot}}$ on $m_{\textrm{min}}$ — see Table 2. Explicitly, the mass distribution in Model A and Model B takes the form

where $C\left(m_{1}\right)$ is chosen so that the marginal distribution is a power law in $m_{1}$: $p\left(m_{1}|m_{\textrm{min}}{},m_{\textrm{max}}{},\alpha,\beta_{q}\right)=m_{1}^{-\alpha}$.

Model A fixes $m_{\textrm{min}}{}=5\mathit{M_{\odot}}$ and $\beta_{q}=0$, whereas Model B fits for all four parameters. Equation 2 implies that the conditional mass ratio distribution is a power-law with $p(q\mid m_{1})\propto q^{\beta_{q}}$. When $\beta_{q}=0$, $C(m_{1})\propto 1/(m_{1}-m_{\textrm{min}}{})$, as assumed in Abbott et al. 2016c; Abbott et al. 2017e.

Model C (from Talbot & Thrane 2018) further builds upon the mass distribution in Equation 2 by allowing for a second, Gaussian component at high mass, as well as introducing smoothing scales $\delta m$, which taper the hard edges of the low- and high-mass cutoffs of the primary and secondary mass power-law. The second Gaussian component is designed to capture a possible build-up of high-mass black holes created from pulsational pair instability supernovae. The tapered low-mass smoothing reflects the fact that parameters such as metallicity probably blur the edge of the lower mass gap, if it exists. Model C therefore introduces four additional model parameters, the mean, $\mu_{m}$, and standard deviation, $\sigma_{m}$, of the Gaussian component, $\lambda_{m}$, the fraction of primary black holes in this Gaussian component, and $\delta m$ the smoothing scale at the low mass end of the distribution.

The factors $A$, $B$, and $C$ ensure each of the power-law component, Gaussian component, and mass ratio distributions are correctly normalized. $S$ is a smoothing function which rises from zero at $m_{\min}$ to one at $m_{\min}+\delta m$ as defined in Talbot & Thrane 2018. $\Theta$ is the Heaviside step function. Models A, B, and C are displayed with a selection of parameters for demonstration purposes in the left panels of Figure 1.

II.4 Parameterized Spin Models

The black hole spin distribution is decomposed into independent models of spin magnitudes, $a$, and orientations, $t$. For simplicity and lacking compelling evidence to the contrary, we assume both black hole spin magnitudes in a binary, $a_{i}$, are drawn from a beta distribution (Wysocki et al. 2018):

To describe the spin orientation, we assume that the tilt angles between each black hole spin and the orbital angular momentum, $t_{i}$, are drawn from a mixture of two distributions: an isotropic component, and a preferentially aligned component, represented by a truncated Gaussian distribution in $\cos t_{i}$ peaked at $\cos t_{i}=1$ (Talbot & Thrane 2017)

We choose to parameterize the cosine of the tilt angles, rather than the angles themselves. This choice prompts the selection of a Gaussian (or uniform) model, rather than a wrapped distribution which would be more appropriate for an angular variable. An example of the Mixture distribution is displayed in the lower right panel of Figure 1.

The parameter $\zeta$ denotes the fraction of binaries which are preferentially aligned with the orbital angular momentum; $\zeta=1$ implies all black hole spins are preferentially aligned and $\zeta=0$ is an isotropic distribution of spin orientations. The typical degree of spin misalignment is represented by the $\sigma_{i}$. For spin orientations we explore two parameterized families of models:

The Gaussian model is motivated by formation in isolated binary evolution, with significant natal misalignment, while the mixture scenarios allow for an arbitrary combination of this scenario and randomly oriented spins, which arise naturally in dynamical formation.

II.5 Redshift Evolution Models

The previous two subsections described the probability distributions of intrinsic parameters $p(\xi)$ (i.e. masses and spins) that characterize the population of BBHs. In addition, we also measure the value of one extrinsic parameter of the population: the overall merger rate density $R$. The models described in the previous two subsections assume that the distribution of intrinsic parameters is independent of cosmological redshift $z$, at least over the redshift range accessible to the LIGO and Virgo interferometers during the first two observing runs ($z\lesssim 1$). However, we consider an additional model in which the overall event rate evolves with redshift. We follow Fishbach et al. 2018 by parameterizing the evolving merger rate density $R(z)$ in the comoving frame by

where $R_{0}$ is the rate density at $z=0$. In this model, $\lambda=0$ corresponds to a merger rate density that is uniform in comoving volume and source-frame time, while $\lambda\sim 3$ corresponds to a merger rate that approximately follows the star-formation rate in the redshift range relevant to the detections in O1 and O2 (Madau & Dickinson 2014). Various BBH formation channels predict different merger rate histories, ranging from rate densities that will peak in the future ($\lambda<0$) to rate densities that peak earlier than the star-formation rate ($\lambda\gtrsim 3$). These depend on the formation rate history and the distribution of delay times between formation and redshift. In cases where we do not explicitly write the event rate density as $R(z)$, it is assumed that the rate density $R$ is constant in comoving volume and source-frame time.

The general model family, including the distributions of masses, spins and merger redshift, is therefore given by the distribution

where $t$ is the time in the source-frame, so that Eq. 8 can be written equivalently in terms of the merger rate density:

II.6 Hierarchical Population Model

We perform a hierarchical Bayesian analysis, accounting for measurement uncertainty and selection effects (Loredo 2004; Abbott et al. 2016c; Wysocki et al. 2018; Fishbach et al. 2018; Mandel et al. 2018; Mortlock et al. 2018). We model the occurrence rate of events through a Poisson process with a mean dependent on the parameter distribution of the binaries While this assumption is embedded (Farr et al. 2015) in the selection of events used in this work, studies of event count per time do not show significant evidence for deviations from Poissonian statistics (Abbott et al. 2018).. The likelihood of the observed GW data given the population hyperparameters $\theta$ that describe the general astrophysical distribution, $dN/d\xi dz$, is given by the inhomogeneous Poisson likelihood:

In order to calculate the expected number of detections $\mu(\theta)$, we must understand the selection effects of our detectors. The sensitivity of GW detectors is a strong function of the binary masses and distance, and also varies with spin. For any binary, we define the sensitive spacetime volume $VT(\xi)$ of a network with a given sensitivity to be

where $p(\xi|\theta)$ describes the underlying distribution of the intrinsic parameters. We performed large scale simulation runs wherein the spacetime volume in the above equation is estimated by Monte-Carlo integration (Tiwari 2018) — these runs are restricted to have no BH less massive than 5 $\mathit{M_{\odot}}$. We then use a semi-analytic prescription, calibrated to the simulation results, to derive the $\langle VT\rangle{}_{\theta}$ for specific hyper-parameters.

Allowing the merger rate to evolve with redshift, the expected number of detections is given by

If the merger rate does not evolve with redshift, i.e., $R(z)=R_{0}$, this reduces to $\mu(\theta)=R_{0}\langle VT\rangle{}_{\theta}$.

We note that the hyperparameter likelihood given by Eq. II.6 reduces to the likelihood used in the O1 mass distribution analysis (Abbott et al. 2016c, Eq. D10 of), which fit only for the shape, not the rate / normalization of the mass distribution, if one marginalizes over the rate parameter with a flat-in-log prior $p(R_{0})\propto 1/R_{0}$ (Fishbach et al. 2018; Mandel et al. 2018). For consistency with previous analyses, we adopt a flat-in-log prior on the rate parameter throughout this work.

II.7 Statistical Framework and Prior Choices

In practice, we sample the likelihood $\mathcal{L}(d_{n}|\xi,z)$ using the parameter estimation pipeline LALInference (Veitch et al. 2015). Since LALInference gives us a set of posterior samples for each event, we first divide out the priors used in the individual-event analyses before applying Eq. II.6 (Hogg et al. 2010; Mandel 2010) (see Appendix C).

Where not fixed, we adopt uniform priors on population parameters describing the models. Unless otherwise noted, for the event rate distribution we use a log-uniform distribution in $R_{0}$, bounded between $[10^{-1},10^{3}]$. While this is a different form than the priors adopted in Abbott et al. 2018, we note that similar results are obtained on the rates (see Sec. IV), indicating that the choice of prior does not strongly influence the posterior distributions. We provide specific limits on all priors when the priors for a given model are introduced. Unless otherwise stated all posterior credible intervals are 90% intervals, symmetric in the quantiles around the median. The MCMC-based analyses presented in this work have approximately $10^{4}$ effective samples, after thinning by their autocorrelation time.

The normalization factor of the posterior density in Bayes’ theorem is the evidence — it is the probability of the data given the model. We are interested in the preferences of the data for one model versus another. This preference is encoded in the Bayes factor, or the ratio of evidences. The odds ratio is the Bayes factor multiplied by their ratio of the model prior probabilities. In all cases presented here, the prior model probabilities are assumed to be equal, and odds ratios are equivalent to Bayes factors.

We often present the posterior population distribution (PPD) of various quantities. The PPD is the expected distribution of new mergers conditioned on previously obtained observations. It integrates the distribution of values (e.g., $\xi$, such as the masses and spins) conditioned on the model parameters (e.g, the power law index) over the posteriors obtained for the model parameters:

It is a predictor for future merger values $\xi_{\textrm{new}}$ given observed data $\xi_{\textrm{observed}}$ and factors in the uncertainties imposed by the posterior on the model parameters. Note that the PPD does not incorporate the detector sensitivity, and therefore is not a straightforward predictor of the properties of future observed mergers.

III The Mass Distribution

For context, Figure 4 in Abbott et al. 2018 illustrates the inferred masses for all of the significant BBH observations identified in our GW surveys in O1 and O2. Despite at least moderate sensitivity to total masses between 0.1 – 500 $\mathit{M_{\odot}}$, current observations occupy only a portion of the binary mass parameter space. Notably, we have not yet observed a pair of very massive (e.g., 100 $\mathit{M_{\odot}}$) black holes, a binary which is bounded away from equal mass in its posterior, or a binary with a component mass confidently below 5 $\mathit{M_{\odot}}$. In our survey, we also find a preponderance of observations at higher masses: six with significant posterior support above $30\mathit{M_{\odot}}$. In this section, we attempt to reconstruct the binary black hole merger rate as a function of the component masses using parameterized models. Table 2 summarizes the mass models adopted from Section II.3 and the prior distributions for each of the parameters in those models. We present results for three increasingly general mass and spin models, the most complex of which ranges over the full set of model parameters in Section II with the exception of rate dependence of rate on redshift. The interdependence of the mass and redshift distribution is explored more fully in Section IV.

Figure 2 shows our updated inference for the compact binary primary mass $m_{1}$ and mass ratio $q$ distributions for several increasingly general population models. In addition to inferring the mass distribution, all of these calculations self-consistently marginalize over the parameterized spin distribution presented in Section V and the merger rate. Figures 3 and 4 show the posterior distribution on selected model hyperparameters.

All models feature a parameter, $m_{\textrm{max}}$, which defines a cutoff of the power law. However, the interpretation of that parameter within Model C is not a straightforward comparison with Models A and B, due to the presence of the Gaussian component at high mass and the large value of the power-law spectral index. Instead, to compare those two features, we compute the 99th percentile of the mass distribution inferred from the model PPDs (see Equation 15). Model A obtains ${44.0}$ $M_{\odot}$, Model B obtains ${41.8}$ $M_{\odot}$, and Model C obtains $41.8$ $M_{\odot}$. Therefore, all models self-consistently infer a dearth of black holes above $\sim 45$ $\mathit{M_{\odot}}$. This is determined by the lower limit for the mass of the most massive black hole in the sample because $m_{\textrm{max}}$ can be no smaller than this value. Similarly, the models which allow $m_{\textrm{min}}{}$ to vary (B and C) disfavor populations with $m_{\textrm{min}}{}$ above $\simeq 9M_{\odot}$. This parameter is close to the largest allowed mass for the least massive black hole in the sample, for similar reasons.

The lower limits we place on $m_{\textrm{min}}{}$ are dominated by our prior choices that constrain $m_{\textrm{min}}{}\in\,\mathit{M_{\odot}}$ (see Table 2). For example, in Figure 3, the posterior on $m_{\textrm{min}}{}$ becomes flat as $m_{\textrm{min}}{}$ approaches the prior boundary at $5\,\mathit{M_{\odot}}$. Given current sensitivities, this is to be expected (Littenberg et al. 2015; Mandel et al. 2015). In the inspiral-dominated regime, the sensitive time-volume scales as $VT\sim m^{15/6}$ (Finn & Chernoff 1993); extending our inferred mass distributions and merger rates into the possible lower black hole mass gap from $3$–$5\,\mathit{M_{\odot}}$ (Özel et al. 2010; Farr et al. 2011b; Kreidberg et al. 2012) yields an expected number of detected BBH mergers $\lesssim 1$. Thus, we are unable to place meaningful constraints on the presence or absence of a mass gap at low black hole mass.

Models B and C also allow the distribution of mass ratios to vary according to $\beta_{q}$. In these cases the inferred mass-ratio distribution favors comparable-mass binaries (i.e., distributions with most support near $q\simeq 1$), see panel two of Figure 2. Within the context of our parameterization, we find $\beta_{q}={{6.9}^{+4.6}_{-5.7}}$ for Model B and $\beta_{q}={{4.5}^{+6.6}_{-5.2}}$ for Model C. These values are consistent with each other and are bounded above zero at 95% confidence, thus implying that the mass ratio distribution is nearly flat or declining with more extreme mass ratios. The posterior on $\beta_{q}$ returns the prior for $\beta_{q}\gtrsim 4$. Thus, we cannot say much about the relative likelihood of asymmetric binaries, beyond their overall rarity.

III.2 Comparison with Theoretical and Observational Models

Previous modeling of the primary mass distribution with a power law distribution (Abbott et al. 2016c) was last updated with the discovery of GW170104 (Abbott et al. 2017e). This analysis measured spectral index of the the power law to be $\alpha={2.3}^{+1.3}_{-1.4}$ at 90% confidence assuming a minimum black hole mass of 5 $\mathit{M_{\odot}}$ and maximum total mass of 100 $\mathit{M_{\odot}}$. None of our models directly emulate this one, but Model A is the closest analog. When allowing $m_{\textrm{max}}$ to vary, 100 $\mathit{M_{\odot}}$ is strongly disfavored. As a consequence of the lower $m_{\textrm{max}}$, the power law index inferred is also shallower than previously obtained (Fishbach & Holz 2017), but remains consistent with the previous distribution.

In Figure 5, we highlight the two mass gaps predicted by models of stellar evolution: the first gap between $\sim{}2$ and $\sim{}5$ M⊙ and the second between $\sim{}50$ M⊙ and $\sim{}150$ M⊙, compared against the observed black holes. A set of tracks (Spera & Mapelli 2017) relating the progenitor mass and compact object is also shown for reference purposes. The tracks are subject to many uncertainties in stellar and binary evolution, and only serve as representative examples. We discuss some of those uncertainties in the context of our results below.

The minimum mass of a black hole and the existence of a mass gap between neutron stars and black holes (lower gray shaded area, right panel of Figure 5) are currently debated. Claims (Özel et al. 2010; Farr et al. 2011b) of the existence of a mass gap between the heaviest neutron stars ($\sim{}2$ M⊙) and the lightest black holes ($\sim{}5$ M⊙) are based on the sample of about a dozen X-ray binaries with dynamical mass measurements. However, Kreidberg et al. 2012 suggested that the dearth of observed black hole masses in the gap could be due to a systematic offset in mass measurements. Moreover, only a subset of theoretical models (e.g., the “rapid” model in Fryer et al. 2012) reproduce this gap in stellar modelling. We can see in Figure 5 that none of the observed binaries sit in this gap, but the sample is not sufficient to definitively confirm or refute the existence of this mass gap.

From the first six announced BBH detections, Fishbach & Holz 2017 argued that there is evidence for missing black holes with mass greater than $\gtrsim{}40$ M⊙. The existence of this second mass gap — see the upper grey shaded area in the right panel of Figure 5 between $\sim{}50$ M⊙ and $\sim{}150$ M⊙ — has been further explored by Talbot & Thrane 2018; Wysocki et al. 2018; Bai et al. 2018; Roulet & Zaldarriaga 2019. This gap might arise from the combined effect of pulsational pair instability (Barkat et al. 1967; Heger et al. 2003; Woosley et al. 2007; Woosley 2017) and pair instability (Fowler & Hoyle 1964; Ober et al. 1983; Bond et al. 1984) supernovae. Uncertainties in stellar evolution models (e.g. stellar winds, rotation) and in the treatment of the final outcomes of (pulsational) pair instability lead to a range of possible low-mass edges for the upper mass gap as well as the shape and abundance in a putative build-up. Predictions for the maximum mass of black holes born after pulsational pair-instability supernovae are $\sim{}50\mathit{M_{\odot}}$ (Belczynski et al. 2016b; Spera & Mapelli 2017). Our inferred maximum mass is consistent with these predictions.

IV Merger Rates and Evolution with Redshift

As illustrated in previous work (Abbott et al. 2016e; Abbott et al. 2018; Fishbach & Holz 2017; Wysocki et al. 2018; Fishbach et al. 2018), the inferred binary black hole merger rate depends on and correlates with our assumptions about their intrinsic mass (and to a lesser extent, spin) distribution. In the most recent catalog of GW BBH events (Abbott et al. 2018), we infer the overall BBH merger rate for two fixed-parameter populations. The first of these populations follows the power-law model given by Equation 2 with $\alpha=2.3$, $\beta_{q}=0$, $m_{\textrm{min}}{}=5\mathit{M_{\odot}}$, and $m_{\textrm{max}}{}=50\mathit{M_{\odot}}$. The second population follows a distribution in which both black hole masses are independently drawn from a flat-in-log distribution:

subject to the same mass cutoffs $5\mathit{M_{\odot}}<m_{2}<m_{1}<50\mathit{M_{\odot}}$ as the fixed power-law population. Both the power-law and flat-in-log populations assume an isotropic and uniform-magnitude spin distribution ($\alpha_{a}=\beta_{a}=1$). These two fixed-parameter populations are used to estimate the population-averaged sensitive volume $\langle VT\rangle{}$ with a Monte-Carlo injection campaign as described in Abbott et al. 2018, with each population corresponding to a different $\langle VT\rangle{}$ because of the strong correlation between the mass spectrum and the sensitive volume. Under the assumption of a constant-in-redshift rate density, these $\langle VT\rangle{}$ estimates yield two different estimates of the rate: $57^{+40}_{-25}$ Gpc-3 yr-1for the $\alpha=2.3$ population, and $19^{+13}_{-8.2}$ Gpc-3 yr-1for the flat-in-log population (90% credibility; combining the rate posteriors from the two analysis pipelines).

In these calculations, we first maintain the assumption in Abbott et al. 2018 that the merger rate is uniform in comoving volume and source-frame time, as discussed in Section II. We then relax this assumption and consider a merger rate that evolves in redshift according to Equation 7, fitting the mass distribution jointly with the rate density as a function of redshift.

We first consider the case of a uniform in volume merger rate, and examine the effects of fitting the rate jointly with the distribution of masses and spins. The first column in Figures 3 and 4 shows the results of self-consistently determining the rate using the models for the mass and spin distribution described in the previous two sections.

IV.2 Evolution of the Merger Rate with Redshift

As discussed in the introduction, most formation channels predict some evolution of the merger rate with redshift, due to factors including the star-formation rate, time-delay distribution, metallicity evolution, and globular cluster formation rate (Dominik et al. 2013; Belczynski et al. 2016a; Mandel & de Mink 2016; Rodriguez & Loeb 2018). Therefore, in this section, we allow the merger rate to evolve with redshift, and infer the redshift evolution jointly with the mass distribution. For simplicity, we adopt the two-parameter Model A for the mass distribution and fix spins to zero for this analysis. As discussed in Section III, the additional mass and spin degrees of freedom have only a weak effect on the inferred merger rate. We assume the redshift evolution model given by Equation 7. Because massive binaries are detectable at higher redshifts, the observed redshift evolution correlates with the observed mass distribution of the population, and so we must fit them simultaneously. However, as in Fishbach et al. 2018, we assume that the underlying mass distribution does not vary with redshift. We therefore fit the joint mass-redshift distribution according to the model:

Note that this model assumes that the merger rate density increases or decreases monotonically with redshift over the sensitive range $z<1$. If the merger rate follows the star formation rate, we expect the rate to peak around $z\sim 2$, which is currently far beyond the horizon redshift for BBH detections.

Marginalizing over the two mass distribution parameters and the redshift-evolution parameter, the merger rate density is consistent with being constant in redshift ($\lambda=0$), and in particular, it is consistent with the rate estimates recovered under the different mass distribution models in subsection IV.1 above. However, we find a preference for a merger rate density that increases at higher redshift ($\lambda\geq 0$) with probability $0.93$. This implies that models that predict a constant, or slightly decreasing merger rate with redshift, such as certain models of primordial black holes (Mandic et al. 2016), are disfavored. This preference for a merger rate that increases with increasing redshift becomes less significant when GW170729 is excluded from the analysis, because this event likely merged at redshift $z\gtrsim 0.5$, close to the O1-O2 detection horizon. Although GW170729 shifts the posterior towards larger values of $\lambda$, implying a stronger redshift evolution of the merger rate, the posterior remains well within the uncertainties inferred from the remaining nine BBHs. When including GW170729 in the analysis, we find $\lambda=8.4^{+9.6}_{-9.5}$ at 90% credibility, compared to $\lambda=2.3^{+9.9}_{-10.9}$ when excluding GW170729 from the analysis. With only 10 BBH detections so far, the wide range of possible values for $\lambda$ is consistent with most astrophysical formation channels. The precision of this measurement will improve as we accumulate more detections in future observing runs and may enable us to discriminate between different formation rate histories or time-delay distributions (Sathyaprakash et al. 2012; Van Den Broeck 2014; Fishbach et al. 2018).

V The Spin Distribution

The GW signal depends on spins in a complicated way, but at leading order, and in the regime we are interested in here, some combinations of parameters have more impact on our inferences than others, and thus are measurable. One such parameter is $\chi_{\textrm{eff}}$. For binaries which are near equal mass, we can see from Equation 1 that only when black hole spins are high and aligned with the orbital angular momentum $\chi_{\textrm{eff}}$ will be measurably greater than zero. Figure 5 in Abbott et al. 2018 illustrates the inferred $\chi_{\textrm{eff}}$ spin distributions for all of the BBHs identified in our GW surveys in O1 and O2. Only GW170729 and GW151226 show significant evidence for positive $\chi_{\textrm{eff}}{}$; the rest of the posteriors cluster around $\chi_{\textrm{eff}}{}=0$.

Despite these degeneracies, several tests have been proposed to use spins to constrain BBH formation channels (Vitale et al. 2017; Farr et al. 2017; Farr et al. 2018; Stevenson et al. 2017a; Talbot & Thrane 2017; Gerosa & Berti 2017; Wysocki et al. 2018; Gerosa et al. 2018). Drawing upon these methods, we now seek to estimate the black hole spin magnitude and misalignment distributions, under different assumptions regarding isotropy or alignment.

We examine here the individual spin magnitudes and tilt distributions. Throughout this section, when referring to the parametric models, we also allow the merger rate and population parameters describing the most general mass model to vary (Model C, see Table 2). Changing the parameterization of the mass model does not significantly change our inferences about the spin distribution. However, to account for degeneracies between mass and spin that grow increasingly significant for longer, low-mass signals (Baird et al. 2013), we must consistently model the mass and spin distributions together. See Table 6 for a summary of the models and priors used in this Section.

We also compute the posterior distribution for the magnitude of black hole spins from $\chi_{\textrm{eff}}{}$ measurements by modeling the distribution of black hole spin magnitudes non-parametrically with five bins, assuming either an isotropic or perfectly aligned population following Farr et al. 2018. We show in the bottom panel of Figure 8 that under the perfectly aligned scenario there is preference for small black hole spin, inferring 90% of black holes to have spin magnitudes below $0.6^{+0.24}_{-0.28}$. However, when spins are assumed to be isotropic the distribution is relatively flat, with 90% of black hole spin magnitudes below $0.8^{+0.15}_{-0.24}$. Thus, the non-parametric analysis produces conclusions consistent with our parametric analyses described above. These conclusions are also reinforced by computing the Bayes factor for a set of fixed parameter models of spin magnitude and orientation in Appendix B. There we find that the very low spin magnitude model is preferred by a log Bayes factor of 1 or greater in most mass and spin orientation configurations tested (see Figure 13 and Table 8 for details).

Figure 9 shows the inferred distribution of the primary spin tilt for the more massive black hole. These results were obtained without including the effects of component spins on the detection probability: see Appendix A for further discussion. In the Gaussian model ($\zeta=1$), all black hole spin orientations are drawn from spin tilt distributions which are preferentially aligned and parameterized with $\sigma_{i}$. In that model, the $\sigma_{i}$ distributions do not differ appreciably from the their flat priors. As such, the inferred spin tilt distribution are influenced by large $\sigma_{i}$ and the result resembles an isotropic distribution. The Mixture distribution does not return a decisive measurement of the mixture fraction, obtaining $\zeta={{0.6}^{+0.4}_{-0.5}}$. Since the Gaussian model is a subset of the Mixture model, we can compare preferences via the Savage-Dickey ratio. The log Bayes factor for $\zeta=1$ is $\ln\textrm{BF}=0.15$, indicating virtually no preference for any particular orientation distribution. While we allow both black holes to have different typical misalignment, the inference on the second tilt is less informative than the primary. The inferred distribution for $\cos t_{2}$ is similar to $\cos t_{1}$, but also closer to the prior.

The mixture fraction distribution is also modelled with the fixed parameter models in Appendix B. The fixed magnitude distributions considered in Appendix B prefer isotropic to aligned, but the preference is weakened for distributions concentrated at lower spins. A few exceptions occur for the very low spin fixed mass ratio models, with aligned models being slightly preferred.

In general, we are not able to place strong constraints on the distribution of spin orientations. We elaborate in Appendix B.3 on how our black hole spin measurements are not yet informative enough to discern between isotropic and aligned orientation distribution via $\chi_{\textrm{eff}}$.

V.2 Interpretation of Spin Distributions

The spins of black holes are affected by a number of uncertain processes which occur during the evolution of the binary. As a consequence, the magnitude distribution is difficult to predict from theoretical models of these processes alone. While the spin of a black hole should be related to the rotation of the core of its progenitor star, the amount of spin which is lost during the final stages of the progenitor’s life is still highly uncertain. While we have modeled the spins independently, correlations from binary evolution and stellar collapse are possible (Belczynski et al. 2017; Gerosa et al. 2018; Qin et al. 2018; Postnov & Kuranov 2019; Arca Sedda & Benacquista 2019). The core rotational angular momentum before the supernova can be changed from the birth spin of the progenitor by several processes (Langer 2012; de Mink et al. 2013; Amaro-Seoane & Chen 2016). Examples include mass transfer (Shu & Lubow 1981; Packet 1981), or tidal interactions (Petrovic et al. 2005), as well as internal mixing of the stellar layers across the core-envelope boundary via magnetic torquing (Spruit 2002; Maeder & Meynet 2003) and gravity waves (Talon & Charbonnel 2005; Talon & Charbonnel 2008; Fuller et al. 2015). In principle, an off-center supernova explosion could also impart significant angular momentum and tilt the spin of the remnant into the collapsing star (Farr et al. 2011a).

Once a black hole is formed, however, changing the spin magnitude is more difficult due to limitations on mass accretion rates affecting how much a black hole can be spun up (Thorne 1974; Valsecchi et al. 2010; Wong et al. 2012; Qin et al. 2019). Once the binary black hole system is formed, the spin magnitudes do not change appreciably over the inspiral (Farr et al. 2014).

No BBH detected to date has a component with confidently high and aligned spin magnitude. The results in the previous section imply that black holes tend to be born with spin less than our PPD bound of $0.55$, or that another process (e.g., supernova kicks or dynamical processes involved in binary formation) induces tilts such that $\chi_{\textrm{eff}}$ is small.

The possibility of a spin magnitude distribution that peaks at low spins incurs a degeneracy between models that is not easily overcome: when the spin magnitudes are small enough models produce features which cannot be distinguished within observational uncertainties.

VI Discussion and Conclusions

We have presented a variety of estimates for the mass, spin, and redshift distributions of BBH, based on the observed sample of 10 BBH and generic phenomenological population models motivated by electromagnetic observations and theory. Some model independent features are evident from the observations. Notably, no binary black holes more massive than GW170729 have been observed to date, but several binaries have component masses likely between $20-40M_{\odot}$. No highly asymmetric (small $q$) system has been observed. Only two systems (GW151226 and GW170729) produce a $\chi_{\textrm{eff}}$ distribution which is confidently different from zero; conversely, most BH binaries are consistent with $\chi_{\textrm{eff}}$ near zero. These features drive our inferences about the mass and spin distribution.

Despite exploring a wide range of mass and spin distributions, we find the BBH merger rate density is $R={{64.0}^{+73.5}_{-33.0}}\,\mbox{Gpc}^{-3}\,\mbox{yr}^{-1}$ for Model A and is within $R={{53.2}^{+55.8}_{-28.2}}\,\mbox{Gpc}^{-3}\,\mbox{yr}^{-1}$ for Models B and C. This result is consistent with the fixed model assumptions reported in the combined O1 and O2 observational periods (Abbott et al. 2018). We find a significant reduction in the merger rate for binary black holes with primary masses larger than $\sim 45M_{\odot}$. We do not have enough sensitivity to binaries with a black hole mass less than 5 $\mathit{M_{\odot}}$ to be able to place meaningful constraints on the minimum mass of black holes. We find mild evidence that the mass distribution of coalescing black holes may not be a pure power law, instead being slightly better fit by a model including a broad gaussian distribution at high mass. We find the best-fitting models preferentially produce comparable-mass binaries (i.e., $\beta_{q}>0$ is preferred).

The mass models in this work supersede results from an older model from O1 which inferred only the power law index (Abbott et al. 2016c; Abbott et al. 2017e). That model found systematically larger values of $\alpha$ than its nearest counterpart in this work, Model A, because the older model used a fixed value for the minimum and maximum mass of 5 and 100 $\mathit{M_{\odot}}$, respectively. This extreme $m_{\textrm{max}}$ is highly disfavored by our current results, and so the older model is also disfavored. Moreover, volumetric sensitivity grows as a strong function of mass. The lack of detections near the older $m_{\textrm{max}}$ drives a preference for a much smaller maximum BH mass in the new models (Fishbach & Holz 2017). A reduced maximum mass is associated with a shallower power-law fit.

Inferring the redshift distribution is difficult with only a small sample of local events (Fishbach et al. 2018). We have constrained models with extreme variation over redshift, favoring instead those which are uniform in the comoving volume or have increasing merger rates with higher redshift. Many potential formation channels in the literature (Belczynski et al. 2014; Rodriguez et al. 2016b; Antonini & Rasio 2016; Mandel & de Mink 2016; Inayoshi et al. 2016; Mapelli et al. 2017; Bartos et al. 2017; Kruckow et al. 2018) produce event rates which are compatible with those from the previous observing runs (Abbott et al. 2018) and this work. It is, of course, plausible that several are contributing simultaneously, and no combination of mass, rate, or redshift dependence explored here rules out any of the channels proposed to date. The next generation of interferometers will allow for an exquisite probe into this dependence at large redshifts (Sathyaprakash et al. 2012; Van Den Broeck 2014; Vitale & Farr 2018).

We have modeled the spin distribution in several ways, forming inferences on the spin magnitude and tilt distributions. In all of our analysis, the evidence disfavors distributions with large spin components aligned (or nearly aligned) with the orbital angular momentum; specifically, we find that 90% of the spin magnitude PPD is smaller than $0.55$. We cannot significantly constrain the degree of spin-orbit misalignment in the population. However, regardless of the mass or assumed spin tilt distribution, there is a preference (demonstrated in Figure 8 and Appendix B) for distributions which emphasize lower spin magnitudes. Our inferences suggest 90% of coalescing black hole binaries are formed with $\chi_{\textrm{eff}}{}<{0.3}$. Low spins argue against so-called second generation mergers, where at least one of the components of the binary is a black hole formed from a previous merger (González et al. 2007; Berti et al. 2007) and possesses spins near 0.7 (Fishbach et al. 2017).

GW170729 is notable in several ways: it is the most massive, largest $\chi_{\textrm{eff}}$, and most distant redshift event detected so far. To quantify the impact it has on our results, where possible, we have presented model posteriors which reflect its presence in or exclusion from the event set. Many of our predictions are robust despite its extreme values — by far, and not unexpectedly, its influence is most significant in the distribution of $m_{\textrm{max}}$. It also impacts our conclusions about redshift evolution, where its absence flattens the inferred redshift evolution.

Recent modelling using only the first six released events (Wysocki et al. 2018; Roulet & Zaldarriaga 2019) have come to similar conclusions about low spin magnitudes and the shape of the power law distribution. The presence of an apparent upper limit to the merging BBH mass distribution was also observed after the first six released events (Fishbach & Holz 2017). An enhancement which will benefit these types of analyses in the future is a simultaneous fit of the astrophysical model and its parameters and noise background model (Gaebel et al. 2019).

Several studies have noted that population features (Mandel & O’Shaughnessy 2010; Stevenson et al. 2015; Fishbach & Holz 2017; Stevenson et al. 2017a; Zevin et al. 2017; Kovetz et al. 2017; Farr et al. 2017; Talbot & Thrane 2017; Fishbach et al. 2017; Gerosa & Berti 2017; Talbot & Thrane 2018; Farr et al. 2018; Barrett et al. 2018; Wysocki et al. 2018; Gerosa et al. 2018) and complementary physics (Abbott et al. 2016f; Zevin et al. 2017; Stevenson et al. 2017a; Chen et al. 2018) will be increasingly accessible as observations accumulate. Additional events will also permit the enhancement of the simple phenomenological models used in this work and comparison with modeling of astrophysical processes. Given the event merger rates estimated here and anticipated improvements in sensitivity (Abbott et al. 2018), hundreds of BBHs and tens of binary neutron stars are expected to be collected in the operational lifetime of second generation GW instruments. Thus, the inventory of BBH in the coming years will enable inquiries into astrophysics which were previously unobtainable.

References

Appendix A Systematics

In this section, we discuss the systematic uncertainties that affect our analysis, and show that they are subdominant to statistical uncertainties. We focus on two major sources of systematic uncertainty. The first of these is introduced by the waveform models that are used to extract the parameters of individual events, and the second is in the estimation of the detection efficiency.

A.2 Selection Effects and Sensitive Volume

In this subsection we detail the various assumptions and possible systematics that enter into our calculation of the detection efficiency. The detectability of a BBH merger in GWs depends on the distance and orientation of the binary along with its intrinsic parameters, especially its component masses. In order to model the underlying population and determine the BBH merger rate, we must properly model the mass, redshift and spin-dependent selection effects, and incorporate them into our population analysis according to Equation II.6. One way to infer the sensitivity of the detector network to a given population of BBH mergers is by carrying out large-scale simulations in which synthetic GW waveforms are injected into the detector data and subsequently searched for. The parameters of the injected waveforms can be drawn directly from the fixed population of interest, or alternatively, the injections can be placed to more broadly cover parameter space and reweighed to match the properties of the population (Tiwari 2018). Such injection campaigns were carried out in Abbott et al. 2018 to measure the total sensitive spacetime volume $\langle VT\rangle{}$ and the corresponding merger rate for two fixed-parameter populations (power-law and flat-in-log). However, it is computationally expensive to carry out an injection campaign that sufficiently covers the multi-dimensional population hyper-parameter space considered in this work. For this reason, for the parametric population studies in this work, we employ a semi-analytic method to estimate the fraction of found detections as a function of masses, spins and redshift (or equivalently, distance).

We therefore pursue two modifications to the raw semi-analytic calculation in order to reduce the bias in our sensitivity estimates and the resulting population estimates. We emphasize that these modifications do not noticeably affect the inferred shape of the population, e.g. the mass power-law slope, but do lead to different rate estimates, reflecting a systematic uncertainty in the inferred merger rate and its evolution with redshift that, given the small number of events and uncertainty in the phenomenological population models, remains subdominant to the statistical uncertainty. This is explicitly shown in the remainder of this section.

The top panel of Figure 11 shows the comparison between the raw semi-analytic $\langle VT\rangle{}$, the calibrated $\langle VT\rangle{}$, and the injection $\langle VT\rangle{}$ across the two-dimensional hyperparameter space of Model A for the mass distribution. We have repeated our mass distribution analysis with different choices of the $\langle VT\rangle{}$ calibration, and found that the effect on the shape of the mass distribution and the overall merger rate $\mathcal{R}$ are much smaller than the differences between Models A, B and C and the statistical errors associated with a small sample of 10 events.

As shown in Figure 11, the main effect of this calibration is to decrease $\langle VT\rangle{}$ by a factor of $\sim 1.6$. Over the relevant part of parameter-space (i.e. the regions of the $\alpha$–$m_{\textrm{max}}{}$ plane that have likelihood support), this factor remains fairly constant, implying that the inferred shape of the mass distribution is not affected by applying the $\langle VT\rangle{}$ calibration, although the overall rate is increased by about a factor of $\sim 1.6$ compared to the raw semi-analytic calculation. We have verified this explicitly by repeating the analysis with and without calibrated $\langle VT\rangle{}$.

For the redshift evolution analysis (Section IV), it is not sufficient to calibrate the mass-dependence of the detection probability; we must verify that the semi-analytic calculation reproduces the proper redshift-dependence. Therefore, we pursue an alternative modification to the raw semi-analytic calculation. In this modification, we replace the single PSD of the raw semi-analytic calculation with a different PSD calculated for the Livingston detector for each five-day chunk of observing time in O1 and O2. We find that this assumption correctly reproduces the redshift-dependent sensitivity empirically determined by the injection campaigns into the GstLAL pipeline for two fixed mass distributions (see Figure 12), whereas adopting different assumptions, such as using the PSDs calculated for the Hanford detector instead of the Livingston detector, or changing the single-detector SNR threshold away from 8, yields curves in Figure 12 that deviate noticeably from the distribution of recovered injections. This modification to the sensitivity calculation is necessary in the redshift analysis because the detection probability can fluctuate significantly at high redshifts $z>0.5$, where there is a very small probability of detection but considerable physical volume. Due to computational cost, the number of detections available at high redshift is insufficient to directly calibrate the redshift-dependent detection probability to injections as we did in the mass distribution section.

We find that between the two methods we use to estimate detection efficiency, the effect on the inferred mass distribution is negligible. However, the second time-varying approach employed in the redshift analysis underpredicts the overall merger rate by $\sim 70\%$ compared to the first calibrated approach (see the bottom right panel of Figure 11). This reflects a systematic uncertainty in the high-redshift detection efficiency and the implied merger rate. When additional detections lead to improved statistical constraints on the merger rate across redshift, it will become increasingly necessary to place a very large number of injections at high redshift and closely spaced in time in order to accurately estimate the high-redshift sensitivity.

We also note that all our calculations of the detection efficiency are based on the IMRPhenomPv2 waveform. Differences between the phasing, and more importantly, the amplitude of the waveform can lead to different SNRs and detection statistics for the same sets of physical parameters. To bound the significance of this effect, we carry out the injection-based $\langle VT\rangle{}$ estimation for both the IMRPhenomPv2 and SEOBNRv2 waveforms and find that for populations described by the two-parameter mass Model A, the waveforms produce $\langle VT\rangle{}$ estimates consistent to 10% across the relevant region of hyperparameter space with high posterior probability. Therefore, compared to the statistical uncertainties, the choice in waveform does not contribute a significant systematic uncertainty for the $\langle VT\rangle{}$ estimation.

Finally, an additional systematic uncertainty we have neglected in the $\langle VT\rangle{}$ and parametric rates calculations is the calibration uncertainty. While the event posterior samples have incorporated a marginalization over uncertainties on the calibration Farr et al. 2015 for both strain amplitude and phase, the $\langle VT\rangle{}$ estimation here does not. The amplitude calibration uncertainty results in an 18% volume uncertainty (Abbott et al. 2018), which is currently below the level of statistical uncertainty in our population-averaged merger rate estimate.

Appendix B Alternative Spin Models

We perform here a number of complementary analyses to reinforce the robustness of the results in Section V, and gauge the effect of fixed parameter choices on spin inferences. Instead of a parameterized model such as those used in Section V, we focus on a few discrete choices of model parameters to reinforce the conclusions in that Section. These choices provide a complementary view to the results presented earlier and also display our current ability (or inability) to measure features in differing parts of the mass and spin parameter space.

We choose a set of specific realizations of the general model described in Section II.2, building on Farr et al. 2017; Tiwari et al. 2018. Four discrete spin magnitude models are considered, the first three being special cases of Equation 4:

Low (L): $p(a)=2(1-a)$, i.e., $\alpha_{a}=1$, $\beta_{a}=2$.

Flat (F): $p(a)=1$, i.e., $\alpha_{a}=1$, $\beta_{a}=1$.

High (H): $p(a)=2a$, i.e., $\alpha_{a}=2$, $\beta_{a}=1$.

Such magnitude distributions are chosen as simple representations of low, moderate and highly spinning individual black holes. The very low (V) population is added to capture the features of an even lower spinning population — this is motivated by the features at low spin of the parametric distribution displayed in Figure 8.

For spin orientations we consider three fixed models representing extreme cases of Equation II.4:

Isotropic (I): $p(\cos t_{i})=1/2$; $-1<\cos t_{i}<1$, i.e., $\zeta=0$.

Aligned (A): $p(\cos t_{i})=\delta(\cos t_{i}-1)$, i.e., $\zeta=1$, $\sigma_{i}=0$.

Restricted (R): $p(\cos t_{i})=1$; $0<\cos t_{i}<1$, this is the same as I, except the spins are restricted to point above the orbital plane.

The isotropic distribution is motivated by dynamical or similarly disordered assembly scenarios, while the aligned one better capture a population of isolated binaries, under the simplifying assumption that the stars remain perfectly aligned throughout their evolution. In order to assess any preferences in the data for binaries with $\chi_{\rm eff}>0$, we introduce the restricted model: it resembles the isotropic distribution, but limits tilt angles to be positive. While we have mathematically defined the R model by assuming tilted spins, the same $\chi_{\rm eff}$ distribution can be generated with nonprecessing spins.

Here we perform our inference entirely through $\chi_{\textrm{eff}}{}$, whose 12 different distributions are illustrated in Figure 13. Since we do not have conclusive results on $\beta_{q}$ from Figure 3, we cannot make a single simplifying assumption on the mass model, which the $\chi_{\textrm{eff}}{}$ distribution depends on. We therefore consider three limiting cases: two of these fix the mass ratio to fiducial values, $q=1$ and $q=0.5$. The third corresponds to a fixed parameter model with $\alpha=1,m_{\textrm{min}}{}=5,m_{\textrm{max}}{}=50$. Figure 13 illustrates the $\chi_{\textrm{eff}}{}$ distributions implied by each of these scenarios.

Following Farr et al. 2017; Tiwari et al. 2018, we calculate the evidence and compute the Bayes factors for each of the zero dimensional spin models. Results are provided in Table 8, with the low and isotropic distribution (LI) as the reference.

Because of degeneracies in the GW waveform between mass ratio and $\chi_{\textrm{eff}}{}$, the choice of mass distribution impacts inferences about spins. This effect explains the significant difference in Bayes factors for the third row in the table. We find again our result moderately favors small black hole spins. The restricted models with $\chi_{\rm eff}$ strictly positive consistently produce the highest Bayes factors. For the small-spin magnitude models we cannot make strong statements about the distribution of spin orientations. Models containing highly spinning components are significantly disfavored, with high or flat aligned spins particularly selected against (e.g., FA and HA are disfavored with Bayes factors ranging in $\left[10^{-11},10^{-6}\right]$ and $\left[10^{-21},10^{-13}\right]$, respectively). As a bracket for our uncertainty on the mass and mass ratio distribution, we evaluated the Bayes factors for the fixed parameter model $\alpha=2.3$, $m_{\textrm{min}}{}=5$, $m_{\textrm{max}}{}=50$. They differ from the third mass model in Table 8 by a factor comparable to unity.

B.2 Spin Mixture Models

The models considered for model selection in Table 8 all assume a fixed set of spin magnitudes and tilts. There is no reason to believe, however, that the Universe produces from only one of these distributions. A natural extension is to allow for a mixing fraction describing the relative abundances of perfectly-aligned and isotropically distributed black holes spins.

We assume that the aligned and isotropic components follow the same spin magnitude distribution. It is possible that black holes with a different distribution of spin orientations would have a different distribution of spin magnitudes, but given our weaker constraints on spin magnitudes, we focus on spin tilts sharing the same magnitude distribution.

We compute the posterior on the fraction of aligned binaries $\zeta$ in the population as per Equation II.4 in the limit ($\sigma_{i}\rightarrow 0$). The models here are subsets of the Mixture distribution, with a purely isotropic being $\zeta=0$, and completely aligned being $\zeta=1$. The prior on the mixing fraction is flat.

All of the models which contain a completely aligned component favor isotropy over alignment. This ability to distinguish a mixing fraction diminishes with smaller spin magnitudes. This is because such spin magnitudes yield populations which are not distinguishable to within measurement uncertainty of $\chi_{\textrm{eff}}$. We do not include the most-favored restricted (R) configuration, but expect that the results would be similar. Coupled with the model selection results in the previous section, this implies that the mixing fraction is not well determined when fixed to the models (low and very low) which are favored by the data (see Figure 13). As stated above, in this case our ability to measure the mixing fraction is negligible.

B.3 Three-bin Analysis of χeff\chi_{\textrm{eff}}

We illustrate here how $\chi_{\textrm{eff}}$ measurements can provide insights into discerning spin orientation distributions. Following Farr et al. 2018, we split the range of $\chi_{\textrm{eff}}{}$ into three bins. One encompasses the fraction of uninformative binaries with $\chi_{\textrm{eff}}$ consistent with zero ($|\chi_{\textrm{eff}}{}|\leq 0.05$); the vertical axis of Figure 14 shows the fraction of binaries lying outside of this bin. The other two capture significantly positive ($\chi_{\textrm{eff}}{}>0.05$), and significantly negative ($\chi_{\textrm{eff}}{}<-0.05$) binaries. The width $0.05$ is chosen to be of the order of the uncertainty in a typical event posterior.

The aligned spin scenario is preferred in the posterior support on the right half of Figure 14: the small fraction of binaries which are informative tend to possess $\chi_{\textrm{eff}}$ greater than zero. Conversely, if the spins are isotropic, there would be no preference for positive or negative $\chi_{\textrm{eff}}{}$, and the posterior in Figure 14 would peak towards the middle. However, of the ten observed binaries, eight are consistent with zero $\chi_{\textrm{eff}}$ and only two are informative, thus demonstrating our ability to distinguish between the two scenarios is weak.

Appendix C Importance Resampling The Single-Event Likelihood

Our hierarchical population analysis uses the individual-event likelihood for each event $n=1,\ldots,N$, $\mathcal{L}\left(d_{n}\mid\xi,z\right)$ (see Section II, Eq. (II.6)). Individual-event analyses report posterior samples drawn a density that is proportional to this likelihood times a prior (Veitch et al. 2015; Abbott et al. 2018). The prior density used is uniform in detector frame masses and proportional to the square of the luminosity distance (Veitch et al. 2015); in terms of the source frame masses and redshift, the prior is

The derivative of the luminosity distance in a spatially flat universe (Hogg 1999) is

where $d_{H}=c/H_{0}$ is the Hubble distance and

Given a set of posterior samples as described above, we can transform them to samples from the likelihood over source frame masses and redshift by importance resampling with weights that are the inverse prior

The integral in Eq. (II.6) may then be approximated as