Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization

I. Introduction

Simulating the final stages of stellar binary evolution is computationally expensive. Additionally, there are significant theoretical uncertainties in key aspects of binary evolution, especially the common envelope phase (Ivanova et al. 2013) and supernova mechanism. For these reasons, populations of compact objects are simulated using population synthesis models (e.g., Dominik et al. 2015; Belczynski et al. 2017; Stevenson et al. 2017a). These are phenomenological models calibrated against a small number of more detailed stellar simulations.

There has been significant work using gravitational-wave data to infer the properties of black hole formation with ensembles of detections. These works range from comparing gravitational-wave data to specific, non-parameterized models (Mandel and O’Shaughnessy 2010; Stevenson et al. 2015; Dominik et al. 2015; Belczynski et al. 2016; Stevenson et al. 2017b; Zevin et al. 2017; Belczynski et al. 2017; Miyamoto et al. 2017; Farr et al. 2017a; Wysocki et al. 2017; Barrett et al. 2017), to attempts to group the data by binning, clustering or Gaussian mixture modeling (Mandel et al. 2017; Farr et al. 2017b; Wysocki 2017), to fitting physically motivated phenomenological population (hyper)parameters (Kovetz et al. 2017; Talbot and Thrane 2017; Fishbach and Holz 2017). In this work, we take the last approach and demonstrate that it is possible to identify physical features in the black hole mass spectrum with an ensemble of detections using phenomenological models, building on work in Kovetz et al. 2017 and Fishbach and Holz 2017.

Previous attempts to determine the binary black hole mass spectrum have employed one or more of these three approaches. Clustering is applied to a binned mass distribution in Mandel et al. 2017 to demonstrate that a mass gap between neutron star and black hole masses can be identified after $O(100)$ observations. In (Zevin et al. 2017; Stevenson et al. 2015; Barrett et al. 2017), the authors compare population synthesis models with different physical assumptions and show that predicted mass distribution can be distinguished using $O(10)$ of observations.

Kovetz et al. 2017 model the lower mass limit of black holes and use a Fisher analysis to demonstrate that it should be possible to identify the presence of the neutron-star black hole mass gap. They also model the distribution of mass ratios and propose a test for detecting primordial black holes. A method to simultaneously estimate the binary black hole mass spectrum and the merger rate is presented in Wysocki 2017. Wysocki 2017 also considers a Gaussian mixture model for fitting the distribution of compact binary parameters.

The rest of the paper is structured as follows. In section II, we introduce the statistical tools necessary to make statements about the black hole population. We then develop our model in section III in terms of population (hyper)parameters by considering current observational constraints and predictions from theoretical astrophysics and population synthesis. In section IV, we perform a Monte Carlo injection study. We consider how many detections will be necessary to identify different features using Bayesian parameter estimation and model selection. We show how the predicted mass distributions differ when using different (hyper)parameterizations. We also explore some of the consequences of using inadequate (hyper)parameterizations. In particular, we show that inadequate (hyper)parameterization can lead to significant bias in the estimate of the merger rate and the predicted amplitude of the stochastic gravitational-wave background (SGWB). Some closing thoughts are provided in section V.

II. Bayesian Inference

A binary black hole system is completely described by 15 parameters, $\Theta$. Recovering these parameters from the observed strain data requires the use of specialized Bayesian parameter inference software, e.g., LALInference (Veitch et al. 2015). The likelihood of a given set of binary parameters is computed by comparing the strain data to the signal predicted by general relativity. For the analysis presented here, the expected signal is calculated using phenomenological approximations to numerical relativity waveforms (Hannam et al. 2014; Schmidt et al. 2015; Smith et al. 2016). For a given set of strain data, $h_{i}$, LALInference returns a set of $n_{i}$ samples, $\{\Theta\}$, which are sampled from the posterior distribution, $p\left(\Theta\middle|h_{i},H\right)$, of the binary parameters, along with the Bayesian evidence, $\mathcal{Z}(h_{i}|H)$, where $H$ is the model being tested.

