Unraveling the origin of black holes from effective spin measurements with LIGO-Virgo

TL;DR

This paper uses LIGO-Virgo gravitational wave data to analyze the effective spin of binary black hole mergers, aiming to distinguish between primordial and astrophysical origins through statistical modeling and future event predictions.

Contribution

It introduces a method to differentiate black hole origins by analyzing effective spin distributions and odds ratios from gravitational wave observations.

Findings

Primordial black holes have a distinct effective spin distribution.

Current data can constrain the mixture of black hole origins.

Future detections will improve discrimination between models.

Abstract

We investigate how to use information on the effective spin parameter of binary black hole mergers from the LIGO-Virgo gravitational wave detections to discriminate the origin of the merging black holes. We calculate the expected probability distribution function for the effective spin parameter for primordial black holes. Using LIGO-Virgo observations, we then calculate odds ratios for different models for the distribution of black holes' spin magnitude and alignment. We evaluate the posterior probability density for a possible mixture of astrophysical and primordial black holes as emerging from current data, and calculate the number of future merger events needed to discriminate different spin and alignment models at a given level of statistical significance.

Tables

Table 1. Table 1 : Selected parameters of the ten BBH mergers events detected during LIGO’s O1 and O2 runs. The parameters are median values, with 90% credible intervals [ 33 ] .

Event		$m_1\left[M_{\odot }\right]$		$m_2\left[M_{\odot }\right]$		${\chi }_{\mathrm{eff}}$
GW150914		${35.6}_{-3.0}^{+4.8}$		${30.6}_{-4.4}^{+3.0}$		$-{0.01}_{-0.13}^{+0.12}$
GW151012		${23.3}_{-5.5}^{+14.0}$		${13.6}_{-4.8}^{+4.1}$		${0.04}_{-0.19}^{+0.28}$
GW151226		${13.7}_{-3.2}^{+8.8}$		${7.7}_{-2.6}^{+2.2}$		${0.18}_{-0.12}^{+0.20}$
GW170104		${31.0}_{-5.6}^{+7.2}$		${20.1}_{-4.5}^{+4.9}$		$-{0.04}_{-0.20}^{+0.17}$
GW170608		${10.9}_{-1.7}^{+5.3}$		${7.6}_{-2.1}^{+1.3}$		${0.03}_{-0.07}^{+0.19}$
GW170729		${50.6}_{-10.2}^{+16.6}$		${34.3}_{-10.1}^{+9.1}$		${0.36}_{-0.25}^{+0.21}$
GW170809		${35.2}_{-6.0}^{+8.3}$		${23.8}_{-5.1}^{+5.2}$		${0.07}_{-0.16}^{+0.16}$
GW170814		${30.7}_{-3.0}^{+5.7}$		${25.3}_{-4.1}^{+2.9}$		${0.07}_{-0.11}^{+0.12}$
GW170818		${35.5}_{-4.7}^{+7.5}$		${26.8}_{-5.2}^{+4.3}$		$-{0.09}_{-0.21}^{+0.18}$
GW170823		${39.6}_{-6.6}^{+10.0}$		${29.4}_{-7.1}^{+6.3}$		${0.08}_{-0.22}^{+0.20}$

Table 2. Table 2 : Natural log Bayes factors for various spin distributions with q = 1 𝑞 1 q=1 . Right : Values reported by LIGO [ 2 ] . Left : Values found in this work.

LIGO	Low	Flat	High
Isotropic	$0.0$	$-0.93$	$-2.07$
Aligned	$-4.12$	$-12.92$	$-32.37$

Equations

$$\chi =\frac{|\vec{S}|}{G{m}^2},$$

$${\chi }_{\mathrm{eff}}=\frac{m_1{\chi }_1\cos {\theta }_1+m_2{\chi }_2\cos {\theta }_2}{m_1+m_2},$$

$$p(\chi )\approx \exp \left[-\frac{{\chi }^2}{2{\sigma }^2}\right].$$

$$p(\chi ;\mu ,\sigma )=\frac{\mathcal{N}(\mu ,\sigma )}{\sqrt{2\pi }\sigma }\exp \left[-\frac{(\chi -\mu {)}^2}{2{\sigma }^2}\right],$$

$$p\left(\{d^i\}\mid \Lambda \right)=\prod _{i=1}^{N_{\mathrm{obs}}}p\left(d^i\mid \Lambda \right),$$

$$p\left(d^i\mid \Lambda \right)=\int \mathrm{d}{\chi }_{\mathrm{eff}}^{i}p\left(d^i\mid {\chi }_{\mathrm{eff}}^i\right)p\left({\chi }_{\mathrm{eff}}^i\mid \Lambda \right)$$

$$p\left(\Lambda \mid \left\{d^i\right\}\right)\propto \left[\prod _{i=1}^{N_{\mathrm{obs}}}\int \mathrm{d}{\chi }_{\mathrm{eff}}^{i}p\left(d^i\mid {\chi }_{\mathrm{eff}}^i\right)p\left({\chi }_{\mathrm{eff}}^i\mid \Lambda \right)\right]p\left(\Lambda \right).$$

$${\chi }_{\mathrm{eff}}=\frac{{\chi }_1^z+q{\chi }_2^z}{1+q},$$

$$p({\chi }_i^z)={\int }_0^1{f}_{{\chi }_i}({\chi }_i)\mathrm{d}{\chi }_i{\int }_{-1}^1{f}_{\cos {\theta }_i}(\cos {\theta }_i)\delta ({\chi }_i^z-{\chi }_i\cos {\theta }_i)\mathrm{d}\cos {\theta }_i,$$

