Reconstructing phenomenological distributions of compact binaries via gravitational wave observations

I Introduction

The Advanced Laser Interferometer Gravitational Wave Observatory (LIGO) Abbott et al. (2015) The LIGO Scientific Collaboration and Virgo Accadia and et al 2012; Acernese et al. 2015 detectors have and will continue to discover gravitational waves (GW) from coalescing binary black holes (BBHs) and neutron stars. Several tens of binary black holes and potentially neutron stars are expected to be seen in O3, LIGO’s next observing run, alone; and several hundreds more detections are expected over the next five years Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration; Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration. Already, the properties of the sources responsible – the inferred event rates, masses, and spins – have confronted other observations of black holes’ masses and spins Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration, challenged previous formation scenarios Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration; Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration, and inspired new models Mandel and de Mink 2016; Marchant et al. 2016; Rodriguez et al. 2016a; Bird et al. 2016 and insights Kushnir et al. 2016; Lamberts et al. 2016 into the evolution of massive stars and the observationally accessible gravitational waves they emit Dvorkin et al. 2016; Abbott et al. (2016b) The LIGO Scientific Collaboration and the Virgo Collaboration. Over the next several years, our understanding of the lives and deaths of massive stars over cosmic time will be transformed by the identification and interpretation of the population(s) responsible for coalescing binaries Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration; Barack et al. 2018; Wysocki et al. 2018a, because measurements will enable robust tests to distinguish between formation scenarios Mandel and O’Shaughnessy 2010 with present Rodriguez et al. 2016b and future instruments Breivik et al. 2016; Nishizawa et al. 2016.

During the first few years of discovery, substantial theoretical modeling challenges and the rapid pace of events suggest that GW observations could soon outpace theory. In this work, we introduce a flexible, concrete, and production-ready approach to infer compact binary merger rate and compact binary distribution, in the context of an (arbitrary) parametrized phenomenological model. We extend or employ previously proposed models Fishbach and Holz 2017; Talbot and Thrane 2018. We are motivated by how constraints on these phenomenological models enable us to address broad astrophysical questions—the mass and spin distribution of neutron stars and black holes, as imparted at their birth; the dominant formation mechanism for compact binaries, such as the role of dynamical versus isolated formation channels for binary black holes. To that end, we provide concrete demonstrations of how a few GW measurements will provide insights that enable sharp discrimination between proposed astrophysical alternatives, or measurements of their parameters. We use simple phenomenological arguments and calculations to characterize the information that these first few hundred observations should provide. Conversely, we provide simple approaches to extend our phenomenological approach in sophistication and complexity as several thousand compact binary mergers provide sharp constraints on their underlying properties. This approach complements inferences that work within a concrete model family as explored in other proof-of-concept investigations (see, e.g., Mandel and O’Shaughnessy 2010; O’Shaughnessy 2013; Stevenson et al. 2015; Belczynski et al. 2016a; Zevin et al. 2017; Barrett et al. 2018; Miyamoto et al. 2017; Wysocki et al. 2018a and references therein).

GW measurements probe only a selection-biased part of the compact binary distribution. Previously reported estimates of the overall compact binary event rate rely on extrapolation away from the observed population, using some fixed model for the compact binary mass distribution Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration. In fact, the compact binary mass distribution and inferred event rate are strongly coupled. This paper provides the first self-consistent approach to infer both the compact binary event rate and parameter distribution; then it describes and explains the expected correlation in an accessible way.

Several recent studies have explored how well GW measurements can constrain the mass and spin distribution of binary black holes O’Shaughnessy 2013; Wysocki 2017; Mandel et al. 2017; Kovetz et al. 2017; Talbot and Thrane 2017; Farr et al. 2017; Fishbach and Holz 2017; Gerosa and Berti 2017; Fishbach et al. 2017; Stevenson et al. 2017; Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration; Vitale et al. 2017a; Farr et al. 2018. Our approach is novel insofar as it reconstructs both the strongly correlated event rate and the parameter distribution, making our method a robust tool to assess astrophysical formation scenarios. In our modeling, we focus on measuring the black hole (BH) spin magnitude and misalignment distribution, as a method to probe the formation scenarios for binary BHs. As first described in Mandel and O’Shaughnessy 2010, GW provide a unique opportunity to distinguish between isolated and dynamic formation mechanisms: measurements of the spin properties of the BHs Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration; Rodriguez et al. 2016b; Vitale et al. 2017b; Stevenson et al. 2017; O’Shaughnessy et al. 2017; Talbot and Thrane 2017. The presence of a component of the BH spins in the plane of the orbit leads to precession of that plane. If suitably massive and significantly spinning, such binaries will strongly precess within the LIGO sensitive band. If BBHs are the end points of isolated binary star systems, they would be expected to contain BHs with spins preferentially aligned with the orbital angular momentum Kalogera 2000; O’Shaughnessy et al. 2017, and therefore rarely be strongly precessing. If, however, BBHs predominantly form as a result of gravitational interactions inside dense populations of stellar systems, the relative orientations of the BH spins with their orbits will be random, and some gravitational wave signals may be very strongly precessing. At this early stage, observations cannot firmly distinguish between these two scenarios, or more broadly other possible BBH formation mechanisms Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration. These include the evolution of isolated pairs of stars Belczynski et al. 2016a; Belczynski et al. 2010; O’Shaughnessy et al. 2012; Dominik et al. 2012; Mandel and de Mink 2016; Marchant et al. 2016, dynamic binary formation in dense clusters Rodriguez et al. 2016a, and pairs of primordial black holes BHs Bird et al. 2016; see, e.g., Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration and references therein. Loosely speaking, however, the isolated evolution and globular cluster formation scenarios are the most well-developed and verifiable using independent observational constraints. More broadly, precise measurements of their properties will provide unique clues into how BHs and massive stars evolve Vitale et al. 2017b; Stevenson et al. 2017; Rodriguez et al. 2016b; Farr et al. 2017; Wysocki et al. 2018b; Breivik et al. 2016; Nishizawa et al. 2016.