The distribution of binary black hole systems observed by current detectors is not representative of the astrophysical distribution of binary black holes. The observing volume of current gravitational-wave detectors is limited by the instruments’ sensitivity. The sensitive volume for a detector to a given binary is primarily determined by the masses of the black holes with spin entering as a higher order effect. More massive systems produce gravitational waves of greater amplitude. However, these more massive systems merge at a lower frequency and, hence, spend less time in the observing band of the detector. Additionally, distant sources undergo cosmological redshift and appear more massive than they actually are. Here, we will deal only with the un-redshifted “source-frame” masses, not the “lab-frame” masses observed by gravitational-wave detectors We note that the source/lab-frame distinction is about cosmological redshift and is not a statement about detectability and/or selection effects.

Accounting for these factors, we calculate $V_{obs}(\Theta)$, the sensitive volume for a binary with parameters $\Theta$, following Abbott et al. 2016c, using semi-analytic noise models corresponding to different sensitivities (Abbott et al. 2016e). The noise, and hence sensitivity, in real detectors is time-dependent and so calculating this volume requires averaging over the observing time to obtain a mean sensitive volume $\langle V_{obs}(\Theta)\rangle$ (Abbott et al. 2016c).

II.2. Population Inference

We are interested in inferring population (hyper)parameters describing the distribution of source-frame black hole masses. The formalism to do this is briefly described below (see e.g., Gelman et al. 2013, chapter 29 for a more detailed discussion of hierarchical Bayesian modeling and Mandel et al. 2014 for a discussion of selection biases).

Hierarchical inference of this kind can be cast as a post hoc method of changing from the prior distribution used in the single event parameter estimation to a new prior, which depends on population (hyper)parameters, $\Lambda$. We marginalise over all of the binary parameters while reweighting the posterior samples by the ratio between our (hyper)parameterized model and the prior used to generate the posterior distribution. This marginalisation integral is approximated by summing over the posterior samples for each event. The $N$ events are then combined by multiplying the new marginalised likelihood for the individual events,

We combine this likelihood with $\pi(\Lambda|H)$, the prior for the (hyper)parameters assuming a model $H$, and the Bayesian evidence for the data given $H$ to obtain the posterior distribution for our (hyper)parameters,

To perform the (hyper)parameter estimation we use the python implementation of MultiNest (Feroz et al. 2009; Buchner et al. 2014). Additionally, we calculate the posterior predictive distribution (PPD) of the binary parameters,

where $\Lambda_{k}$ are the $n_{k}$ (hyper)posterior samples. The PPD shows the probability that a subsequent detection will have parameters $\Theta$ given the previous data, $\{h\}$.

II.3. Model Selection

Model selection is performed in our Bayesian framework by considering Bayes factors,

A large Bayes factor, $BF^{\alpha}_{\beta}\gg 1$, indicates that $H_{\alpha}$ is strongly favored over $H_{\beta}$. We adopt a conventional threshold of $\ln BF=8$ to distinguish between two models.

III. Phenomenology