$$p({\chi }_i^z)=\frac{1}{2}{\int }_{|{\chi }_i^z|}^1{f}_{{\chi }_i}({\chi }_i^z/\cos {\theta }_i)\frac{1}{|\cos {\theta }_i|}\mathrm{d}\cos {\theta }_i.$$

$$p\left({\chi }_{\mathrm{eff}}\right)$$

$$=(1+q){\int }_a^{b}p\left({\chi }_2^z\right)p\left((1+q){\chi }_{\mathrm{eff}}-q{\chi }_2^z\right)\mathrm{d}{\chi }_2^z,$$

$$a$$

$$b$$

Full text

1 Introduction

After the detection of binary black hole (BBH) merger events with LIGO-Virgo [1, 2], the question of the physical origin of the black holes has become somewhat pressing. In particular, there has been some significant interest in the possibility that some or all of the BBH events originate from primordial black holes (PBH) [3, 4], i.e., black holes originating from large over-densities in the very early universe rather than from the collapse of stellar objects (see e.g. [5, 6]). Interestingly, this interpretation is compatible with the notion that these PBH could be the dark matter needed for a consistent picture of the early and large-scale structure of the universe (see e.g. [3]). Whether or not this possibility is ruled out is the focus of intense debate, the key issues at stake including (1) the problem of matter accretion that might produce significant-enough accelerated charged particles to perturb in a measurable, and excessive way the cosmic microwave background photons [7, 8] (see however [9]); (2) the problem of disruption of small-scale structure by the relatively massive black holes that would make up the dark matter [10, 11] (see however [12]); and (3) limits from gravitational lensing of type Ia supernovae [13] (see however [14]).

These constraints notwithstanding, it is known (and it is reviewed below) that PBH that were produced during a radiation-dominated cosmological epoch have low intrinsic spin magnitude. As a result, a generic prediction of the picture where the LIGO-Virgo BBH events are in part or all PBH mergers is that the effective spin parameter (to be defined below, and which depends among other things on the intrinsic black hole spin magnitudes) is very low. Incidentally, it is important to point out that this fact (that binary mergers of PBH should have a low effective spin parameter) does depend on cosmology: it has been shown, for instance, that if the universe went through a matter-domination phase during PBH formation, then in fact the intrinsic spin of the resulting PBHs is generically close to maximal [15]. Furthermore, other possibilities might arise for cosmologies where the early universe was neither matter- nor radiation-dominated at early times, as we explored recently [16, 17]. With this caveat in mind, we shall hereafter assume radiation domination at the time of PBH formation.

The key motivation for the present study is that, so far, 9 out of 10 of the LIGO-Virgo BBH mergers are compatible with very low effective spin parameters. It seems timely, therefore, to assess what the predicted probability density for PBH’s effective spin parameter is, and to compare it with current observations; additionally, we intend to explore how many events it will take to differentiate in a statistically robust fashion different effective spin parameter models. Although at present there is no firm prediction for astrophysical black holes’ effective spin parameter, and there is debate concerning selection and bias effects in the detected BBH events, our study intends to point out that (i) current observations of the effective spin parameter are largely compatible with a dominant PBH component in the BBH mergers, and that (ii) in the future, the effective spin parameter distribution could help discriminating PBH from “ordinary”, astrophysical black holes.

Aspects of the question we address here have been considered in the recent literature, with different assumptions and methods. Soon after the detection of the first four BBH mergers [1, 18, 19, 20, 21], ref. [22] argued that information on the effective spin could be used to distinguish between aligned versus isotropic angular distributions; specifically, [22] concluded that as long as the black hole spin values are not intrinsically small (which, however, might well be a distinct possibility, as pointed out e.g. by ref. [23]) then an aligned angular distribution is strongly disfavored. Additionally, ref. [22] also showed how with relatively few additional events, the odds ratio would conclusively point in one direction or another (i.e., isotropic or aligned).

Ref. [24] reiterated how scenarios considered in the formation of astrophysical black hole binaries naturally lead to isotropic (for dynamical capture) or near-aligned (for common envelope evolution) black hole spins. They also showed that Bayesian statistics allows one to distinguish, at a given confidence level, which fraction of the binaries are preferentially aligned versus isotropically distributed. A similar analysis was conducted in ref. [25], with the additional points that a discrimination between isotropic and aligned spin distributions might be possible even regardless of the intrinsic spin magnitude distribution. They also showed that once an aligned or isotropic spin distribution is established, it is possible to re-construct the spin magnitude distribution with a high degree of confidence. The possibility of disentangling the existence of sub-populations of binary black holes with different spin orientations was explored in ref. [26, 27], which also pointed out how a “pure” distribution would be statistically preferred with relatively few events (see also ref. [28]).

In this study, we utilize the intrinsic spin distribution of PBHs from the results of ref. [29] to calculate (to our knowledge for the first time) the predicted prior probability distribution for the effective spin parameter. Ref. [29] assumes that there is no correlation between the overdensity leading to the formation of the PBH and its spin, and that the probability density for the spin distribution as a function of the overdensity is flat. A recent study, ref. [30], challenges these assumptions, and finds a peaked distribution for the PBH spin parameter, which critically depends on the width of the power spectrum peak giving rise to the PBH, and on the relative abundance of PBH. As a result, the spin distribution is significantly narrower than what was predicted in ref. [29] (see also refs. [31, 32]). In what follows, we compare the resulting spin distribution for a variety of assumptions for PBH formation as outlined in ref. [30], and compare it to our benchmark choice which reflects the results of ref. [29] (and which, as we explain below, can be seen as a limiting case of the setup of ref. [30]).