This paper is organized as follows. In Sec. II we describe our techniques to infer compact binary populations, building upon inferences about parameters of individual events. Unlike prior work, we simultaneously reconstruct the event rate, mass distribution, and spin (vector) distribution. In Sec. III, we demonstrate our our population inference strategy with two examples. In the first, we perform a full end-to-end analysis of a synthetic GW data generated from a synthetic population of astrophysically distributed sources. In the second, using a tool to mimic how well we could constrain parameters of a candidate GW signal, we perform a large-scale investigation into how well GW measurements could constrain the mass and spin distribution of binary black holes. We find that the mass and spin distribution can be tightly constrained with only a few tens of events. By virtue of explicitly exploiting only some of the available information, our estimates are necessarily conservative. In Sec. IV, we apply our method to the currently reported binary black hole population. For simplicity, assuming the reported events to date represent a fair sample of the results of LIGO’s first two observing runs (O1 and O2), we corroborate previous results, finding black hole spins are likely small and that the black hole mass spectrum may have an upper bound. Due to small BH spins, except for GW151226, we can extract no information about typical BBH spin-orbit misalignments. We emphasize our demonstration uses a nonfinal sample for LIGO’s O2 survey: depending on that survey’s results, applying our methods to final O2 results could produce substantially different astrophysical conclusions. In Sec. V we briefly discuss the accuracy to which population parameters can be determined, and the surprisingly significant role of waveform systematics in the near future. After summarizing our conclusions in Sec. VI, we supply three appendixes. In Appendix A, we describe a robust, extensible procedure for generating synthetic posterior distributions for proposed GW events. This open-source procedure could be widely used to assess the viability of GW measurements to distinguish between proposed astrophysical channels. A subsequent short Appendix B describes how to generate synthetic populations of selection-biased GW sources using this procedure. Next, in Appendix C, following on and extending previous work, we use toy models for both the measurement process and source population to illustrate how well GW observations will constrain the mass and spin distribution of compact binaries, likely providing robust insights into compact object formation (e.g., BH natal spins and maximum masses) and binary formation mechanisms (e.g., dynamical over isolated).

II Method

A coalescing compact binary in a quasicircular orbit can be completely characterized by its intrinsic parameters, namely its individual 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). In this work, we will also use the total mass $M=m_{1}+m_{2}$ and mass ratio $q$ defined in the following way:

We will also refer to two other commonly used mass parametrizations: the chirp mass $\mathcal{M}_{\text{c}}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ and the symmetric mass ratio $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$. With regard to spin, we define an effective spin Damour 2001; Racine 2008; Ajith et al. 2011, which is a combination of the spin components along the orbital angular momentum direction $\hat{L}$, in the following way:

where $\bm{S}_{1}$ and $\bm{S}_{2}$ are the spins on the individual BH. We will also characterize BH spins using the dimensionless spin variables

We will express these dimensionless spins in terms of Cartesian components $\chi_{i,x},\chi_{i,y},\chi_{i,z}$, expressed relative to a frame with $\hat{z}=\hat{L}$ and (for simplicity) at the orbital frequency corresponding to the earliest time of astrophysical interest (e.g., an orbital frequency of $\simeq 10\,{\rm Hz}$).

When necessary, compact binary parameters are inferred through the use of Bayesian analysis via Rapid parameter Inference on gravitational wave sources via Iterative Fitting (RIFT) Lange et al. 2018, which reproduces the results of standard Monte Carlo techniques described in Abbott et al. (2016c) The LIGO Scientific Collaboration and the Virgo Collaboration; Veitch et al. 2015 and references therein. For any event, fully characterized by parameters $x$, we can compute the (Gaussian) likelihood function $p(d|x)$ for detector network data $d$ containing a signal by using waveform models and an estimate of the (approximately Gaussian) detector noise on short timescales (see, e.g., Veitch et al. 2015; Abbott et al. (2016c) The LIGO Scientific Collaboration and the Virgo Collaboration; Abbott et al. (2016d) The LIGO Scientific Collaboration and the Virgo Collaboration and references therein). In this expression $x$ is shorthand for the set of 15 parameters needed to fully specify a quasicircular BBH. The posterior probability distribution is therefore $p(x|d)\propto p(d|x)p(x)$, where $p(x)$ is the prior probability of finding a merger with different masses, spins, and orientations somewhere in the universe. These parameters $x$ can and are often described with alternate coordinate systems. We sometimes refer to the source luminosity distance $d_{L}$ or equivalently its source redshift $z$, and to the detector-frame or redshifted masses $m_{i,z}=m_{i}(1+z)$. (To distinguish from the detector-frame masses, we will sometimes refer to $m_{i}$ as the source-frame binary masses.) LIGO-Virgo analyses have adopted a fiducial prior $p_{\rm ref}(x)$ that is uniform in orientation, in luminosity distance cubed, in redshifted mass, in spin direction (on the sphere), and, importantly for us, in spin magnitude Veitch et al. 2015; Abbott et al. (2016c) The LIGO Scientific Collaboration and the Virgo Collaboration. Using standard Bayesian tools Abbott et al. (2016c) The LIGO Scientific Collaboration and the Virgo Collaboration; Veitch et al. 2015, one can produce a sequence of independent, identically distributed samples $x_{n,s}$ ($s=1,2,\ldots,S$) from the posterior distribution $p(x|d)$ for each event $n$; that is, each $x_{n,s}$ is drawn from a distribution proportional to $p(d_{n}|x_{n})p_{\rm ref}(x_{n})$. Typical calculations of this type provide $\lesssim 10^{4}$ samples Abbott et al. (2016c) The LIGO Scientific Collaboration and the Virgo Collaboration; Veitch et al. 2015 from which the posterior probability distribution is inferred.

For other examples involving purely synthetic observing scenarios, we perform this procedure with a familiar Fisher matrix approximation for the form of $p(d|x)$ as a function of $x$ Cutler and Flanagan 1994; Poisson and Will 1995; Cho et al. 2013; see Appendix A for details.

II.2 Population inference

