Constraining the Black Hole Initial Mass Function with LIGO/VIRGO Observations

TL;DR

This study uses N-body simulations to analyze the initial mass distribution of black holes and compares the results with LIGO/Virgo gravitational wave data, suggesting a specific power-law spectrum with an upper mass cutoff.

Contribution

It introduces high-resolution dynamical simulations of black hole clusters to constrain the initial mass function using gravitational wave observations.

Findings

A power-law initial mass spectrum with index -2.35 is favored.

An upper mass cutoff around 50 solar masses fits the data.

The merger rate likely increases with redshift.

Abstract

Prior to the detection of black holes (BHs) via the gravitational waves (GWs) they generate at merger, the presence of BHs was inferred in X-ray binaries, mostly via dynamical measurements, with masses in the range between $\sim 5-20~M_\odot$. The LIGO discovery of the first BHs via GWs was surprising in that the two BHs that merged had masses of $35.6^{+4.8}_{-3.0}$ and $30.6^{+3.0}_{-4.4}\,M_\odot$, which are both above the range inferred from X-ray binaries. With 20 BH detections from the O1/O2 runs, the distribution of masses remains generally higher than the X-ray inferred one, while the effective spins are generally lower, suggesting that, at least in part, the GW-detected population might be of dynamical origin rather than produced by the common evolution of field binaries. Here we perform high-resolution N-body simulations of a cluster of isolated BHs with a range of initial mass spectra and upper mass cut-offs, and study the resulting binary mass spectrum resulting from the dynamical interactions. Our clusters have properties similar to those of the massive remnants in an OB association $\sim 10 \, \mathrm{Myr}$ after formation. We perform a likelihood analysis for each of our dynamically-formed binary population against the data from the O1 and O2 LIGO/Virgo runs. We find that an initial mass spectrum $M_{\rm BH}\propto M^{-2.35}$ with an upper mass cutoff $M_{\rm max}\sim 50M_\odot$ is favored by the data, together with a slight preference for a merger rate that increases with redshift.

Equations

$$p\left(M\mid \mathbf{d}\right)\propto p\left(\mathbf{d}\mid M\right)p\left(M\right).$$

$$p\left(\mathbf{d}\mid M\right)=\prod _{i=1}^{N_{\det }}p\left(d_i\mid M\right).$$

$$\frac{\mathrm{d}N}{\mathrm{d}V\mathrm{d}t}\propto {\left(1+z\right)}^{\lambda }.$$

$$\begin{gathered} p\left(d_i\mid M\right)= \\ \frac{\int \mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}zp\left(d_i\mid m_1,m_2,z\right)p\left(m_1,m_2,z\mid M\right)}{\int \mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}z{P}_{\det }\left(m_1,m_2,z\right)p\left(m_1,m_2,z\mid M\right)}. \end{gathered}$$

$$m_1,m_2,z\sim p\left(d_i\mid m_1,m_2,z\right)p_{\mathrm{PE}}\left(m_1,m_2,z\right),$$

$$\begin{gathered} \int \mathrm{d}{m}_1\mathrm{d}{m}_2\mathrm{d}zp\left(d_i\mid m_1,m_2,z\right)p\left(m_1,m_2,z\mid M\right) \\ \propto ⟨\frac{p\left(m_1,m_2,z\mid M\right)}{p_{\mathrm{PE}}\left(m_1,m_2,z\right)}⟩, \end{gathered}$$

Full text

gravitational waves — black hole physics — methods: numerical — binaries: general

††journal: ApJL

1 Introduction

The existence of black holes (BHs) is one of the primary predictions of the Theory of General Relativity. Prior to their direct discovery via the gravitational waves they generated in a merger event (Abbott et al., 2016a), their presence was inferred via dynamical mass measurements in X-ray binaries (see i.e. Wiktorowicz et al. 2014 for a summary). The values of the inferred masses vary between $\sim 4-5M_{\odot}$ to about 20 $M_{\odot}$, marking a clear separation with the inferred neutron star masses, for which the largest measurement to date is 1.96 $M_{\odot}$ (Demorest et al., 2010).

The discovery of the first binary black hole merger via the gravitational waves generated at the time of coalescence led to a mass measurement for the BH components of the merging binary: $35.6^{+4.8}_{-3.0}$ and $30.6^{+3.0}_{-4.4}\,M_{\odot}$. The large BH masses, both well above the maximum value measured to date in X-ray binaries, came as a surprise (Abbott et al., 2016b). The discovery triggered an intense debate in the literature on the formation pathway of this BH binary.