To these ends, in this study we first explore, in sec. 2, the theoretical prediction for the effective spin parameter distribution for PBH. We outline the assumed astrophysical black holes spin magnitude distribution we consider, and we review LIGO-Virgo observations; in sec. 3 we compare the odds ratios for the models we consider and study the favored “mixture” of different such models; we then forecast how future events will inform both the odds ratios and the inference of the relative level of mixing of different models. Finally, in sec. 4 we discuss our results and present our conclusions.

2 Effective Spin Distribution

The spin magnitude of a Kerr black hole (BH) is commonly defined via a dimensionless spin parameter $\chi$,

[TABLE]

where $\vec{S}$ and $m$ are the spin and mass of the BH, respectively. One of the most important parameters that LIGO can infer from the gravitational waveform is the effective spin parameter $\chi_{\rm eff}$, defined as:

[TABLE]

where $\theta_{i}=\cos^{-1}(\vec{L}\cdot\vec{S_{i}})$ is the tilt angle between the spin $\vec{S}_{i}$ and the orbital angular momentum vector $\vec{L}$. As apparent from its definition, the parameter $\chi_{\rm eff}$ is a quantity sensitive to both the spin alignment of the two black holes with their orbit (angular momentum of the binary) before the merger, and to the magnitude of the individual spins. $\chi_{\rm eff}$ is a dimensionless number ranging from $-1$ to $1$, where for $\chi_{\rm eff}=1$ the spins of both black holes are perfectly aligned with their orbit, and $\chi_{\rm eff}=-1$ the spins are perfectly anti-aligned. Values of $\chi_{\rm eff}\approx 0$ can stem from one or both of the following physical situations: (i) the black hole intrinsic spins are anti-aligned with each other, or (ii) the magnitude of the intrinsic effective spin parameters, $\chi_{i}\ll 1$. There could be a third possibility, that both spins are perpendicular to the orbit, but this is somewhat less likely and less physically motivated.

First, let us gain some intuition about the effective spin distribution for a few binary black hole (BBH) formation channels. One possible formation channel for BBH is from massive isolated binaries through common envelope evolution, where the intrinsic spin is generally aligned along the same direction as the orbital angular momentum, meaning that $\chi_{\rm eff}\approx 1$. On the other hand, there exist dynamical scenarios where we expect most BBHs to have spins largely uncorrelated with their orbit meaning that $\chi_{\rm eff}\approx 0$. This is the case, for instance, for BBHs formed dynamically in dense stellar environments, and it is also the case for PBHs, which additionally are predicted to have small intrinsic spins. It is important to notice that a key consequence, and possible signature, of any isotropic formation mechanism is that the distribution of $\chi_{\rm eff}$ is symmetric around zero, regardless of the spin magnitude distribution [25].

2.1 LIGO-Virgo effective spin measurements

We list in Table 1 the relevant observed properties of the 10 BBH merger events we consider in our study: the masses of each individual black hole $m_{1}$ and $m_{2}$ (columns 2 and 3), and the corresponding dimensionless effective spin $\chi_{\rm eff}$ (column 4). It is worth pointing out that before the first gravitational wave detection, LIGO was expecting 33-100 more NS-NS binary events compared with BH-BH binaries [34] (see however [35, 36] for different earlier predictions closer to what currently observed); however, LIGO’s O1 and O2 run showed that the rate of BH-BH binaries is an order of magnitude greater than the NS-NS binaries. Furthermore, the range of black hole masses was expected to be from $\sim 5\mathrm{M}_{\odot}$ to $\sim 15\mathrm{M}_{\odot}$ [37]. As a result, the first detection (GW150914) came somewhat as a surprise, because previously-known black holes were significantly lighter than the inferred masses, among various reasons.

As evident in Table 1, the majority of black holes are over $25\mathrm{M}_{\odot}$ with the heaviest being $50\mathrm{M}_{\odot}$ (GW170729). The nature of this new population of heavy stellar-mass black holes is still debated in the literature [38, 39, 40, 41]; notice that a reason for observing more massive binary mergers over lighter ones could partly be explained as a selection bias [42], since more massive BBH mergers produce a louder signal, and therefore the accessible space-time volume is larger than for lighter systems. It has been proposed that there exists a mass gap ($\sim 50\mathrm{M}_{\odot}-150\mathrm{M}_{\odot}$) due to pair-instability supernovae for stellar black holes [43, 44, 38] and therefore, that the black hole masses cannot be arbitrary large. There have also been claims of a cutoff at high masses in the current detections made by LIGO to date [45, 46, 47, 48].

The last column of Table 1 shows the most interesting parameter for this work: the effective spin $\chi_{\rm eff}$, defined by Eq. (2.2). The listed 10 observed events appear to disfavor high spin magnitude aligned with the orbital angular momentum, unlike the large spin values (near to the maximum possible value) observed in the majority of black holes in X-ray binaries [49, 50]. Most events are consistent with $\chi_{\rm eff}=0$, with two exceptions: GW170729 and GW151226. These two events show evidence of positive, but relatively small, $\chi_{\rm eff}$ values. Fig. 1 shows the posterior distribution of $\chi_{\rm eff}$ for the ten events observed by LIGO’s O1 and O2 run and the prior assumed by the LIGO Collaboration111https://dcc.ligo.org/LIGO-P1800370/public.