Ultimately we are interested in determining the likelihood of the astrophysical BBH population having a given merger rate $\mathcal{R}$ and obeying a given parametrization $\Lambda$, given the data for $N$ detections, $\mathcal{D}=(d_{1},\ldots,d_{N})$. This likelihood, $\mathcal{L}(\mathcal{R},\Lambda)\equiv p(\mathcal{D}\mid\mathcal{R},\Lambda)$, is that of an inhomogeneous Poisson process

Using Bayes’ theorem, $p(\mathcal{R},\Lambda\mid\mathcal{D})\propto p(\mathcal{R},\Lambda)\,\mathcal{L}(\mathcal{R},\Lambda)$, one may obtain a posterior distribution on $\mathcal{R}$ and $\Lambda$, after assuming some prior $p(\mathcal{R},\Lambda)$. To avoid computing the normalization constant, we instead draw samples from the posterior distribution via Goodman and Weare’s affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler Goodman and Weare 2010, as implemented in the Python package emcee Foreman-Mackey et al. 2013.

II.3 Estimate for V​TVT

Current LIGO-Virgo search sensitivity is well approximated by a familiar approximation: a source will typically be detected if the estimated signal to noise (SNR) of the second-most-sensitive detector is greater than $8$; see, e.g., Abadie et al (2010) The LIGO Scientific Collaboration and the Virgo collaboration and references therein. Using this approximation, one can directly evaluate the characteristic volume within which a source will be detected Finn and Chernoff 1993; for nonspinning BH binaries, this estimate is in reasonable agreement with detailed calculations of search sensitivity Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration. In this work, we therefore adopt the same approximation. Specifically, we estimate the orientation-averaged sensitive 3-volume $V$ to which a search is sensitive by the integral Abbott et al. (2016a) The LIGO Scientific Collaboration and the Virgo Collaboration; O’Shaughnessy et al. 2010a

where $p(\lambda\mid\Lambda)$ is the probability density function for a random binary in the Universe to have intrinsic parameters $\lambda$. In this expression, $\Lambda$ denotes the parameters that characterize the distribution from which all coalescing binaries are drawn. To calculate the horizon distance $D_{h}$ and hence $V$ for each combination of candidate binary parameters, we use the IMRPhenomD gravitational waveform approximation Husa et al. 2016; Khan et al. 2016.

The procedure described above allows us to estimate $V$ for any nonprecessing binary. Fig. 1 shows this estimate as a function of the component masses, based on a single LIGO detector operating at O1 sensitivity. Motivated by LIGO observations to date, however, we assume black holes will not be rapidly spinning. In these circumstances, spin has at best a modest impact on the sensitive volume; further complications due to precession would be expected to be smaller still Brown et al. 2012; O’Shaughnessy et al. 2010b.

Though we pursue a semianalytic estimate for $VT$ and hence the expected number of GW-detected events, detailed analysis of gravitational wave searches in real data with synthetic sources can evaluate $\mu$ and hence the search sensitivity directly Biswas et al. 2009; Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration; Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration; Tiwari 2018. Such an approach will be particularly necessary when search selection biases (e.g., due to detector noise non-Gaussianity) cause the search sensitivity threshold to deviate away from the simple SNR threshold described here.

II.4 Examples of phenomenological population models

Motivated by binary neutron star observations as well as the desire to reproduce arbitrary substructure and features in the mass distribution, we will also examine Gaussian mass distributions in component mass $m_{i}$

which is characterized by its mean value $\overline{m}$ and variance $\sigma_{m}$. In this work, we will typically explore the special case of $p(m_{1},m_{2})=p_{G}(m_{1})p_{G}(m_{2})$ and apply this distribution to the case of binary neutron stars, where the narrow width $\sigma$ relative to the mean $\overline{m}$ implies the distribution has effectively no support for undesirable regions (e.g., $m<0$). Finally, for complete generality, we also discuss mixtures of mass distributions, including Gaussian mixture models as previously employed in Wysocki 2017:

This latter approach allows complete generality and, with suitable smoothing priors on $w$, the ability to reproduce arbitrarily complicated mass distributions and circumvent systematic limitations due to our choice of model. In particular, these more generic models would allow us to reproduce features previously proposed in the literature, including overabundances at specific masses near the pair-instability supernova threshold Fraley 1968; Fryer et al. 2001; Woosley et al. 2002; Woosley et al. 2007; Kasen et al. 2011; Belczynski et al. 2016b.

For binary black hole spins, we adopt a simple flexible phenomenological model for each BH spin magnitude $\chi_{i}$: a beta distribution,

with unknown shape parameters $\alpha_{\chi_{i}}$ and $\beta_{\chi_{i}}$ ($i=1,2$). This tractable two-parameter distribution allows us to fit to the observed mean and variance—all that the sparse sample of existing observations will allow. In this work, we for simplicity assume both black hole spins are drawn from the same distribution and $\chi_{\rm max}=1$. Likewise, for simplicity we adopt the unphysical but easily described parametrization of the spin-orbit misalignment $\theta_{i}=\arccos\hat{\mathbf{L}}\cdot\hat{\mathbf{S}}_{i}$ proposed by Talbot and Thrane Talbot and Thrane 2017: a unimodal distribution based on a Gaussian in $\cos\theta$ that smoothly deforms into a uniform distribution in the limit of large $\sigma_{\chi_{i}}$:

When using this model, we assume the polar angles $\phi_{i}$ of each spin vector relative to the orbital angular momentum direction $\hat{L}$ are uniformly distributed between $0,2\pi$. In this work, we assume BH spins are drawn from the same spin misalignment distribution $\sigma_{\chi_{1}}=\sigma_{\chi_{2}}$. In this approach, as in our parameter inference, all spins are assumed specified at a gravitational wave frequency $f_{\rm ref}=20\,{\rm Hz}$. No compelling reason exists that astrophysical formation processes should cause binaries of different masses and spins to be drawn from a single, universal misalignment distribution at an arbitrary reference frequency $f_{\rm ref}$; see, e.g., Wysocki et al. 2018a; Rodriguez et al. 2018 for more detailed models. That said, this phenomenological approach is qualitatively consistent with the kinds of misalignments produced by binary SN natal kicks (e.g., $1-\cos\theta_{i}\lesssim 0.1$ for BH natal kicks of order $50\,{\rm km/s}$ O’Shaughnessy et al. 2017), allowing us a simple way to characterize whether observations support or disfavor plausible amounts of spin-orbit misalignment.

