Measuring the Binary Black Hole Mass Spectrum with an Astrophysically Motivated Parameterization

A binary black hole system is completely described by a set of 15 parameters, Θ. 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, {Θ}, which are sampled from the posterior distribution, $p({\rm{\Theta }}| {h}_{i},H)$, 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 observing volume of a gravitational-wave detector 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 unredshifted "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(Θ), the sensitive volume for a binary with parameters Θ, following Abbott et al. (2016c), using semianalytic 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}_{\mathrm{obs}}({\rm{\Theta }})\rangle$ (Abbott et al. 2016c).

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, Λ. We marginalize 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 marginalization integral is approximated by summing over the posterior samples for each event. The N events are then combined by multiplying together the new marginalized likelihoods for the individual events,

$$\begin{eqnarray}&&{ \mathcal L }(\{h\}{}_{i=1}^{N}| {\rm{\Lambda }},H)\propto \displaystyle \prod _{i}^{N}\displaystyle \sum _{j}^{{n}_{i}}\displaystyle \frac{\pi ({{\rm{\Theta }}}_{j}^{i}| {\rm{\Lambda }},H)}{\pi ({{\rm{\Theta }}}_{j}^{i}| \mathrm{LAL})}.\end{eqnarray} \tag{ 1 }$$

Here $\pi ({\rm{\Theta }}| \mathrm{LAL})$ is the prior probability distribution used for single-event parameter estimation. The distribution $\pi ({{\rm{\Theta }}}_{j}^{i}| {\rm{\Lambda }},H)$ is the probability of a binary having parameters Θ in our model; see Section 3. We do not model any black hole parameters other than source-frame mass, and so Θ can be replaced by m = (m₁, m₂) in Equation (1) and all following equations. The prior distribution of source-frame masses in LALInference is a convolution of a prior that is uniform in component mass between limits determined by the reduced order model (Smith et al. 2016) and the redshift distribution that is uniform in luminosity distance extending out to 4 Gpc. This distance distribution is converted to redshift assuming ΛCDM cosmology using the results from the Planck 2015 data release (Planck Collaboration et al. 2016).