It is important to note that the LIGO Collaboration used Bayesian statistics to analyze the data and to infer the source properties of all ten BBH gravitational wave events [33]. This means that one needs to properly define prior probability density distributions. While ideally the conclusions should be robust and fairly independent under the choice of different priors, if the data are only mildly informative, priors could influence the statistical inference on the source parameters (see ref. [51] for a discussion on this point). An analysis of the importance and effect of the choice of priors on the first three LIGO events has been carried out in ref. [24]. One should also bear in mind that there certainly exist selection bias effects for the posterior distribution of $\chi_{\rm eff}$. For instance, sources with positive $\chi_{\rm eff}>0$ have a more clear signal, due to longer time orbiting before merging, therefore allowing to better constrain the waveform (see e.g. refs. [52, 53]).

2.2 PBH spin distribution

Given a spin distribution for the intrinsic spin of individual primordial black holes, and the assumption of isotropy in the spin-orbit alignment, one can calculate the distribution of $\chi_{\rm eff}$. We follow here ref. [29] in assuming that the distribution function for the intrinsic spin magnitude of a single PBH can be closely approximated by the Gaussian functional form

[TABLE]

The parameter $\sigma$ is, in principle, calculable given the spectrum of density perturbations leading to PBH formation in the early universe. Absent this, and in view of the fact that there might be circumstantial evidence for a PBH origin of at least some of the detected BBH merger events [54], here we infer the value of $\sigma$ directly from observations. To this end, we define a probability distribution for $\chi$ as a half Gaussian

[TABLE]

entertaining a possible non-zero value for $\mu$. In Eq. (2.4) $\mathcal{N(\mu,\sigma)}$ is the appropriate normalization constant. Notice that the function is only defined in the interval $\chi\in[0,1]$ and zero otherwise (see the definition of $\chi$ in Eq.(2.1)). Our goal is to investigate the probability distribution of $\chi_{\rm eff}$ resulting from an isotropic spin orientation distribution and the intrinsic spin distribution of Eq. (2.4) and to infer the posterior probability density for the parameters $\mu$ and $\sigma$ from the LIGO-Virgo data.

We utilize here a hierarchical Bayesian analysis: we assume that the individual BHs in the binary are coming from a primordial population which is described by some hyper-parameters $\Lambda$. Then, we use the 10 LIGO-Virgo events to derive the posterior distribution for $\Lambda$. Our approach is analogous to refs. [22, 25, 46, 55]. Our goal is thus to find $p\left(\Lambda\mid d\right)$, the probability distribution of the parameters $\Lambda$ given the data $d$, where $\Lambda$ simply represents here the parameters that describe the PBH population $\Lambda=\left\{\mu,\sigma\right\}$.

Assuming that the events are independent of each other, we can combine the individual likelihoods to build a joint likelihood

[TABLE]

where

[TABLE]

is the likelihood function for $i$th event and $p\left(d\mid\chi_{\rm eff}\right)$ is the marginal likelihood, meaning that it has been marginalized over all parameters but $\chi_{\rm eff}$. Since the LIGO-Virgo collaboration have not yet released the marginalized or full likelihoods to the general public, but have rather provided the posterior distributions which we show in Fig. 1, we need to re-weight the posterior distribution of $\chi_{\rm eff}$ to obtain the likelihoods. The last term in the integral, $p\left(\chi_{\rm eff}\mid\Lambda\right)$, is the probability of measuring $\chi_{\rm eff}$ given the parameters of our model $\Lambda$. In our case this distribution of $\chi_{\rm eff}$ has been derived in App. A.1 and the result is given in Eq. (A.4). Finally, using Bayes’ theorem we obtain the posterior distribution of the parameters $\Lambda$ as

[TABLE]

Here, $p\left(\Lambda\right)$ is the prior choice for the parameters $\Lambda$, which we take here to be uninformative (i.e. we use a flat prior for both $\mu$ and $\sigma$).

The posterior distribution for the hyper-parameters $\mu$ and $\sigma$ describing the putative PBH population given the 10 LIGO-Virgo events is shown in the left panel in Fig. 2. The distribution for $\mu$ is almost flat, except for $\mu>0.15$, a slightly disfavored range of values. In contrast, the distribution for $\sigma$ offers more information: one can clearly see a peak around $\sim 0.3$. Notice that $\mu$ and $\sigma$ are anti-correlated, as expected: given a half Gaussian with zero mean and fixed width, one can find an approximately equivalent distribution with a negative mean and larger width; conversely, a distribution with a positive mean will correspond to one with mean zero and a narrower width.

Finally, in the right panel of Fig. 2 we show the posterior distribution of $\sigma$ when we set $\mu=0$. This choice is motivated by the analytical findings of ref. [29] (see their Eq.(17) and figure 2), and is approximately valid even in the scenario discussed in ref. [30], where the preferred value of $\mu$ is non-zero but extremely small (see their Eq.(8.7) and Fig.7).

Based on the 10 LIGO-Virgo events under consideration here, the best fit value for $\sigma$ is $0.27^{+0.14}_{-0.15}$, which, remarkably, is inconsistent with zero at 90% confidence level. Notice that the distribution does not change much compared to the general case where we marginalize over $\mu$, indicating that the data is largely insensitive to the value of $\mu$.