II.5 Useful phenomenological parameters

Observations will constrain combinations of these phenomenological parameters which reflect clear physical features in the observed (selection-biased) distribution of binary black holes. We can better characterize what we learn from GW observations early on by adopting coordinates conforming to these features.

In the context of our fiducial single-component model, we adopt a reference mass $m_{1}=m_{\rm ref}=15M_{\odot}$ and characterize the overall event rate not by its normalization, which depends on unobserved binaries with high and low masses, but by the event rate ${\cal R}p(m_{\rm ref})$ of binaries whose primary $m_{1}$ has a mass comparable to GW151226 Abbott et al. 2016a. We identify other natural coordinates for the distribution of $m_{1}$ via its detection-weighted cumulative distribution ${\cal P}(<m_{1})$:

For BH spins, closed-form expressions for the appropriate mean values and variances are generally not available for arbitrary selection biases $VT$; however, to the extent that $VT$ depends only weakly on BH spin, our model for BH spins and misalignments [Eqs. (10,11)] implies that

for our fiducial case where both BH spins are drawn from the same distributions; in these expressions, $\Sigma_{\chi}^{2}$ refers to the variance of the one-dimensional $\chi$ distribution, while $\bar{\chi}$ refers to its mean.

II.6 Interpreting results: Posterior predictive distributions and revised priors

If we ask any question about compact binary properties $x$ rather than model hyperparameters $\Lambda$, the only quantity that appears in our posterior inferences $p(\Lambda|\{d\})$ informed by our observations $\{d\}$ is the posterior predictive distribution $p_{\rm ppd}(x|\{d\})$:

The posterior predictive distribution (PPD) encodes our best estimates of the properties of any randomly selected future binary, based on observations to date and accounting for our initial prior knowledge about $\Lambda$. Unlike the model parameters themselves, which may be highly degenerate and lack physical meaning, the PPD provides an unambiguous estimate for how likely different binary parameters are, given our knowledge. Note that by design, the PPD is a probability distribution and, folding in all uncertainties, does not have an error estimate.

As events accumulate, we can use posterior constraints $p(\Lambda|\{d\}_{k})$ on model hyperparameters $\Lambda$ based on the first $k=1\cdots N$ observations to provide a nuanced, observationally revised perspective on future measurements $k>N$. These prior insights can be particularly powerful when individual future measurements are only weakly informative about certain binary parameters such as the mass ratio or spin; see, e.g., Vitale et al. 2017c; Williamson et al. 2017 for examples.

To be concrete, our usual population inferences are performed using a single fiducial choice of reference prior $p_{\rm ref}(x)=p(x|\Lambda_{\rm ref})$: the posterior is $p(x|d_{k},\Lambda_{*})=p(d_{k}|x)p(x|\Lambda_{*})/\int p(d|x)p(x|\Lambda_{\rm ref})$. We exploit prior measurements via

In this expression, the numerator $\int d\Lambda p(x|\Lambda)p(\Lambda|\{d_{k}\})$ is the posterior predictive distribution described above.

III Controlled tests with synthetic populations and measurements

To demonstrate our method can infer population parameters, we perform several validation studies using toy models which mimic key features of real gravitational wave observations. These completely controlled illustrations also let us highlight what can be inferred and why about the mass and spin distribution, within the context of our approach. Finally, these examples allow us to demonstrate how population inference can strongly inform the interpretation of individual future GW observations.

For each component of a binary neutron star (BNS), observations of galactic pulsars suggest that the component masses are drawn from a Gaussian distribution with mean $1.33M_{\odot}$ and standard deviation $0.09M_{\odot}$ Özel and Freire 2016. Observations of pulsars and theoretical models of pulsar spin-down suggest that if both NS are not recycled, then their dimensionless spins will be small [$\simeq\mathcal{O}(0.05)$]. Under the assumption that NS spins are parallel to their orbital angular momentum, we construct a synthetic population drawn from this phenomenological model; construct synthetic observations for each binary, recovering 13 synthetic sources based on a three-detector advanced LIGO/Virgo network using a threshold set by the second-most-sensitive detector’s recovered amplitude; perform full GW inference on each source using RIFT Lange et al. 2018; and, with the resulting posterior distributions, use the techniques of Sec. II to infer the underlying NS mass and spin distribution. In our reconstruction, we assume both components of a NS binary are independently drawn from a Gaussian distribution with unknown mean and variance; and with spins $\chi_{i,z}$ drawn from a beta distribution with unknown mean and variance, such that $|\chi_{i,z}|\leq 0.05$.

Fig. 2 shows the synthetic measurements used as inputs in our calculation. These synthetic measurements incorporate significant uncertainty in each source’s redshift, which contributes to the overall uncertainty in each binary’s chirp mass. For each neutron star in our synthetic population, we use the APR4 equation of state to calculate each neutron star’s tidal deformability $\lambda_{i}=\lambda(m|\text{APR4})$. We generate and recover our synthetic sources with IMRPhenomD_NRTidal Dietrich et al. 2019. Fig. 3 compares our recovered NS mass and spin distribution. When inferring source parameters, our waveform model and parameter inferences include the effects of NS tides, treating each NS tidal deformability $\lambda_{i}$ as a free parameter. Despite considerable uncertainties in each measurement, each BNS observation constrains that binary’s chirp mass reasonably well, to an accuracy $\sigma_{\mathcal{M}_{\text{c}}}\simeq 0.05M_{\odot}$, dominated by uncertainty in source redshift. Because GW measurements are only weakly informative about the mass ratio, these measurements each constrain the total mass to be $m_{1}+m_{2}\simeq 2^{6/5}\mathcal{M}_{\text{c}}$ to an accuracy $\sigma_{\mathcal{M}_{\text{c}}}2^{6/5}$; averaging all such observations, we can deduce the mean NS mass $\bar{m}$. With $n=13$ such measurements, we expect to constrain the mean mass of the population to a 1 standard deviation accuracy $\sqrt{\sigma_{\mathcal{M}_{\text{c}}}^{2}2^{12/5}/4+\sigma^{2}}/\sqrt{n}\simeq 0.027M_{\odot}$, which compares favorably to $0.02M_{\odot}$, the standard deviation of our Bayesian estimate for $\bar{m}$ . (A similar analysis shows that we constrain the NS population standard deviation $\sigma_{m}$ almost entirely through these one-dimensional chirp mass constraints.) Because GW measurements have a smaller statistical uncertainty than the astrophysical population width in total mass, the accuracy to which we constrain the mean NS mass is dominated by a simple frequentist error estimate ($\sigma/\sqrt{n}$), allowing us to reliably project the information we will extract about NS masses from future GW observations.