Broadly speaking, most formation avenues can be classified within one of two channels: isolated binary evolution, in which two massive stars evolve till their death while remaining gravitationally bound (e.g. Podsiadlowski et al., 2003; Dominik et al., 2013; Belczynski et al., 2016; Marchant et al., 2016; Mandel & de Mink, 2016), and dynamical formation by gravitational capture in dense environments, where binaries are being formed from isolated BHs as a result of frequent dynamical interactions (e.g. Sigurdsson & Hernquist, 1993; Portegies Zwart & McMillan, 2000; O’Leary et al., 2006; Miller & Lauburg, 2009; Mapelli et al., 2013; Leigh et al., 2014; Ziosi et al., 2014; Rodriguez et al., 2016; Antonini et al., 2016; Chatterjee et al., 2017; Rodriguez & Loeb, 2018; Samsing, 2018; Samsing & D’Orazio, 2018; Generozov et al., 2018; Antonini et al., 2018; Banerjee, 2018; Fragione & Kocsis, 2018; Fragione et al., 2018; Di Carlo et al., 2019; Ye et al., 2019). The theoretically-predicted rates are rather uncertain for both scenarios: the models explored by Belczynski et al. (2016) yield rates which vary between $\sim 6-1000$ Gpc*-3* yr*-1*. More recent, state-of-the art estimates of the rates of dynamical formation in Globular Clusters yield a range of $4-18$ Gpc*-3* yr*-1*(Rodriguez & Loeb, 2018). Both these rates are compatible with the current observationally-determined value by LIGO, of 9.7-101 Gpc*-3* yr*-1* (The LIGO Scientific Collaboration & the Virgo Collaboration, 2018). Both channels can in principle contribute to the observed population, something that can be tested since many more mergers are going to be detected in the future.

To date, after the first two observing runs of LIGO/Virgo, there have been 10 BBH mergers reported (The LIGO Scientific Collaboration & the Virgo Collaboration, 2018). While the smallest masses (4 out of 20) fall within the upper range of the masses inferred for X-ray binaries, the other 16 are all bigger, with the largest being $50.6^{+16.6}_{-10.2}M_{\sun}$. While the distributions are clearly not disjoint, there is a marked preference in GW-detected BHs for larger masses than in those found via X-ray binaries, raising the question of whether the two observed populations are dominated by the same progenitor population.

An independent piece of evidence which raises the same question is constituted by the measured spins111 What GWs measure is the so-called effective spin, i.e. the mass-weighed projection of the spins onto the orbital angular momentum. (see i.e. Farr et al. 2017; Roulet & Zaldarriaga 2019): generally high and aligned with the orbital angular momentum in the XRBs (especially the persistent ones), and typically low and isotropic in the GW-detected BHs. This trend is consistent with studies suggesting that, while isolated binaries preferentially yield BHs with spins aligned with the orbital angular momentum (i.e. Kalogera 2000), dynamically formed BHs have no preferred direction for alignment (i.e. Portegies Zwart & McMillan 2000).

In this paper we investigate the possibility that the observed BH population is dominated by a dynamical formation channel complementary to the dynamical channel discussed in the context of globular clusters. More specifically, we generate small clusters of BHs with properties appropriate for the massive remnants of an OB association $\sim 10\,\mathrm{Myr}$ after its formation, following the supernovae explosions of its massive constituents. Therefore, albeit dynamically formed, the binary population discussed here does not necessarily need to originate in globular clusters, but could also be associated with field stars. Using high-resolution N-body simulations of these clusters with a range of mass spectra, we explore the dependence of the BH binary mass distribution on the mass spectrum and upper mass cut-off of the individual BHs (§2). We perform a statistical comparison with the data from the O1 and O2 LIGO/Virgo runs to study the consistency between the data and the simulation results, and assess statistical preferences towards an initial BH mass spectrum (§3). We summarize and conclude in §4.

2 The BBH mass spectrum from dynamical interactions