2.3 Benchmark spin models for astrophysical BH

We are interested in comparing the $\chi_{\rm eff}$ distribution for PBH, discussed above, with what predicted in the case of astrophysical black holes. Given the current status of observations and the output of population synthesis codes, it is presently unwarranted to try to reproduce specific binary black hole spin distribution reflective of given astrophysical formation processes. Rather, it has become somewhat customary in the literature to adopt simplified benchmark models for the alignment and intrinsic spin distributions of astrophysical black holes, following what proposed in ref. [25], and endorsed and utilized by the LIGO collaboration [2] and by others (see e.g. [55]). We shall assume that the merging black holes have equal mass (see a discussion of the effect of unequal mass mergers on $\chi_{\rm eff}$ in the Appendix), and that the distribution for the spin magnitude is statistically independent of that for the spin alignment.

Noting that the spin directions for isolated binary black holes are thought to be dominantly aligned (see e.g. [56, 57, 58]), we choose a distribution for the spin direction with perfect alignment as an extreme case. We note that this assumptions reflects under any circumstances an extreme, idealized situation for a variety of reasons: for instance, there exists evidence for spin-orbit misalignment in black hole X-ray binaries [59], and effects from the supernova explosion could also contribute to tilt the spin-orbit angle (natal kicks) [60]. We study the systematic effects of relaxing the assumption of perfect alignment in the Appendix, see in particular Fig. 10, bottom panel.

We parameterize the astrophysical spin magnitude following the spin distribution proposed by [25] and used by LIGO [2] and [55]. The models consist of $3$ different spin magnitude distributions:

• •

a low (intrinsic) spin distribution $p(\chi)=2(1-\chi)$ (L),

• •

a flat spin distribution $p(\chi)=1$ (F), and

• •

a high spin distribution $p(\chi)=2\chi$ (H).

We reproduce the distributions in the left panel of Fig. 3, together with the PBH intrinsic spin distribution for PBH for the central value of $\sigma=0.27$ and $\mu=0$ inferred above.

Also following ref. [25], we consider two spin-orbit distribution orientations: aligned and isotropic. Notice that the tilt angle is an excellent tracer of BBH formation channels [61], with the aligned distribution expected for isolated binary formation channel [44, 27, 37, 62], under the simplifying assumption that the binaries remain perfectly aligned throughout their evolution (an assumption that could be violated by effects such as supernova natal kicks, although mass transfer and tidal interactions might work in the opposite direction and tend to re-align the binary). The isotropic distribution is motivated by dynamical formation mechanisms in dense stellar environments or similarly disordered assembly scenarios [63, 64, 65], as well as by what expected for PBH [66].

The final prior distribution for $\chi_{\rm eff}$ for the various models under consideration is shown in the right panel of Fig. 3, where we do not include the HA model which is already strongly excluded by data. Notice that since the intrinsic spin distribution is positive-definite, “aligned” models do not allow for negative values of $\chi_{\rm eff}$; finally, also notice how the prior distribution for “isotropic” models is symmetric in $\chi_{\rm eff}$.

3 Analysis and Results

We present here our results for the odds ratios of the different prior distributions for $\chi_{\rm eff}$ outlined in the previous section, as well as the posterior probability density functions for a “mixed” scenario with PBH providing a fraction $f$ of the BBH events. We then discuss how, under different assumptions, such odds ratios will evolve with additional events in the future, and how knowledge of which fraction of the events originates from which prior distribution will change with greater statistics (sec. 3.2).

3.1 Comparing models to observations: odds ratios and mixture

We confront here the prior distributions obtained in the previous section with data by calculating odds ratios, which quantify the statistical support for a model over another, allowing us to compare models and giving us a statistically-motivated selection criterion. Fig. 4 shows the odds ratios between all possible models and the reference low-intrinsic-spin, isotropic (LI) model; what we show is therefore defined as $p(d|\mathrm{M})/p(d|\mathrm{LI})$ for model $\mathrm{M}$ given the 10 events $d$. Here, we have already assigned equal probability to the prior probability distributions of each and every model (i.e., we assume that all models are equally likely a priori). This implies that the odds ratio and the Bayes factor are equivalent.

Fig. 4 shows how all models with aligned spin distribution are significantly disfavored with respect to the isotropic ones. Also if we compare the “favored” aligned model, which is the one corresponding to a low intrinsic spin distribution (LA), with the least favored isotropic model, which is the one with high intrinsic spin distribution (HI), the odds ratio in favour of HI is still very large, at $36:1$.

In addition to showing a strong statistical evidence for isotropy of spin and orbit orientations, the data favour models with small intrinsic spin magnitude distributions and heavily disfavour those with high spin. The two best models are the PBH and LI with the PBH model slightly preferred over the LI with an odds ratio of $3:2$.

Notice that the intrinsic spin distribution we assume for PBH was calculated in ref. [29] by integrating the probability density $P(\chi,\delta)$, with $\chi$ the intrinsic spin and $\delta$ the overdensity giving rise to the PBH formation, over $\delta$. Critically, ref. [29] assumes no correlation between $\delta$ and $\chi$. When this assumption is relaxed, one generally finds much narrower intrinsic spin distributions [30]. In the notation of ref. [30], the probability distribution for $\chi$ in the case of PBH depends on the parameter $\nu$, defined as the ratio of the critical collapse overdensity and the variance of the overdensity at horizon crossing, and on the parameter $\gamma$ which effectively measures the width of the PBH mass function, with $\gamma=1$ for a monochromatic power spectrum.