The measurement accuracy for GW measurements of BNS has been long known Poisson and Will 1995, and their implications for astrophysics (e.g., mass and BNS spin distributions) have been immediately apparent; see, e.g., O’Shaughnessy et al. 2014; Hannam et al. 2013; Zhu et al. 2018 and references therein. We provide the first end-to-end demonstration of how well binary NS population parameters can be measured, using a detailed waveform model at a level where waveform systematics should not dramatically impact the mass, spin, or tidal parameter inferences being performed. By contrast, many previous studies focusing on NS tidal deformation have demonstrated that waveform systematics could bias inferences Wade et al. 2014; Favata 2014; Lackey and Wade 2015, if not controlled. Only recently have systematic errors between waveform models diminished enough to enable consistent infererence; see, e.g., The LIGO Scientific Collaboration et al. 2019.

III.2 BBH mass and (precessing) spin distribution

To be consistent with the priors adopted in other work Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration, we express our results after reweighting to correspond to a Jeffries prior on the rate [$\pi(\mathcal{R})\propto\mathcal{R}^{-1/2}$]. Even with only 25 events drawn from a preferentially low-spin population, our calculations show that GW measurements should strongly constrain the mass and spin distribution of binary black holes

With 25 events, our population model has enough information to produce strong constraints on the underlying phenomenological distributions, even for parameters such as spin which are weakly constrained by individual measurements. Fig. 7 illustrates how informative these constraints can be about the spin distribution. This figure compares the true marginal distribution of $q,\chi_{\text{eff}}$ for the BH-BH population to our best (posterior predictive) estimate of that distribution. Even with only a few tens of detections, the estimate traces the general structure of the true distribution. In particular, we can clearly and unambiguously identify that a bias in the $\chi_{\text{eff}}$ distribution toward positive values suggests an underlying tendency toward alignment. Of course, our synthetic observations were intentionally drawn from the model family we use to fit it; in general, the underlying astrophysical distribution may have a form outside the model family we adopt, introducing small biases into our interpretation. Nonetheless, our analysis substantially generalizes previous proof-of-concept demonstrations on how well BH measurements can measure BH spin distributions, not being limited to a single spin magnitude, a discrete and restrictive family of orientation distributions, or similar strong prior adopted in previous investigations Stevenson et al. 2017; Vitale et al. 2017b.

Even with only 25 events, we strongly constrain the BH spin distribution, in both magnitude and orientation (Fig. 8). As described in Appendix C in greater quantitative detail, these two constraints are easily understood. For this synthetic analysis, the upper limit on spin follows from the $\chi_{\rm eff}$ distribution of recovered sources. Since our synthetic observations included no events with large $\chi_{\rm eff}$, we can be confident BH spins are not extremely large, since by chance we ought to have found one large value of $\chi_{\rm eff}$ out of $25$, even allowing for uncertainty in how they are oriented. Similarly, because our synthetic population is preferentially aligned ($\sigma=0.4$), the recovered population shown in Fig. 2 has a $\chi_{\rm eff}$ distribution biased toward positive values. Using Eq. (14) for $\overline{\chi}_{\rm eff}$, the bias in $\chi_{\rm eff}$ inevitably implies $\cos\theta$ is preferentially positive and, as described in Appendix C, allows us to limit $\sigma$.

In this analysis, we employ conservative synthetic posteriors which assume only the chirp mass, mass ratio, and effective spin can be constrained with GW measurements. Precessing, coalescing binaries can produce a rich symphony of gravitational waves just prior to and during merger, reflecting complex binary dynamics and strong-field multimodal radiation. Given the high expected event rate in ongoing gravitational wave surveys, we expect that future observations will provide clear examples of precessional dynamics, if nature produces them, and that these measurements will allow us to much more sharply constrain the BH spin distribution. However, for massive BH binaries, model systematics complicate attempts to measure BH parameters, including spin. We will conduct full end-to-end calculations with synthetic data and state of the art models in future work.

IV Analysis of reported observational results

Fig. 9 shows our best estimates for the merger rate of BH-BH binaries of different masses, inferred within the context of the model described in Table 1 and demonstrated on synthetic data in Sec. III.2. Naturally, we estimate an overall BH-BH merger rate and mass distribution consistent with previously reported results Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration. Using a Jeffries’ prior for the merger rate, we find $\mathcal{R}=122^{+291}_{-96}\,{\rm Gpc}^{-3}{\rm yr}^{-1}$ based on O1. For O2, we find uncertainty in the event rate is reduced by roughly a factor of 2, both through reduced Poisson error (e.g., six instead of three events) and through sharper constraints on the mass distribution (e.g., reducing prospects for a large maximum mass). Our result for O1 is more conservative (wider) than the power-law result reported previously in Abbott et al. Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration, $97^{+135}_{-67}\,{\rm Gpc}^{-3}{\rm yr}^{-1}$, because we employ a more flexible model and therefore incorporate more model systematics, notably including the correlation between event rate and mass spectrum and also the impact of the upper mass cutoff. Conversely, if we employ consistent assumptions, we arrive at the same answers previously reported for O1 Abbott et al. (2016) The LIGO Scientific Collaboration and the Virgo Collaboration. As we adopt a merger rate model that reduces to previously investigated power laws, by design we reproduce the analysis reported in Fishbach and Holz 2017: the events reported during O2 suggest the absence of very massive BHs in the observable population. While our assumptions about the mass distribution model have modestly changed relative to Fishbach et al. Fishbach and Holz 2017, we reproduce their results when adopting the same inputs and mass model. For this reason our inferences about the mass spectrum exponent $\alpha_{m}$ are considerably wider than prior work which does not take a possible upper mass cutoff into account. Even with the small sample publicly reported so far, our analysis corroborates the analysis in Fishbach and Holz 2017 that O2-scale GW measurements could be weakly informative about the maximum mass of coalescing BHs.