In this section, we develop a parameterization of the black hole mass spectrum using predictions from astrophysics theory, population synthesis models, and electromagnetic observations. In this way, we can relate gravitational-wave measurements to stellar astrophysics. For low-mass systems, the parameter that most strongly affects the observable gravitational waveform is a combination of the component masses known as the chirp mass, $\mathcal{M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$. For high-mass systems, the waveform is primarily determined by the total mass of the system. The mass ratio is more difficult to determine due to covariances between the mass ratio and the spin of the black holes.

The canonical assumed distribution of black hole masses is a power law distribution in the primary mass between some maximum and minimum masses. This power law distribution has three typical parameters: the spectral index $\alpha$, the minimum mass $m_{\min}$, and the maximum mass $m_{\max}$. The distribution of secondary mass is typically taken to be flat between $m_{\min}$ and $m_{1}$. We take this as the starting point for our parameterization.

Here, $m_{\min}$ is the upper limit of the NS-BH mass gap and $m_{\max}$ is the lower limit of the upper mass gap. The variable $m_{\text{PI}}$ is the mass above which stars undergo PISN leaving no remnant. Here, we assume $\alpha>1$ and $m_{\text{PI}}\gg m_{\max}$.

III.2. Low-Mass Binaries

The smaller sensitive volume for lower-mass binaries means that it is more difficult to probe the low-mass end of the black hole mass spectrum with gravitational-wave detections. Previous analyses of the black hole mass spectrum from gravitational-wave detections have assumed that the black hole mass spectrum has a sharp cut-off at some minimum mass $m_{\min}$. However, this overestimates the number of low-mass black holes if the distribution of low-mass black holes in merging binaries is the same as that in low-mass X-ray binaries (Özel et al. 2010). Population synthesis models also generically predict that the primary mass distribution peaks above the minimum mass.

We replace the step function at the low-mass end of the black hole mass spectrum with a smoothing function, $S(m,m_{\min},\delta m)$, which rises from zero at $m_{\min}$ to one at $m_{\min}+\delta m$,

Our model of the distribution of the primary mass can be summarized as

encodes the power-law distribution with a smooth turn on at low mass and

III.3. Mass Ratio

Previous analyses by the LIGO/Virgo scientific collaborations have assumed that the secondary mass is distributed uniformly between a lower limit set by $m_{\min}$ and an upper limit of $m_{1}$. This is motivated by observations of the stellar initial binary population, (e.g., Kroupa et al. 2013 and Belloni et al. 2017). In contrast to this, population synthesis models typically predict that the distribution of mass ratios should be biased towards equal mass binaries, (e.g., Belczynski et al. 2017). We model the distribution of the mass ratio as a power-law with spectral index $\beta$ as in Kovetz et al. 2017; Fishbach and Holz 2017. For a mass ratio distribution peaked at equal masses, $\beta>0$. We also impose the same smoothing at the lower limit as we apply to the primary mass. This allows us to write down the conditional probability distribution for secondary masses given a primary mass,

III.4. Summary

A table listing the (hyper)parameters and their physical meaning is provided in Tab. 1. Including a factor of $V_{obs}$ to account for selection biases, the probability of detecting a mass pair given our (hyper)parameters, $\Lambda=\{\alpha,m_{\min},m_{\max},\delta m,\lambda,m_{\text{pp}},\sigma_{\text{pp}},\beta\}$ and under model $H$, is

III.5. Other Effects

We do not expect either of these mechanisms to significantly affect the position and shape of an excess due to PPSN, although they complicate the interpretation of the maximum black hole mass. It is possible that black holes formed through repeated mergers could be identified on a case by case basis. For example, black holes formed by a binary black hole merger event are expected to have large dimensionless spins, $a\sim 0.7$ for equal mass non-spinning pre-merger black holes (Scheel et al. 2009).

IV. Monte Carlo Study

To enforce selection effects, we keep only binaries with optimal matched filter signal to noise ratio, $\rho>8$, in a single Advanced LIGO detector operating at design sensitivity (Abbott et al. 2016e). We generate a set of 200 events for our simulated universe. Each signal is then injected into a three detector LIGO-Virgo network with all detectors operating at their design sensitivities. Fig. 2 shows the distribution of primary masses and mass ratio in our simulated universe before (dashed) and after (solid) accounting for observation bias. The blue histogram indicates the injected values.

Using the recovered posterior distributions for the injected events, we employ the statistical methods described above for each of our models. The Bayes factors comparing $H_{5}$ to the others are enumerated in Table 3. In Table 3 we also give an approximate number of events needed to reach our threshold $\ln BF=8$, assuming linear growth of $\ln BF$ with number of detections. We consider two cases. “Cosmic” assumes zero measurement error. All uncertainty comes from cosmic variance. “Design” uses posterior samples obtained through running LALInference for a three detector network operating at design sensitivity. Including measurement errors reduces our resolving power between any pair of models by a significant factor for all the models. Unless otherwise specified we will refer to the Design Bayes factors. Below, we consider the effect of each of the modifications on the mass distribution model.

After 200 events, our model without the variable upper-mass cut-off, $H_{0}$, is disfavored with a log Bayes factor of $\sim 160$. We determine how many events are necessary to surpass the threshold of $\ln BF^{5}_{0}=8$ by considering subsets of our injection set. After 20 detections $\ln BF^{5}_{0}\sim N(\mu=13.3,\sigma=2.4$). Here, $N(\mu,\sigma)$ denotes a normal distribution with mean $\mu$ and variance $\sigma^{2}$.

Given this, we expect to be able to identify an upper-mass cut-off in the mass distribution after $\lesssim 20$ events. We note that $\mu$ grows linearly with number of events, this scaling is used in Table 3 to approximate the number of events to reach $\ln BF=8$. This is consistent with a similar study by Fishbach and Holz 2017.

IV.2. PPSN Peak

After 200 events, the posterior distribution on $\lambda$, the fraction of black holes formed through PPSN, is shown in Figure 4. We measure the maximum posterior probability point and 95% highest density confidence interval (HDI) to be $\lambda\sim 0.11^{+0.07}_{-0.04}$ (all future confidence regions will be 95% HDI unless specified) and disfavor $\lambda=0$ at $\gtrsim 3\sigma$. Correspondingly, $\ln BF^{5}_{2}=4.9$ which is moderate evidence for the existence of the PPSN peak, but below our threshold for a confident detection.

IV.3. Mass Ratio

The posterior distribution on $\beta$ is shown in Figure 4. After 200 events, the $1\sigma$ and $2\sigma$ confidence intervals on $\beta$ span $1.1$ and $2.2$ respectively. We disfavor $H_{3}$ with a log Bayes factor of $\ln BF^{5}_{3}=7.1$ after our 200 injections, just below our threshold of 8.

IV.4. Low Mass

IV.5. Mass Distribution Recovery

As a qualitative measure of the difference between the inferred mass distributions, we plot the posterior predictive distribution for the primary mass given the binaries in our injection study for our models in figure 6. The dashed black line indicates the injected distribution. We can see the effect of the different (hyper)parameterizations.

IV.6. Impact on the Merger Rate

The majority of binary black hole mergers are not individually resolvable by Advanced LIGO/Virgo. Using a (hyper)parameterization which does not accurately describe the true distribution leads to a biased estimate of the fraction of mergers which are individually resolvable and hence the merger rate (Abadie et al. 2010; Abbott et al. 2017a). Compact binary coalescences are a Poisson process which can be described by a merger rate $R(\Lambda)$. For a detector with time-independent sensitivity and a model of the distribution of binary black hole systems, the merger rate can be inferred from: the number of observed events $N$, the sensitive volume of our detectors $V(\Lambda)$, and the observation time $T$,

and $V_{obs}(\Theta)$ is the sensitive volume to a given binary introduced in section Sec. II.

IV.7. Impact on the Stochastic Background

The cross-correlation method is expected to require years of observation before the background can be resolved. Recently, a method involving searching directly for the stochastic background due to binary black hole mergers has been introduced in Smith and Thrane 2017. This method is expected to be able to detect this component of the background using days of data. Since this method relies on the rate of binary black hole mergers rather than $\Omega_{\text{GW}}$ it will be more sensitive to the black hole mass function than cross-correlation searches.

V. Discussion

We highlight several other interesting results that can be obtained using 200 detections at design sensitivity:

We will be able to identify the presence of an excess due to PPSN at $\sim 3\sigma$ and constrain the fraction of black holes forming through PPSN to within $\sim 0.05$ at 95% confidence.

We will be able to measure the power-law index on the mass ratio to within $\sim\pm 1$.

Detailed measurement of the low-mass end of the mass distribution will most likely require 1000s of detections and may have to wait for future detectors, e.g., the proposed Einstein Telescope (Punturo et al. 2010) or Cosmic Explorer (Abbott et al. 2017g).

We demonstrate that neglecting the presence of either a cut-off or a mass peak can lead to a mis-recovery of the astrophysical distribution of black holes in merging binaries. For example, the higher sensitivity of current detectors to high-mass binaries means that in order to fit the upper mass range well, the low-mass distribution is biased. This leads to incorrect estimates of the total binary black hole merger rate and the predicted amplitude of the SGWB. The amplitude of the SGWB is also sensitive to the distribution of mass ratios.

Our analysis assumes that a clear distinction can be made between binary black hole systems and other compact binaries. In reality, if there is not a well-defined mass gap between neutron stars and black holes, it will be non-trivial to distinguish between binary black hole, neutron star-black hole, and binary neutron star systems (Yang et al. 2017). Although, differences in, e.g., the spins of the component objects may enable this distinction (Littenberg et al. 2015). Our framework can be naturally expanded to include these other classes of compact binaries.

References