Figure 5 of ref. [30] shows that the relevant range for the parameter $\nu$ for LIGO-sized PBH is around $\nu\sim 6$. While the range for the parameter $\gamma$ depends on the PBH mass function, mass functions peaked around a few solar masses typically have values of $\gamma\sim 0.85...0.88$ (see their sections 7.2 and 7.3), although a broader range is possible. Here, we take as benchmark cases $(\nu=6,\ \gamma=0.8)$ (which we indicate in Fig. 4 as PBH 2) and $(\nu=6,\ \gamma=0.88)$ (PBH 3 in Fig. 4), the latter leading to the narrowest possible prior distribution for $\chi_{\rm eff}$, and the former with a broader distribution. We show the prior distribution for the two models in the bottom panel of Fig. 9.

We find odds ratio of $0.44$ and $0.37$ respectively for PBH 2 and 3, indicating (since both odds ratios are smaller than 1) a statistical preference for the LI model as well as for our benchmark PBH model over these narrower intrinsic PBH spin distribution. As expected, our current ability to distinguish between models with very low spin distribution is very limited with the available data. Similar results and conclusions have been found in [22, 55]. For reference and to summarize our findings, we list in Table 2 the natural log Bayes factors for various models compared to the benchmark LI model.

The actual LIGO BBH population likely reflects a mixed population of two different models (or more). Here, we are going to assume that the mixture is our PBH model with the second one any of the following models: LI, LA, FI, FA and HI. In Fig. 5 we show the posterior probability density for the fraction of the BBH mergers coming from PBH, which we indicate with $f$, where $f=1$ means that all the events are from PBH and $f=0$ means the opposite, i.e. all the events are coming from the second model and none from PBH. Notice that if we allow a mixed model we find a statistical preference for the majority of events coming from a PBH population. For the case where the mixed model consists of PBH and FA (dashed orange line), the $f$ distribution peaks around $\sim 0.9$, therefore favoring a scenario where 9 events come from a PBH population and 1 event from the FA population. This is somewhat expected because there clearly is one event (GW170729) that could have come more likely from a population favoring large $\chi_{\rm eff}$ values such as what predicted in the FA prior distribution rather than PBH. A similar conclusion can be drawn when the second mixture model is LA (dashed green line). In this second case, there are two events that could be ascribed to a LA distributions: GW170729 and GW170729. This is why the $f$ posterior distribution peaks around $\sim 0.85$, a little lower than the case of FA. For the case of a mixture of FI (solid orange line) or HI (solid pink line) with PBH, the probability distribution for $f$ is flatter, but we still can conclude that more than half of the events are coming from a PBH population. Lastly, in the case of a mixture model consisting of PBH and LI (solid green line) the distribution is almost half PBH and half LI with a slight preference for the PBH model, as expected from the odds ratio between the two models, which have comparably similar $\chi_{\rm eff}$ prior distributions.

3.2 Future events

To test the sensitivity of our setup to the different benchmark models under consideration, we simulate future events for each of the six population models under consideration: PBH, LI, LA, FI, FA and HI. We generate mock observations following the same approach as in refs. [25, 22]. First, We approximate each LIGO event as a Gaussian with the same mean value and $90\%$ credible interval; Second, we draw a value of $\chi_{\rm eff}^{\mathrm{true}}$ from the population’s distribution we want to simulate; third, we generate an observation from the distribution $\chi_{\rm eff}^{\mathrm{obs}}\sim\mathcal{N}(\chi_{\rm eff}^{\mathrm{true}},\sigma^{\mathrm{unc}})$, where $\sigma^{\mathrm{unc}}$ is a random uncertainty from one of LIGO’s ten events. Finally, the we calculate the posterior probability as $\sim\mathcal{N}(\chi_{\rm eff}^{\mathrm{true}},\sigma^{\mathrm{unc}})p_{\mathrm{FI}}(\chi_{\rm eff})$.

Assuming that all events are coming from the same population, we simulate 200 events for three possible “true” scenarios: FI, LI or PBH. Fig. 6 shows the dependence of the odds ratio with respect to the number of extra events from LI (left) and FI (right) populations. Notice that the odds ratios for the models already start at different values, because the current 10 LIGO events are included: the starting point for each model is thus just the odds ratio from Table 2.

For the fully LI-simulated population, our results show that with only 10 extra events the FA and LA models might be disfavoured at more than the $5\sigma$ level, and that 75 extra events are needed to reject the the HI model at the same confidence level. The LI and PBH models would be disfavored at the same level with even more events, reaching close to a $5\sigma$ level with 200 extra events.

Under the assumption that the true population is FI, the evolution of the odds ratios is shown in the right panel of Fig. 6. With less than 10 extra events we find that it would be possible to heavily disfavor both the FA and LA models; interestingly, the PBH model can also be rejected at the $5\sigma$ level with only 50 extra events. The entire 200 extra events would allow for a $5\sigma$ rejection threshold for the LI model and $3\sigma$ for the HI model. Finally, if the true population is that of PBHs, as in Fig. 7, all models except LI can be rejected at $5\sigma$ level with only 75 extra events. We find that 200 extra events would be needed to discriminate the LI model over the PBH one at the $4\sigma$ level.

Once again, we emphasize that there is no reason to assume all events observed by LIGO are coming from the same single population. In what follows we therefore consider a mixed population consisting of half PBH and half of a second model among LI, FI, HI or LA. We simulate up to 500 mixed events for each of these mixed populations.