We combine this likelihood with $\pi ({\rm{\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,

$$\begin{eqnarray}&&p({\rm{\Lambda }}| \{h\}{}_{i=1}^{N},H)=\displaystyle \frac{{ \mathcal L }(\{h\}{}_{i=1}^{N}| {\rm{\Lambda }},H)\pi ({\rm{\Lambda }}| H)}{{ \mathcal Z }(\{h\}{}_{i=1}^{N}| H)},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{ \mathcal Z }(\{h\}{}_{i=1}^{N}| H)=\int d{\rm{\Lambda }}\,{ \mathcal L }(\{h\}{}_{i=1}^{N}| {\rm{\Lambda }},H)\pi ({\rm{\Lambda }}| H).\end{eqnarray} \tag{ 3 }$$

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,

$$\begin{eqnarray}p(m| \{h\}{}_{i=1}^{N},H) & = & \displaystyle \int \,d{\rm{\Lambda }}\pi (m| {\rm{\Lambda }},H)p({\rm{\Lambda }}| \{h\}{}_{i=1}^{N},H)\\ & \approx & \displaystyle \frac{1}{{n}_{k}}\displaystyle \sum _{k}^{{n}_{k}}\pi (m| {{\rm{\Lambda }}}_{k},H),\end{eqnarray} \tag{ 4 }$$

where Λ_k are the n_k (hyper)posterior samples. The PPD shows the probability that a subsequent detection will have parameters Θ given the previous data, {h}.

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

$$\begin{eqnarray}&&{\mathrm{BF}}_{\beta }^{\alpha }=\displaystyle \frac{{ \mathcal Z }(\{h\}{}_{i=1}^{N}| {{\rm{H}}}_{\alpha })}{{ \mathcal Z }(\{h\}{}_{i=1}^{N}| {{\rm{H}}}_{\beta })}.\end{eqnarray} \tag{ 5 }$$

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

The observation of binary black hole mergers through the detection of gravitational waves revealed the presence of a previously unobserved population of black holes with mass ∼30 M_⊙ (Abbott et al. 2016a, 2017a, 2017b). Since gravitational-wave detectors can observe more massive binaries at greater distances, binaries containing larger black holes are preferentially detected over less massive systems. Fishbach and Holz (2017) note that, given the observed rate of binaries with mass ∼30 M_⊙, it is somewhat surprising that we have not seen more massive black holes. They propose that this is due to a cutoff in the black hole mass spectrum around this mass. By comparing the Bayesian evidence for m_max = 41 M_⊙ and m_max = 100 M_⊙ using the first four events, they find tentative support for a cutoff. They further show that it will be possible to identify the presence of an "upper mass gap" with a Bayes factor of ≳150 (ln BF ≈ 5) using 10 detections, and the cutoff mass can be measured with ∼40 detections.

The theoretical motivation for such a cutoff is pulsational pair-instability supernovae (PPSNe) (Heger & Woosley 2002; Woosley & Heger 2015), whereby large amounts of matter are ejected prior to collapse to form a black hole. The expected result of this process is that all stars with initial mass 100 M_⊙ ≲ M ≲ 150 M_⊙ form black holes with masses ∼40 M_⊙. Stars with 150 M_⊙ ≲ M ≲ 250 M_⊙ are expected to undergo pair-instability supernovae (PISNe) and leave no remnant. Hence, we expect a gap in the black hole mass spectrum between ∼40 M_⊙ and ∼250 M_⊙ along with an excess of black holes at some mass m_pp ∼ 40 M_⊙. The excess is due to the stars with 100 M_⊙ ≲ M ≲ 150 M_⊙, which become PPSNe. The size, position, and shape of the peak are determined by the unknown details of PPSNe. While the observation of cutoff near 40 M_⊙ could be interpreted as evidence for PPSNe, the additional observation of a peak would provide a smoking-gun signature that the highest mass stellar binaries are reduced in mass via PPSNe.

Thus, we extend our description of the upper end of the black hole mass spectrum to allow for the possibility of an excess due to PPSNe. We model this as a normal distribution with unknown mean m_pp ∈ [25, 100]  M_⊙ and variance σ_pp < 5 M_⊙. It is also necessary to introduce a mixing fraction parameter, λ, which describes the proportion of binaries that are drawn from the normal distribution.

We expect that any mass peak due to PPSNe should be near the high-mass cutoff, which corresponds to m_max ≈ m_pp. Also, assuming that the power law is otherwise a good description of the black hole mass spectrum, we expect that the number of black holes in the PPSNe peak should be no more than the number of black holes that would have formed had the power law distribution continued to the upper limit of the upper mass gap, determined by the onset of pair-instability supernovae, m_PI ∼ 150 M_⊙. We impose this condition by requiring that the extrapolated area that would be under the power-law curve is less than the area contained within the Gaussian. This amounts to a restriction on the allowed values of λ,

$$\begin{eqnarray}&&\lambda \leqslant \displaystyle \frac{{\int }_{{m}_{\max }}^{{m}_{\mathrm{PI}}}{m}^{-\alpha }}{{\int }_{{m}_{\min }}^{{m}_{\mathrm{PI}}}{m}^{-\alpha }}=\displaystyle \frac{{m}_{\mathrm{PI}}^{1-\alpha }-{m}_{\max }^{1-\alpha }}{{m}_{\mathrm{PI}}^{1-\alpha }-{m}_{\min }^{1-\alpha }}\approx {\left(\displaystyle \frac{{m}_{\min }}{{m}_{\max }}\right)}^{\alpha -1}.\end{eqnarray} \tag{ 6 }$$

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_PI is the mass above which stars become PISNe, leaving no remnant. Here, we assume α > 1 and m_PI ≫ m_max.

To hone our intuition, we can plug in plausible values of α, m_min, m_max, and m_PI to determine a typical value of λ. For example, if we set m_min = 5 M_⊙, m_max = 50 M_⊙, and α = 2, we find λ ∼ 0.1. If we measure a peak consistent with these values, we anticipate that the position and width of the peak can inform our physical understanding of this mechanism. If λ is measured to be inconsistent with this constraint it could indicate either that the extrapolation of the power law is not a valid assumption or that the peak is not entirely due to PPSNe. We note that the fraction of observed black holes that formed through PPSNe will be larger than λ since the more massive black holes are observable out to a greater distance.

The smaller sensitive volume for low-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 it has a sharp cutoff 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 fiducial black hole mass spectrum with a smoothing function, S(m, m_min, δm), which rises from unity at m_min to unity at m_min + δm,

$$\begin{eqnarray}&&S(m,{m}_{\min },\delta m)={(\exp f(m-{m}_{\min },\delta m)+1)}^{-1}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&f(m,\delta m)=\displaystyle \frac{\delta m}{m}-\displaystyle \frac{\delta m}{m-\delta m}.\end{eqnarray} \tag{ 8 }$$

We note that δm = 0 recovers the step function used in previous analyses. Since the mass distribution is expected to be an increasing function, the peak of p(m₁) occurs below m_min + δm. We expect m_min ∼ 5 M_⊙ and δm ≳ 3 given that the black hole mass spectrum inferred from electromagnetic observations peaks at 8 M_⊙, with no black holes less massive than 5 M_⊙.

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

$$\begin{eqnarray}&&p({m}_{1}| {\rm{\Lambda }})=(1-\lambda ){p}_{\mathrm{pow}}({m}_{1}| {\rm{\Lambda }})+\lambda \,{p}_{{pp}}({m}_{1}| {\rm{\Lambda }})\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray*}&&{p}_{\mathrm{pow}}({m}_{1}| {\rm{\Lambda }})\propto {m}_{1}^{-\alpha }S({m}_{1},{m}_{\min },\delta m){ \mathcal H }({m}_{\max }-{m}_{1})\end{eqnarray*}$$

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

$$\begin{eqnarray*}&&{p}_{{pp}}({m}_{1}| {\rm{\Lambda }})\propto \exp \left(-\displaystyle \frac{{({m}_{1}-{m}_{{pp}})}^{2}}{2{\sigma }_{{pp}}^{2}}\right)S({m}_{1},{m}_{\min },\delta m).\end{eqnarray*}$$

encodes the peak from PPSNe.

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₁. This is motivated by observations of the stellar initial binary population (e.g., Kroupa et al. 2013; Belloni et al. 2017). In contrast to this, population synthesis models typically predict that the distribution of mass ratios should prefer equal-mass binaries (e.g., Belczynski et al. 2017). We model the distribution of the mass ratio as a power law with spectral index β as in Fishbach & Holz (2017) and Kovetz et al. (2017). For a distribution of mass ratio peaked at equal masses, β > 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,

$$\begin{eqnarray}&&p({m}_{2}| {m}_{1},{\rm{\Lambda }})={\left(\displaystyle \frac{{m}_{2}}{{m}_{1}}\right)}^{\beta }S({m}_{2},{m}_{\min },\delta m){ \mathcal H }({m}_{1}-{m}_{2}).\end{eqnarray} \tag{ 10 }$$

A table listing the (hyper)parameters and their physical meaning is provided in Table 1. Including a factor of V_obs to account for selection biases, the probability of detecting a mass pair given our (hyper)parameters, Λ = {α, m_min, m_max, δm, λ, m_pp, σ_pp, β} and under model H, is

$$\begin{eqnarray}&&\pi (m| {\rm{\Lambda }},H)\propto p({m}_{1}| {\rm{\Lambda }},H)p({m}_{2}| {m}_{1},{\rm{\Lambda }},H){V}_{\mathrm{obs}}(m).\end{eqnarray} \tag{ 11 }$$

We consider six different models for our (hyper)prior distribution, each corresponding to decreasingly stringent physical assumptions. These different prior assumptions can be tested with our Bayesian framework using Bayes factors. In our first model, H₀, we take the power-law distribution with maximum mass, m_max = 100 M_⊙. Since the injected data set considered in Section 4 has a minimum mass of 3 M_⊙, we allow m_min to vary, rather than fixing it at m_min = 5 M_⊙ as in previous analyses. In H₁ we introduce a uniform prior on m_max. In order to determine the relative importance of the different features, we switch to models that include all the effects described above except one. In H₂, H₃, and H₄ we do not include the Gaussian component, mass ratio, and low-mass smoothing respectively. Finally, in H₅ we include all of these effects. These prior choices are summarized in Table 2.

As with any phenomenological model, our model has limitations. If the proposed mass gaps exist in the population of black holes formed as the endpoint of stellar evolution there may still be black holes found in these gaps. The remnant of the binary neutron star merger GW170817 has mass M_rem ≲2.8 M_⊙ (Abbott et al. 2017d). It is not clear whether this object is a neutron star or a black hole. Similarly, the remnant from binary black hole mergers such as GW150914 is more massive than the suggested upper mass limit due to PPSNe, ${M}_{\mathrm{rem}}\,={62}_{-4}^{+4}\,{M}_{\odot }$ (Abbott et al. 2016f). Both of these objects lie within the proposed mass gaps. If either of these mergers happened in a dense environment such as a globular cluster, it is possible that such objects could merge with a new companion (Heggie 1975; Rodriguez et al. 2017). Similarly, primordial black holes (Hawking 1971) are not bound by the limitations of stellar evolution.

We do not expect either of these mechanisms to significantly affect the position and shape of an excess due to PPSNe, 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 ∼ 0.7 for equal-mass non-spinning pre-merger black holes (Scheel et al. 2009).

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

Given this, we expect to be able to identify an upper mass cutoff in the mass distribution after ≲20 events. We note that μ grows linearly with the 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 & Holz (2017).

After 200 events, the posterior distribution on λ, the fraction of black holes formed through PPSNe, 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.04}^{+0.07}$ (all future confidence regions will be 95% HDI unless specified) and disfavor λ = 0 at ≳3σ. Correspondingly, ${\mathrm{lnBF}}_{2}^{5}=4.9$, which is moderate evidence for the existence of the PPSN peak, but below our threshold for a confident detection.

Figure 3 shows the one- and two-dimensional posterior distributions for the position and width of the PPSN peak. We measure ${m}_{{pp}}={34.4}_{-1.2}^{+1.0}{M}_{\odot }$, ${\sigma }_{{pp}}={1.2}_{-1.2}^{+0.9}{M}_{\odot }$. We note that the peak and width of the distribution are covariant, with smaller values of m_pp requiring a larger σ_pp. This is unsurprising because the black holes with the highest mass, ∼40 M_⊙, must be accounted for. The posterior distribution on λ shows no significant correlation with m_pp or σ_pp.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

The (hyper)parameters describing the low-mass end of the distribution are more difficult to measure than the high-mass (hyper)parameters because the observation bias favors high-mass systems. After 200 events, there is no evidence for or against the low-mass smoothing described by δm. This is unsurprising since only nine injected binaries have m₁ ≲ 8 M_⊙, above which H₄ and H₅ are identical. Measuring the same events with improved strain sensitivity would not improve our sensitivity to δm. We are limited by cosmic variance: Cosmic ${\mathrm{lnBF}}_{4}^{5}=1$. The two-dimensional posterior distribution on m_min and δm, Figure 5, shows the correlation between low minimum masses and long turn-on lengths.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

In order to accommodate the lack of black holes with m ≳ 40 M_⊙, α is overestimated in model H₀. This leads to an overly steep inferred distribution and an overestimate of the total merger rate (see Section 4.6). Models H₁ and H₂, which do not include the Gaussian component, favor a mass spectrum that is less steep than the injected distribution. This is manifested as a positive gradient after accounting for observational bias. If we do not fit the power-law index of the mass ratio, we overestimate the number of heavy black holes. This bias is seen as an overestimate of λ in H₃ and maximum mass in H₁.

The majority of binary black hole mergers are not individually resolvable by Advanced LIGO/Virgo. Using a (hyper)parameterization that does not accurately describe the true distribution leads to a biased estimate of the fraction of mergers that are individually resolvable and hence the merger rate (Abadie et al. 2010; Abbott et al. 2017a). Compact binary coalescences are a Poisson process that can be described by a merger rate R(Λ). 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(Λ), and the observation time T:

$$\begin{eqnarray}&&R({\rm{\Lambda }})=\displaystyle \frac{N}{V({\rm{\Lambda }})T},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&V({\rm{\Lambda }})=\int d{\rm{\Theta }}\,\pi ({\rm{\Theta }}| {\rm{\Lambda }}){V}_{\mathrm{obs}}({\rm{\Theta }}),\end{eqnarray} \tag{ 13 }$$

and V_obs(Θ) is the volume sensitive to a given binary as introduced in Section 2.

To illustrate the dependence of R on the mass distribution model, we calculate the posterior distribution for the inferred estimate of the merger rate for the models described in Table 2. For each model, we compute the posterior distribution for R. During the first observing run of Advanced LIGO, ∼3 binary black hole mergers were identified in ∼48 days joint observing time (Abbott et al. 2016c). We use these values, i.e., N = 3 and T = 48/365 yr, to normalize our rate estimates. We neglect the (currently large) Poisson uncertainty in the arrival rate since this will be small once 200 detections have been made. Figure 7 shows the posterior distribution for the merger rate for each of our models. We note that if we assume that the power-law mass distribution extends out to 100 M_⊙ (the dashed–dotted line), we overestimate the merger rate by a factor of 2–3. This is due to the much larger α required to be consistent with the lack of detections at high masses. A similar result is obtained in Wysocki (2017) by simultaneously fitting the merger rate and power-law spectral index.

Zoom In Zoom Out Reset image size

Unresolvable mergers are widely believed to make the dominant contribution to the SGWB (Abbott et al. 2016g, 2017e, 2017f). The SGWB is typically characterized by the ratio of the energy density of the Universe in gravitational waves to the energy required to close the Universe, Ω_GW. The frequency of current detectors that is most sensitive to the SGWB is ∼30 Hz; this frequency corresponds to the inspiral phase of all binaries relevant to this work. The energy density due to binary black hole mergers depends on the distribution of chirp masses and the merger rate (Zhu et al. 2011),

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{\rm{GW}}}\sim \langle {{ \mathcal M }}^{5/3} \rangle R.\end{eqnarray} \tag{ 14 }$$

As seen above, cutting off the mass distribution around 40 M_⊙ leads to a reduction in the merger rate; however, this is accompanied by an increase in ${{ \mathcal M }}^{5/3}$. Overall, this leads to a reduction of ∼10% in the expected SGWB as seen in Figure 8. We also observe that relaxing the assumption that the secondary mass is uniformly distributed leads to a further ∼10% reduction in Ω_GW. This is because the chirp mass is maximized for equal-mass binaries for a given primary mass. These reductions are smaller than the current uncertainty on the amplitude of the background due to Poisson uncertainty in the observed merger rate. The current method of searching for this background is by cross-correlating the strain data from the two LIGO detectors; this method will take more time to resolve a weaker background.

Zoom In Zoom Out Reset image size

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 & 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 Ω_GW it will be more sensitive to the black hole mass function than cross-correlation searches.