Motivated by the considerations of Sec. 1, here we set to perform high-resolution N-body simulations of a cluster of BHs, with the goal of exploring the dependence of the mass spectrum and the orbital parameters of the dynamically-formed binaries on the mass spectrum of the isolated BH population. We note that the mass spectrum of BH binaries produced as a result of dynamical interactions in clusters has a long history in the literature, predating the era of gravitational waves. A preferential tendency for dynamically formed binaries to have heavier mass has been noted in a number of works (O’Leary et al., 2006; Miller & Lauburg, 2009; Ryu et al., 2016; Rodriguez et al., 2018; Di Carlo et al., 2019). Our fewbody simulations, while not including the effects of the cluster potential, they allow us to accurately follow binary formation, since in dense environments this is dominated by 3-body scatterings, as well as perform a large number of Monte Carlo realizations.

We consider a cluster of 20 BHs, with an initial binary fraction equal to zero in order to purely explore the binary properties due to dynamical formation. Their positions in the cluster are initially distributed randomly in a sphere of radius 0.1 pc. Astrophysically, this configuration can be thought of as representing the remnants of an OB association (typically comprised of $\sim 10-100$ stars), still confined within the dense core of a molecular cloud222 Dynamically, this population however bears resemblance with a core-collapsed cluster (see e.g. Sigurdsson & Hernquist 1993; Rodriguez et al. 2016). (e.g. Zhou et al. 1994). From a numerical point of view, we note that the particular number of 20 was chosen as a ’sweet spot’: large enough for obtaining a reasonably well sampled binary mass distribution, but small enough to enable the running of a large number of realizations with high numerical accuracy. However, we ran several additional simulations with different numbers of BHs in order to verify that the shape of the binary distribution remains statistically the same as the number of BHs is varied. Similarly, we chose the size of the initial spatial domain after verifying that it was large enough that the results for the binary mass distributions were converged as the region size was varied.

The BHs were assigned a mean speed of 5 km/s, as typical of the velocity dispersions observed in low-mass star clusters (e.g. Harris, 1996). We follow their evolution using our code SpaceHub (see Wang et al. 2019 for details), which employs the ARCHAIN algorithm (Mikkola & Merritt, 2008) to accurately trace the motion of tight binaries with arbitrarily large mass ratios and eccentricities, and a chain structure to reduce the round-off errors from close encounters. Binaries are detected in the simulations when the following conditions are verified: (a) the BH masses $m_{1}$ and $m_{2}$ are gravitationally bound to one another; (b) the system $(m_{1}+m_{2})$ is gravitationally unbound in the potential of the remaining BHs; (c) the binary has traveled a large enough distance from the original BH cluster. We chose the distance to be 20 times the size of the original cluster, after verifying that the simulation results are converged for that value. The initial BH masses are drawn from a distribution with $M_{\rm min}=5M_{\odot}$ and $M_{\rm max}$ varying within the range $40-50\,M_{\odot}$ in steps of $5\,M_{\rm\odot}$. The mass spectrum of BH remnants, in addition to depending on the evolutionary model, is also strongly dependent on metallicity, varying from an almost flat distribution at solar metallicities to an almost linear dependence on the main sequence mass at low metallicities and for high-mass progenitor stars (e.g. Spera et al. 2015). Therefore, we explored a variety of mass spectra, ranging from a flat distribution to a power-law $M^{-\alpha}$ with index $\alpha=4$ in steps of 0.5. Additionally, we also investigated the particular case of $M^{-2.35}$, corresponding to a BH mass spectrum reflective of the initial mass function of the massive progenitor stars (Salpeter, 1955).

We show some representative results in Fig. 1 ($\alpha=1$), Fig. 2 ($\alpha=2$), and Fig. 3 ($\alpha=3$), all with $M_{\rm max}=50\,M_{\odot}$. In the left panels we show the mass ratio $q$, while in the right panels the merger time $\tau_{\rm GW}$ calculated according to Peters’ formulae (Peters, 1964), both versus the total mass of the binaries. The flatter mass spectrum $(\alpha=1)$ yields a large fraction of of binaries with total mass $\sim 40-65\,M_{\odot}$. Interestingly, not only is there an exceedingly large number of very massive binaries formed, but they are also the ones which are more tightly bound, hence resulting in shorter merger times due to GW emission (note a similar result found by Ryu et al. (2016) in the context of the formation of the first X-ray binaries). The most massive binaries also tend to have higher mass ratios. Note that the kink at $\sim 55\,M_{\odot}$ corresponds to the mass of the binary formed by the most and the least massive BHs, which results in the mass distribution to undergoe a change of slope as it passes through that point.