In the upper panels of Fig. 8, we show that it is possible to infer the relative fraction of PBH $f$ (which is given on the horizontal axis) with high confidence with 200 events for the case of PBH-LA and PBH-HI mixture models. To be able to discriminate the mixture fraction between PBH and FI, more events are needed than in the previous two cases, as shown in the bottom left panel of Fig. 8. After detecting 500 events the value of $f$ can be determined fairly precisely. Additional events are needed when the two prior distributions for the two models are similar, such as for the case of a mixed PBH-LI population, bottom right in Fig. 8. This notwithstanding, the value of $f$ peaks at the correct value of $\sim 0.5$ but the spread is still substantial, even with 500 extra events.

Run O3 of LIGO-Virgo has already commenced, and the projected inferred rate of BBH mergers from the previous runs is around $9.7-101\,\rm{Gpc}^{-3}\,\rm{y}^{-1}$ [33]. At this moment, the number of putative candidate run O3 BBH merger events is 10, from a period of approximately 1.5 months222https://gracedb.ligo.org/latest/, which would imply an approximate total number of events per year of around 80.

4 Discussion and Conclusions

In this study we considered how measurements of the effective spin parameter $\chi_{\rm eff}$ provide information on the origin of merging binary black holes observed with gravitational wave telescopes. To this end, we presented a calculation of the prior distribution for primordial black holes as well as for a few representative benchmark distributions possibly indicative of what expected in simple astrophysical black hole binary populations.

In the case of PBHs, following ref. [29] we assumed no correlation between the overdensity leading to the black hole formation in the early universe and the intrinsic black hole spin; the resulting intrinsic spin distribution is a positive-definite half-Gaussian with zero mean; to fix the width of the intrinsic spin distribution we calculated the prior distribution for the effective spin parameter $\chi_{\rm eff}$ and we utilized 10 LIGO-Virgo measurements of $\chi_{\rm eff}$ to determine a best-fit intrinsic spin width.

We then proceeded to compare odds ratios for the current set of 10 measurements for $\chi_{\rm eff}$ for the various models under consideration. We showed that with current data the $\chi_{\rm eff}$ measurements have a marginal, and not highly statistical significant preference for a dominant population of primordial black holes over the best fitting astrophysical model. We also calculated the posterior probability for the relative fraction of primordial versus astrophysical BBH events, finding that there is a preference for a scenario where one or two events originate from a population with preferentially aligned spin-orbit distributions, and the remaining eight-to-nine events from a population with an isotropic spin-orbit distribution and low intrinsic spin.

We then assessed the number of future events needed to disentangle, at a given significance level, BBH from different populations, assuming that all binaries originate from the same “true” distribution. Generally, non-isotropic spin distributions are highly disfavored even when only considering the current 10 events. In addition, even for isotropic alignment distributions, it will quickly become possible to distinguish models with large intrinsic spin from those with low intrinsic spin magnitude distributions. Assuming that the entirety of the merging black holes have a primordial origin, we anticipate that distinguishing at 3$\sigma$ their $\chi_{\rm eff}$ distribution from a low intrinsic spin, aligned spin distribution will require on the order of 100 additional events.

Finally, we illustrated the number of events necessary to acquire information on the relative fraction of primordial versus non primordial binary black holes. Once again, if the population of non-primordial black holes has a preferentially aligned spin-orbit distribution, such fraction can be pinpointed with relatively few additional events; the posterior distribution for the relative fraction of events from populations with isotropic spin-orbit distribution shows a peak generally pointing to a systematically larger-than-true fraction of primordial black holes, but eventually converging with large-enough statistics to the “true” value.

In the Appendix, we discuss the systematics associated with three key assumptions used in our analysis: (1) the width $\sigma$ of the primordial black hole intrinsic spin distribution, (2) the mass ratio we use to calculate the $\chi_{\rm eff}$ prior distribution, and (3) the assumption of perfect alignment between spin and orbit.

Clearly, as an increasingly large statistics of BBH mergers becomes available, the question of the origin of the population of merging black holes will come into sharper focus not only with studies of the effective spin parameter, but using all other pieces of information, including but not limited to the mass distribution of events, the correlation between mass and spin, localization information, etc. Interesting questions relating for instance to the degree to which BBH accrete will also be tackled, including the issue of whether there could be substantial accretion for PBH and how much accretion would spin up individual PBHs [67]: possibly, if PBHs do significantly accrete, the heavier ones would have both higher spin and higher $\chi_{\rm eff}$. Interestingly, this seems to be the case for the most massive LIGO event (GW170729) and for the recent, claimed detection events GW151216 [68] and GW170403 [69]. The central issue of observational bias in inferring the origin of merging BBH will also be helped greatly by benefiting from increased statistics. Finally, we also expect that theoretical and observational progress will lead to more realistic models for the expected distributions for the intrinsic spin and spin-orbit correlation of different populations of astrophysical black holes, allowing firmer statements on the origin of merging black holes than what is possible with the benchmark models currently in use.

Acknowledgments

We are grateful to Yonatan Kahn and to Benjamin V. Lehmann, John Tamanas and other members of the UCSC particle theory group for helpful comments. This work is partly supported by the U.S. Department of Energy grant number de-sc0010107.

Appendix A Priors