As demonstrated in several previous investigations Farr et al. 2017; Wysocki et al. 2018a, we know that BHs in merging binaries likely have low typical spin. For example, based on the distribution of $\chi_{\rm eff}$, Farr et al. Farr et al. 2017 argued that several members of a discrete array of candidate spin orientations (aligned or isotropic) and magnitude distributions are inconsistent with observations to date, and that BH spins were likely randomly oriented or small. Later, Wysocki and collaborators Wysocki et al. 2018a demonstrated that, if binary black holes arose from isolated binaries whose spins were weakly misaligned by SN natal kicks, then only relatively small BH natal spins were consistent with observations available at the time. As shown in Fig. 10, with more events available to our analysis, and using much more flexible models, we can draw sharper and more generic conclusions about the BH spin distribution, even using only six reported events. First and foremost, exactly as seen with synthetic data, the absence of large $\chi_{\rm eff}$ allows us to with increasing confidence bound above the fraction of BHs in merging binaries that have large spin. Too, because collectively the observed population distribution of $\chi_{\rm eff}$ remains nearly symmetrically distributed around zero, we can with increasing confidence bound the fraction of binaries that are preferentially aligned and with modest spin. With at least one BH known to have spin (GW151226) and for simplicitly assuming the BH spin and mass distribution are uncorrelated, we are led to weakly disfavor scenarios where BHs are preferentially aligned (i.e., small $\sigma$ is disfavored). We emphasize, however, that this conclusion is driven by the absence of strong support for any spin in all but one binary (GW151226). We would arrive at the same nominal conclusion for a comparable number of random draws from a binary population model with perfectly aligned binaries with small BH spins. Future and more informative observations of BH binaries could significantly alter this conclusion.

V Discussion

In this work, we present concrete examples for how well just a handful of GW measurements can improve our phenomenology of the BH mass and spin distribution. Our examples include real observational data from LIGO’s O1 and (an incomplete sample from) O2 observing run, suggesting current observations could be on the cusp of constraining BH spins and maximum masses. We provide simple estimates to understand how well these parameters have been constrained, allowing the reader to extrapolate to larger sample sizes. For example, in the absence of positive support for spin, the upper limit on BH spin will decrease rapidly, allowing us to place strong upper limits for (or enable discovery of) BH natal spin.

Because each empirical marginal distribution possesses an infinite number of degrees of freedom, any phenomenological parametrization such as our own can quickly be exhausted by the data O’Shaughnessy 2013, particularly when the population must reproduce multiple observational features. In the short run, therefore, we anticipate a fully generic and regularized infinite-dimensional approach will soon be required to adequately reproduce the thousands of events that even the current generation of instruments will discover. A fully generic approach, however, can easily be misled, not least because GW measurements are subject to many subtle strong-field systematics due to model incompleteness. For example, a waveform approximation widely used for rapid parameter inference of binary black holes (IMRPv2 Hannam et al. 2014) omits astrophysically critical degrees of freedom—the calculation allows for only one precessing spin instead of the two necessary to fully describe the dynamics—and demonstrably has systematic errors large enough to shift posterior distributions for O3-scale events by an appreciable fraction of their statistically expected extent Williamson et al. 2017; Lange et al. 2018. To illustrate the pernicious impact of these systematic biases, we can consider a simple order-of-magnitude estimate: a single quantity, with intrinsic Gaussian distribution of mean $\mu$ and width $\sigma$, being observed multiple times by an apparatus with a (Gaussian, random) measurement error $\Delta x$ and bias $\delta x$. The bias will be important when it influences our best estimate of the average (i.e., when $\delta x\gtrsim\sqrt{\sigma^{2}+\Delta x^{2}}/\sqrt{N}$). Applying this order-of-magnitude approach to GW measurements, we expect that after only a few tens of binary mergers, these modeling systematics will progressively contaminate the interpretation of coalescing binaries, as posterior biases in each event become reflected in biases in the inferred population distribution. Waveform systematics will be even more important because BH spins appear to be small: greater accuracy is needed to separate the secular effects of spin. In this work, when carrying out a full parameter inference, we use the newly developed RIFT parameter inference engine Lange et al. 2018 to produce posteriors. We will discuss the impact of waveform systematics on BH spin misalignment measurements in future work.

VI Conclusions

We have introduced a flexible, ready-to-use, and self-consistent parametric method to estimate the compact binary merger rate as a function of binary parameters, specifically emphasizing mass and spin. Unlike prior work, our procedure self-consistently estimates the merger rate and binary parameter distribution, accounting for statistical sampling error, measurement error, and selection bias. Using this procedure, we show by example that only a handful of NS-NS and BH-BH measurements can enable strong constraints on their respective populations via GW observations alone. Even in the astrophysically likely scenario of small BH spin, we emphasize that just a few measurements will enable sharp constraints on the BH spin distribution. Interpreting current observations, we show that GW measurements are already beginning to place astrophysically interesting constraints on the spin of BHs. We reproduce prior results about the lack of reported BHs at high mass and its implications for the BH mass spectrum. Finally, particularly in our appendix, we explain how to extrapolate toward the measurement prospects available in the very near future.