As the mass spectrum steepens to $M^{-2}$, the distribution becomes more apparently dominated by lower-mass binaries, but a tail in the high mass range still remains. For $\alpha=3$, the binary mass distribution becomes clearly peaked towards low masses, and the high-mass tail of the distribution becomes vanishingly small. The fraction of BHs which end up as ejected binaries is also dependent on the mass spectrum. In particular, we found this fraction to be 0.11% for $\alpha=1$, 0.09% for $\alpha=2$ and 0.065% for $\alpha=3$. Correspondingly, the relative fraction of those binaries which merge within the Hubble time is 0.092%, 0.069%, and 0.048% for $\alpha=1,2,3$, respectively.

Of particular relevance for the LIGO/Virgo results is the fact that the more massive binaries tend to be the more tightly bound; the heaviest objects are in fact the ones with the largest cross-sections for encounters, hence they undergo the most scatterings and end up as the hardest binaries. This tendency is especially pronounced for shallower BH mass spectra, when the number of massive BHs is not much smaller than the number of lighter BHs. For very steep slopes of the mass spectrum, the interaction probability becomes dominated by the number of small BHs, which hence have much higher chances of interacting, and thus forming tight binaries.

3 Statistical comparison with the LIGO/Virgo data from the O1/O2 runs

We now wish to compare the models described in the previous section to the 10 binary black hole mergers observed by LIGO and Virgo in GWTC-1 (The LIGO Scientific Collaboration & the Virgo Collaboration, 2018)333Our analysis can be found at https://github.com/farr/ClusterBHGW.. The analysis is very similar to Farr et al. (2017): here we have a collection of zero-parameter models that predict the mass distribution of merging black holes. Unlike in Farr et al. (2017), the detectability of mergers is a strong function of mass, so we must account for selection effects; see Mandel et al. (2018) and references therein.

We are comparing to a dataset, $\mathbf{d}$, consisting of a catalog of detections, $i=1,\ldots,N_{\mathrm{det}}=10$. The Bayseian posterior probability of a particular mass model, $M$, is given by

[TABLE]

$p\left(\mathbf{d}\mid M\right)$ is sometimes known as the "Bayes factor;" it is the likelihood of model $M$ given the observed data. $p(M)$ is the prior probability of model $M$, which we are free to assign based on our experience and intuition; we describe our model priors below as we discuss Figure 4.

We assume that the noise realization in the LIGO and Virgo detectors is statistically independent for each event, so that

[TABLE]

Each model makes predictions about the masses of the merging binaries. In principle each model also makes predictions about the redshift distribution of merging binaries, but we leave study of this prediction to future work. Instead, we select only the mergers whose time to merger is $t_{\mathrm{GW}}<10^{10}\,\mathrm{yr}$ and, following Fishbach et al. (2018), impose a parameterized redshift distribution of events corresponding to a volumetric merger rate in the comoving frame of

[TABLE]

Setting $\lambda=3$ gives a merger rate that approximately tracks the star formation rate; setting $\lambda=0$ gives a merger rate that is constant in the comoving frame (Fishbach et al., 2018). The likelihood of the data depends on the masses and redshifts of the merging systems, which are subject to selection effects, so we have (Farr et al., 2017; Mandel et al., 2018)

[TABLE]

Here the numerator is the likelihood of the LIGO data given the masses and redshifts predicted by the model $M$, and the denominator is the correction for the selection function and gives the average detectability for the model population.

The denominator is independent of the data, $d_{i}$, and common to all events. We use a Monte-Carlo estimate of the integral (Farr, 2019) obtained by generating synthetic merger events and detecting them using an analytic estimate of the LIGO/Virgo O1+O2 sensitivity (Abbott et al., 2016c).

The numerator can also be computed via Monte-Carlo using parameter estimation samples from The LIGO Scientific Collaboration & the Virgo Collaboration (2018). Those are drawn from a posterior density that incorporates the likelihood and a prior, $p_{\mathrm{PE}}\left(m_{1},m_{2},z\right)$ (where ’PE’ stands for parameter estimation)

[TABLE]

so the likelihood integral that we need can be computed via

[TABLE]

where the final average is taken over the PE samples. We use a Gaussian kernel density estimator in a two-dimensional, un-constrained parameter space $\left(x,y\right)=\left(\log m_{1},\log m_{2}-\log\left(m_{1}-m_{2}\right)\right)$ to smooth the distribution over masses predicted by each model when computing the model likelihood.