In this Appendix we derive the probability distribution for $\chi_{\rm eff}$ for general given spin magnitude and angular distributions. Subsequently, we show how the distribution for $\chi_{\rm eff}$ differs when we relax some assumptions that we employed in our calculations, specifically the assumption that the merging black holes have the same mass, and the assumption of perfect alignment between the individual spin and the orbital angular momentum for the “aligned” models.

Let us define the direction of the orbital angular momentum as the direction of the $z$-axis, then $\chi_{\rm eff}$ can be re-written as

[TABLE]

where $q$ is the mass ratio $m_{2}/m_{1}$ such that $0\leq q\leq 1$, and $\chi^{z}_{i}=\chi_{i}\cos\theta_{i}$ is the individual spin component along the $z$-axis. Our goal is to derive the distribution for $\chi^{z}_{i}$ and the probability distribution for the sum $\chi^{z}_{1}+q\chi^{z}_{2}$.

First, the distribution of the product for $\chi^{z}_{i}$ is given by the following integral over probability distributions

[TABLE]

where the probability density functions $f_{\chi_{i}}$ is just the spin magnitude distribution, as given, for example, for PBH in Eq. (2.4), and where $f_{\cos\theta_{i}}$ is the distribution for the cosine of $\theta_{i}$ ($\theta_{i}$ is defined under Eq. (2.2)). For instance, if the angular distribution is isotropic the $p(\chi^{z})$ distribution reads as

[TABLE]

Using the probability distribution for $\chi^{z}_{i}$, and assuming that the individual spins are independent of each other, we find that the probability distribution for $\chi_{\rm eff}$ is given by the following convolution of the distributions

[TABLE]

where the extrema of integration read

[TABLE]

In general, there is no analytic solution for the probability density, and we have to perform the integration numerically to solve for the distribution.

A.1 PBH

In figure 9, we compare the LIGO prior distribution with the results for the prior distribution for $\chi_{\rm eff}$ for our PBH models, for various values of $\sigma=0.1,$ 0.2 and 0.3 (top left) and of the mass ratio $q$ (top right). By allowing the width to vary we get a sense of how sensitive is the PBH distribution with respect to $\sigma$ is.

The top right panel, instead, illustrates how the choice of $q$ has a very marginal impact on the functional shape of the predicted prior distribution. Notice that all events reported by LIGO have a median value of the mass ratio q larger than $0.5$ and larger than $0.2$ at $90\%$ credible interval (figure 5 on ref. [33]), which is why we focus on relatively large mass ratios, $q=0.2$, 0.5 and 1.0.

Finally, the bottom panel shows our benchmark $q=1$ and $\sigma=0.27$ PBH prior distribution with what is predicted in the models of ref. [30]. Notice that as illustrated in Fig. 5 of ref. [30], the relevant value for $\nu\simeq 6$ for masses in the 10 $M_{\odot}$ range, while the choice of $\gamma$ reflects the prediction for a nearly flat power spectrum ($\gamma=0.88$) and a slightly smaller value, which could result from e.g. a broader power spectrum (see e.g. their Fig. 6, left panel).

A.2 Benchmark spin models

Here we study the systematic dependence of the prior distribution for $\chi_{\rm eff}$ on the values of the mass ratio $q$ (Fig. 10, top panels) and on the assumption of perfect alignment (bottom panel).

The top left panel shows how the LI, FI and HI prior distributions change when switching to our benchmark value of $q=1$ (equal mass ratio) to $q=0.5$. The top right panel does the same for the aligned models LA and FA. It is clear that the changes for the isotropic models is very small, while for the aligned models there is a slight but noticeable effect near the peak of the distribution for the FA model, and a slight shift in the peak location for the LA model. For the aligned models, as the mass ratio starts to get more extreme, the distribution for $\chi_{\rm eff}$ begins to resemble the distribution for the spin magnitude.

The bottom panel of Fig. 10 shows how the prior distributions for the LA and FA models change when the assumption of perfect BH spin-orbit is relaxed. If the BBH is formed via classical isolated binary evolution with an initial perfect alignment, large values for the misalignment angle are required to be consistent with the GW151226 event [71], but the natal kicks necessary to explain this large misalignments usually exceed the typical values for binary evolution models. Therefore, the distribution of the tilt angle is taken to be flat from [math] up to $35^{\circ}$, and thus the maximal angle is within the region suggested by ref. [71]. The figure shows that the prior distribution is hardly affected at all, with the only noticeable change being a shift of the peak distribution to slightly lower values of $\chi_{\rm eff}$.

Finally, we consider the scenario of binary BH evolution from isolated massive binary systems studied in ref. [70]; in this scenario $\chi_{\rm eff}$ can take negative values because of effects such as initially misaligned binary components spins; while other effects such as mass exchange, tidal interaction, and common envelope evolution tend to align the spins, the resulting distribution for $\chi_{\rm eff}$ was found to range between -0.2 and 0,8, thus encompassing negative values (see their Figure 11). To mimic this physical situation, here we assumed (following what shown in fig. 8 of ref. [70]) that one black hole has a spin orientation that is isotropic compared to the orbital angular momentum, and has low intrinsic spin magnitude distribution; for the second black hole, we assumed perfect alignment, and a flat intrinsic spin magnitude distribution. The resulting $\chi_{\rm eff}$ distribution is shown in Fig. 10 with a dot-dashed light-blue line; our results reflect qualitatively the range for $\chi_{\rm eff}$ found in ref. [70]. We calculated the resulting odds ratio for this setup, and found that compared to the benchmark LI model, the odd ratio is of -6.04, thus comparable to what we found for the LA model.