The procedure described here assumes all sources have been unambiguously resolved from observational data, omitting any treatment of source significance aside from a naive selection bias. Farr et al. Farr et al. 2015 demonstrated and popularized an approach to self-consistently perform the detection and population inference process, estimating the foreground and background distributions simultaneously; see also Loredo 2004; Buchner et al. 2015; Messenger and Veitch 2013. Recently, Gaebel and collaborators Gaebel et al. 2019 developed a concrete procedure to apply this technique to gravitational wave observations. Owing to many deep similarities between our strategies, we anticipate we will shortly incorporate this technique in our own analysis.

The approach described here also employs several strong assumptions about the (lack of) correlations between model parameters. For example, our fiducial BH model assumes the mass-dependent BH merger rate is independent of redshift; that BH masses and spins are completely independent; and that BH spin misalignment and spin magnitudes are likewise uncorrelated. We will explore more physically motivated correlations in future work.

In the long run, phenomenology is only as sound as the underlying parametrization. Previous analyses have repeatedly shown that adopting an overly restrictive model will produce biased results, as demonstrated by Fishbach et al (with the maximum mass) Fishbach and Holz 2017 and Talbot et al Talbot and Thrane 2018 (with the shape of the maximum mass cutoff). With sufficient data, a suitably regularized infinite-dimensional parametrization will make unintended systematic biases less frequent. Mature methods for infinite-dimensional or nonparametric inference exist Gelman et al. 2013; Orbanz and Teh 2010; Ghosal and van der Vaart 2017, beginning with simple infinite-dimensional parametrizations plus smoothing priors or with Gaussian processes Rasmussen and Williams 2006. Early investigations have applied nonparametric methods to GW population estimates Wysocki 2017; Mandel et al. 2017. However, because the GW signal is so rich, many parameters can be measured for each event, several of which are believed to be correlated in most astrophysical formation scenarios. These correlations should be more sharply identified with strong theoretical priors for the immediate future.

Appendix A Mock posterior populations precessing binaries: Aligned Fisher matrix approach

We test our code using synthetic or “mock” posterior distributions for binary black hole parameters, designed to mimic the results of full end-to-end Bayesian inference on synthetic data. For the mock BBH posterior distributions constructed in this work, we adopt a very simple approximation, motivated by decades of experience suggesting that for short BBH signals the likelihood for gravitational wave signals is nearly Gaussian in three coordinates ($\mathcal{M}_{\text{c}},\eta,\chi_{\rm eff}$) and does not strongly constrain any other degrees of freedom. Specifically, if $\lambda_{0}$ are the true binary parameters and $\rho$ is the true network signal amplitude; if $\Gamma_{ab}=\left\langle\partial_{a}h|\partial_{b}h\right\rangle$ is the Fisher matrix for the binary parameters $\lambda$, evaluated at $\lambda=\lambda_{0}$ and for a signal amplitude $\rho$ using a fiducial detector power spetcrum; and if $p(\lambda)$ is the prior distribution on $\lambda$, then we approximate the posterior distribution by a distribution proportional to

where $\lambda_{*}$ is a fixed random realization from a normal distribution with mean $\lambda_{0}$ and covariance matrix $\Gamma^{-1}$. We generate samples from this distribution via Monte Carlo techniques. We evaluate the approximate Fisher matrix $\Gamma$ using the effective Fisher technique O’Shaughnessy et al. 2014; Cho et al. 2013; Cho and Lee 2014, applied to a nonprecessing binary waveform model assigned the same values of $\mathcal{M}_{\text{c}},\eta,\chi_{\rm eff}$ (i.e., via $\chi_{1,z}=\chi_{2,z}=\chi_{\rm eff}$).

This approximate posterior distribution has several distinct advantages. First and foremost, it captures in $\Gamma_{ab}$ the strong, parameter-dependent, and well-understood correlations between the variables that most significantly impact the GW inspiral signal, while simultaneously populating all intrinsic binary parameters. For example, it captures the shape of the posterior distribution in mass ratio and spin while correctly accounting for parameter boundary effects, as described in Ng et al. 2018. Second, it accounts via $\lambda_{*}$ for the effect of random noise realizations, which impact the best-fitting parameters associated with each set of synthetic data. By including an explicit prior $p_{\rm ref}(\lambda)$, it allows us to carefully adopt fiducial prior assumptions, which have a substantial impact on inferred binary masses and spins.

A ready-to-use implementation of this algorithm is available. See https://git.ligo.org/daniel.wysocki/synthetic-PE-posteriors.

For simplicity, in this implementation, no cosmological effects are applied. If used unaltered, this approximate posterior applies either if cosmological redshift effects are small compared to the width of the distribution in mass (i.e., bias is small compared to the statistical uncertainty) or if these ambiguity distributions are used to approximate the source-frame ambiguity function. Cosmological effects dominate the accuracy to which a binary neutron star’s chirp mass can be measured; to be used in such a scenario, this approximation must be refined to reflect the significant impact of the sources’ unknown redshift.

Appendix B Mock populations

Appendix C Overview of key phenomenological constraints

Classical frequentist statistical methods provide a quick way to assess how rapidly observations will constrain model hyperparameters. For example, the sample mean of maximum likelihood estimators converges rapidly to the true mean, and (to a first approximation) the sample variance is approximately $\chi^{2}$ distributed. Thus, by adopting the mean and variance of our underlying distributions as coordinates on the space $\Lambda$ of hyperparameters, we can estimate how efficiently observations will constrain them. For example, if we account for measurement error, we can measure the mean spin to an accuracy $\sqrt{V(\chi)+\sigma_{\chi}^{2}}/\sqrt{N}$ where $V(\chi)$ is the variance of the spin magnitude distribution and $\sigma_{\chi}$ is the typical spin measurement accuracy for the mass range of interest [typically $\mathcal{O}(0.3)$]. Because of sharp cutoffs, the maximum and minimum masses have a qualitatively different behavior; see, e.g., Amari and Nagaoka 2007. Both the maximum and minimum masses are best estimated using the most extreme individual event, with an accuracy converging as $1/N$. In our context—the power-law mass distribution—the accuracy with which these maximum masses can be determined scales directly with the number of events in a given region. We therefore expect the maximum mass can be determined to an accuracy of order $m_{\rm max}/N$; the appropriate scale factor can be calibrated to detailed analyses of the kind performed in Sec. III. Similarly, as described below in Appendix C.2, we can use the observed range of $\chi_{\rm eff}$ to constrain spin magnitudes and misalignments.