We begin our analysis by considering, for each of the models described in Sec. 3, a population drawn from a redshift distribution varying from constant ($\lambda=0$), to a rapidly evolving one ($\lambda=6$). We impose uniform prior density in models in $\alpha$ and $\lambda$. The 2D posterior on models in $\alpha$-$\lambda$ space is displayed in Fig. 4, together with the 1D projections on the $\lambda$ and $\alpha$ axis, for the case with $M_{\rm max}=50\,M_{\odot}$, which appears to provide the best match.

The analysis shows that there is a slight preference for an evolving redshift distribution, though the trend is only marginal. This is consistent with the results on redshift evolution from The LIGO Scientific Collaboration et al. (2018). Marginalizing over $\lambda$, the preferred spectral index is found to be $\alpha=2.35$. A rough 1$\sigma$ (68% credible) interval for $\alpha$ is $\alpha=2.35^{+0.55}_{-0.36}{}$.

In Figure 5 we compare our models’ predictions for the distribution of the observed total mass of 10 merging binaries for various values of $M_{\rm max}$, fixing $\lambda=3$ (i.e. a merger rate that tracks the star formation rate), and $\alpha=2.35$. The predictions incorporate both the GW selection function and also observational uncertainties. No model fully reproduces the distribution of observed total masses, though the observations lie within the 1$\sigma$ band of the $M_{\rm max}=50M_{\odot}$, $\alpha=2.35$ model. We remark that our models are minimally parameterized, being dependent only on the two parameters $\alpha$ and $M_{\rm max}$; the shape of the binary mass distribution is then determined by dynamics alone, and hence its resemblance to the observed distribution for some astrophysically interesting sets of parameters ($\alpha=2.35$ and $M_{\rm max}=50\,M_{\odot}$) is especially intriguing.

4 Summary

The discovery of BHs via the GWs emitted when they merge in a binary has confirmed one of the milestone predictions of the Theory of General Relativity, while at the same time opening a new window into our exploration of the Universe. As often happens with new observations, this new window has also raised some questions. In the case of the 20 BHs discovered by the LIGO/Virgo collaboration via their mergers in binaries during the O1/O2 runs, their mass spectrum has been somewhat surprising, being shifted towards larger masses with respect to the BHs whose masses had previously been measured dynamically in X-ray binaries. The measured spins, on the other hand, have been found to be mostly consistent with low and isotropic, unlike the generally high ones measured in X-ray binaries.

In this Letter we have investigated the possibility that X-ray and GW-detected BHs are dominated by different formation channels: isolated binary evolution for the former and (primarily) dynamical formation for the latter. Via high-resolution N-body simulations of a mini-cluster of initially isolated BHs (which can be thought of as the remnants of an OB association), we have shown a tendency for binary BH formation among the heaviest objects in the cluster (see also O’Leary et al. 2006; Rodriguez & Loeb 2018; Di Carlo et al. 2019 for similar trends in globular clusters). The heaviest BH binaries tend to also be the ones which are more tightly bound, hence resulting in shorter merger times, which enhances the probability of being detected via GWs.

While weighed towards larger masses, the precise shape of the mass distribution of the dynamically-formed binary BHs is also reflective of the particular initial BH mass function. We investigated this distribution for a variety of BH mass spectra, from a flat distribution to a powerlaw with index of -4 and extracted the sub-population of dynamically-formed binaries which merge within a Hubble time. We hence performed a Bayesian statistical analysis to compare the likelihood of each of these models to the LIGO/Virgo data from the O1 and O2 observing runs. We found that an initial BH mass spectrum $\propto M^{-2.35}$ is favored by the data, together with a maximum BH mass $M_{\rm max}\sim 50M_{\odot}$. This is consistent with the theoretical upper limit for stellar BH masses, which is set by the occurrence of Pair Instability (Woosley, 2017; Marchant et al., 2018).

A slope $\propto M^{-2.35}$ reflects the initial mass function of massive stars, which is expected at the low metallicities required to form very massive BH remnants (i.e. Spera et al. 2015). Hence our work shows that dynamically formed binaries from low-metallicity stars are reasonably compatible with the binary BH mass distribution from the O1 and O2 LIGO/Virgo observing runs.

As more data is expected to be gathered in the years to come, statistical comparisons with numerical simulations will allow to establish whether the dynamical formation channel is indeed dominant, and to reconstruct the mass spectrum of the initial BHs, thus shedding a new light on massive stars and their evolution.

We thank Johan Samsing for valuable comments on our manuscript. RP acknowledges support from the NSF under grant AST-1616157. The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.