While providing a useful order-of-magnitude estimate into how well we can measure distribution parameters, the simple estimates above become cumbersome when trying to capture correlations between our phenomenological parameters, notably the event rate and mass distribution. Following O’Shaughnessy 2013, we assess how well we can distinguish model hyperparameters from the (expected) log-likelihood as a function of model hyperparameters $\Lambda$ of

where the expectation is performed relative to some reference model characterized by parameters $\Lambda_{*}$ such that $p_{*}(\lambda)\equiv p(\lambda|\Lambda_{*})$ and $\mu_{*}=\mu(\Lambda_{*})$. Rather than work in full generality, we perform a Taylor series expansion of the likelihood around the local maximum, characterizing the second order term by its inverse covariance or Fisher matrix $\Gamma_{ab}$

If $\gamma_{k}$ are eigenvalues of $\Gamma$, then hyperparameters can be measured to an accuracy $1/\sqrt{\gamma_{k}}$, which scales as $1/\sqrt{N}$ for $N$ the number of observed events.

For the power-law model described above, the only two derivatives needed are $\partial_{\ln R}X=1$ and $\partial_{\alpha}X=\partial_{\alpha}\ln{\left\langle VT\right\rangle}_{\alpha}$, the latter of which can be well approximated by $-1$. This term introduces correlations between the rate variable ($\ln{\cal R}$) and shape ($\alpha$). Conversely, using coordinates $\mu$ and $\alpha$ to characterize the observed population, by construction our inferred posterior distribution on the total number and mass distribution are uncorrelated.

Roughly speaking, the effects of measurement error add in quadrature in the Fisher matrix:

We can therefore refine the estimates provided above to incorporate simple estimates of GW measurement errors and their correlations. For the simple power-law estimate described above, however, these measurement errors are relatively small compared to the range of the distribution, unless $\alpha$ is very large.

In the above order-of-magnitude discussion, we have not accounted for parameter-dependent selection bias. To a good first approximation, GW selection bias enters only through the masses, roughly as the (chirp) mass to a power. We can therefore treat the observed population as a (different) power law, which observations constrain to an accuracy loosely characterized by the analysis above.

Therefore, for the power-law mass distribution, we expect the posterior distribution of (log) rate and powerlaw exponent will be correlated and follow a Gaussian distribution characterized by the inverse covariance

relative to the coordinates $(\ln{\cal R},\alpha)$, if we adopt a uniform prior on $\alpha$ and $\ln R$. This expression captures the correlations between the rate and mass ratio seen in our inferences, when only varying the total event rate and mass ratio.

C.2 Semianalytic model for constraints on the spin magnitude and misalignment distribution

In this paper, for the purposes of illustration and as a leading-order approximation suitable for the BH-BH binaries reported to date, we adopt three simplifying approximations: that the sensitive volume depends weakly on spin; that GW measurements will only constrain $\chi_{\rm eff}$; and that the underlying mass and spin distributions of BH-BH binaries are uncorrelated. In this framework of approximations, only $\chi_{\rm eff}$ measurements and hence the underlying $\chi_{\rm eff}$ distribution of the population determines how well we can distinguish between population models via spin measurements. Within this framework, we can simply and largely analytically estimate how much information we gain about the BH spin distribution from repeated measurements.

In our synthetic model (and nature) where BH spins appear to be small, the first few measurements will principally inform our upper limit on the BH spin distribution, via the absence of observations consistent with large $\chi_{\rm eff}$. For example, in our synthetic model, the 90% upper limit expected in 25 events is $\chi_{\rm eff}<0.31$; for our inferred posterior predictive distribution based on all published events, it is $0.19$. In Fig. 11, we use a simple toy model to illustrate how upper limits loosely inform our estimates of the BH spin distribution. In this model, we assume each BH in a binary has a random spin magnitude drawn from a uniform distribution between 0 and $\chi_{\rm max}$, randomly (isotropically) oriented, for binaries with a random mass ratio uniformly drawn between $0.1$ and $1$. This figure shows the cumulative distribution of $\chi_{\rm eff}$ implied by these assumptions, for different choices of $\chi_{\rm max}$. These cumulative distributions are well approximated by analytic expressions for the cumulative distribution of $\chi_{1,z}$ and $\chi_{\rm eff}$ under these assumptions; see Lange et al. 2018 for concrete expressions. For comparison, the vertical shaded regions show the largest values of $\chi_{\rm eff}$ which have significant support in our synthetic sample ($\chi_{\rm eff}\lesssim 0.5$), consistent with the largest plausible spins reported for O1 and O2 events. The lack of support for large $\chi_{\rm eff}$ in any observation to date strongly suggests that BH spins cannot be large. Conversely, an observation of a binary with $\chi_{\rm eff}$ bounded below by $\epsilon$ (e.g., GW151226) implies that a significant fraction of BH spins must be greater than of order $\epsilon$.

We emphasize that we provide these estimates (and perform our calculation within these underlying approximations) to produce a conservative, well-understood benchmark for how well the BH spin distribution can be constrained with present and future GW measurements. Real GW measurements, particularly of low-mass or closer and therefore higher-amplitude BH-BH mergers, will provide additional direct constraints on the other spin degrees of freedom.

Appendix D End-to-end tests of population hyperparameter recovery: PP–PP plots

In addition to our population inference code, we made $P$–$P$ plots for our synthetic parameter estimation code, described in Appendix B, as our population inference tests make use of it. Here we generated $k=1\cdots 1000$ synthetic BBH signals, drawing true values from the prior we used for measuring the posteriors. We repeated the same process just described, making posterior distributions on the intrinsic parameters $\lambda$, and evaluating the marginal cumulative distribution functions at the true values $\lambda^{*}_{k}$. $P$–$P$ plots for some representations of the intrinsic parameters are shown in the bottom panel of Fig. 12.

References