Are merging black holes born from stellar collapse or previous mergers?

TL;DR

This paper investigates whether gravitational-wave data can distinguish black hole mergers originating from stellar collapse versus those from previous mergers, using Bayesian analysis of population models and observed parameters.

Contribution

It introduces a Bayesian model selection method to differentiate black hole formation scenarios based on gravitational-wave observations of mass, spin, and redshift.

Findings

Bayesian analysis can effectively distinguish formation scenarios

Mass and spin distributions differ between first and second-generation black holes

Observations can identify the formation history of merging black holes

Abstract

Advanced LIGO detectors at Hanford and Livingston made two confirmed and one marginal detection of binary black holes during their first observing run. The first event, GW150914, was from the merger of two black holes much heavier that those whose masses have been estimated so far, indicating a formation scenario that might differ from "ordinary" stellar evolution. One possibility is that these heavy black holes resulted from a previous merger. When the progenitors of a black hole binary merger result from previous mergers, they should (on average) merge later, be more massive, and have spin magnitudes clustered around a dimensionless spin ~0.7. Here we ask the following question: can gravitational-wave observations determine whether merging black holes were born from the collapse of massive stars ("first generation"), rather than being the end product of earlier mergers ("second generation")? We construct simple, observationally motivated populations of black hole binaries, and we use Bayesian model selection to show that measurements of the masses, luminosity distance (or redshift), and "effective spin" of black hole binaries can indeed distinguish between these different formation scenarios.

Tables

Table 1. Table 1: Odds ratios from the three O1 observations (GW150914, GW151226 and LVT151012) and from hypothetical observations of the same events at Advanced LIGO design sensitivity. Odds ratios in parentheses were computed omitting all redshift information, i.e. considering the 3-dimensional vector of observables 𝐮 = { M , q , χ eff } 𝐮 𝑀 𝑞 subscript 𝜒 eff {\bf u}=\{M,\,q,\,\chi_{\rm eff}\} .

		$1g+1g$ vs. $2g+2g$	$1g+1g$ vs. $1g+2g$	$1g+2g$ vs. $2g+2g$
O1 LIGO	flat	12.7 (15.8)	2.0 (2.0)	6.4 (7.6)
	log	3.3 (3.5)	0.9 (0.9)	3.5 (3.8)
	power law	0.7 (1.0)	1.3 (1.6)	0.6 (0.6)
Ad. LIGO (design)	flat	30.2 (37.8)	1.4 (3.7)	21.9 (10.11)
	log	4.3 (7.0)	0.6 (1.4)	6.9 (5.1)
	power law	0.6 (1.7)	1.0 (3.8)	0.6 (0.5)

Table 2. Table 2: Number of observations needed to distinguish populations at 5 ​ σ 5 𝜎 5\sigma with 5 % percent 5 5\% , 50 % percent 50 50\% and 95 % percent 95 95\% probability. The “true” model is marked by a T in the column header. For instance, in column 1g1gT /2g2g we compare models 1g+1g and 2g+2g when observations are drawn from the 1g+1g catalog.

$N_{\mathrm{obs}}$ at $5\sigma$		1g1gT/2g2g			2g2gT/1g1g			1g1gT/1g2g			1g2gT/1g1g			1g2gT/2g2g			2g2gT/1g2g
		$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$	$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$	$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$	$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$	$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$	$\mathrm{\hspace{0.33em}5}\%$	$\mathrm{\hspace{0.17em}50}\%$	$\mathrm{\hspace{0.17em}95}\%$
LIGO O1	flat	27	53	100	31	57	103	40	76	143	44	80	146	50	105	204	77	133	233
	log	27	52	94	25	50	94	30	58	106	29	56	106	42	86	165	59	104	182
	power law	14	29	57	19	35	64	7	17	34	13	23	41	31	61	114	35	64	117
Ad. LIGO	flat	23	46	86	26	50	91	45	82	146	37	73	139	37	83	170	73	122	206
	log	20	41	79	24	45	83	41	73	132	33	66	122	26	56	112	48	81	138
	power law	20	39	72	18	37	70	10	21	40	11	22	41	15	31	61	20	37	67
A+	flat	18	39	75	22	43	79	46	83	149	34	69	136	34	80	165	75	123	211
	log	16	34	65	19	38	69	41	73	131	30	62	120	22	51	107	50	81	136
	power law	17	35	67	17	34	65	10	22	41	10	21	40	12	27	52	20	35	61
Voyager	flat	6	15	33	10	21	39	34	69	128	27	62	122	13	36	80	36	61	102
	log	4	11	25	8	17	32	25	53	102	20	51	101	8	23	51	26	44	73
	power law	5	13	26	7	16	31	9	19	37	7	18	36	4	11	24	12	21	35

Table 3. Table 3: Model comparison tests between populations characterized by different merger generations and mass distributions using the Advanced LIGO sensitivity curve. For each injected distribution and model comparison we report the preferred population in the limit where N obs → ∞ → subscript 𝑁 obs N_{\rm obs}\to\infty (in practice we use N obs = 10 3 subscript 𝑁 obs superscript 10 3 N_{\rm obs}=10^{3} ). The true population (T) is correctly identified whenever it is among those tested. While most of the comparisons correctly identify the merger generation (rows denoted by a check mark ✓), in some cases making the wrong assumption on the underlying mass distribution prevents a correct identification (rows denoted by a cross ✗). In one case (second row) we obtained odds ratios consistent with one even when N obs → ∞ → subscript 𝑁 obs N_{\rm obs}\to\infty , so that no conclusions can be drawn and the comparison is marked as “not significant.” In all other cases the behavior of 𝒪 ​ ( N obs ) 𝒪 subscript 𝑁 obs \mathcal{O}(N_{\rm obs}) is qualitatively similar to Fig.  6 , i.e. the odds ratio grows exponentially until populations can be distinguished at 5 ​ σ 5 𝜎 5\sigma .

Injection		Test		Preferred
flat 1g+1g	flat 1g+1g	vs	flat 2g+2g	flat 1g+1g	✓T
	log 1g+1g	vs	log 2g+2g	not significant
	power law 1g+1g	vs	power law 2g+2g	power law 1g+1g	✓
flat 2g+2g	flat 1g+1g	vs	flat 2g+2g	flat 2g+2g	✓T
	log 1g+1g	vs	log 2g+2g	log 2g+2g	✓
	power law 1g+1g	vs	power law 2g+2g	power law 1g+1g	✗
log 1g+1g	flat 1g+1g	vs	flat 2g+2g	flat 1g+1g	✓
	log 1g+1g	vs	log 2g+2g	log 1g+1g	✓T
	power law 1g+1g	vs	power law 2g+2g	power law 1g+1g	✓
log 2g+2g	flat 1g+1g	vs	flat 2g+2g	flat 1g+1g	✗
	log 1g+1g	vs	log 2g+2g	log 2g+2g	✓
	power law 1g+1g	vs	power law 2g+2g	power law 1g+1g	✗
power law 1g+1g	flat 1g+1g	vs	flat 2g+2g	flat 1g+1g	✓
	log 1g+1g	vs	log 2g+2g	log 1g+1g	✓
	power law 1g+1g	vs	power law 2g+2g	power law 1g+1g	✓T
power law 2g+2g	flat 1g+1g	vs	flat 2g+2g	flat 1g+1g	✗
	log 1g+1g	vs	log 2g+2g	log 1g+1g	✗
	power law 1g+1g	vs	power law 2g+2g	power law 2g+2g	✓T

Equations

$$\frac{da}{dt}=-\frac{64}{5}\frac{q}{(1+q{)}^2}\frac{M^3}{a^3}\frac{G^3}{c^5},$$

$$t={\int }_a^0\frac{dt}{d{a}^{\prime }}d{a}^{\prime }=\frac{5}{256}\frac{(1+q{)}^2}{q}\frac{a^4}{M^3}\frac{c^5}{G^3}.$$

$$\frac{dn}{dt}=\frac{dn}{da}\frac{da}{dt}\propto \frac{1}{a^4}\propto \frac{1}{t}.$$

$$t_L=\frac{1}{H_0}{\int }_0^z\frac{dz}{(1+z)\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$$t_L(\tilde{z})-t_L(z)=t_D.$$

$${\chi }_{\mathrm{eff}}=\frac{1}{M}\left(\frac{{\mathbf{S}}_1}{m_1}+\frac{{\mathbf{S}}_2}{m_2}\right)\cdot \hat{\mathbf{L}}.$$

$$\mathbf{u}=\{M,q,z,{\chi }_{\mathrm{eff}}\}.$$

$${\rho }^2=4{\int }_0^{\infty }\frac{|\tilde{h}(f)|}{S_n(f)}df,$$

$${\kappa }^{(i)}=P\left(w^{(i)}\right),$$

$${\sigma }_{\mathrm{Detector}}={\sigma }_{\mathrm{LIGO}}\frac{{\rho }_{\mathrm{LIGO}}}{{\rho }_{\mathrm{Detector}}}.$$

$$D_L=\frac{1+z}{H_0}{\int }_0^z\frac{dz}{\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$$\frac{d{D}_L}{dz}=\frac{D_L}{1+z}+\frac{1+z}{H_0\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}}.$$

$${\left(\frac{{\sigma }_z}{D_L}\frac{d{D}_L}{dz}\right)}^2={\left(\frac{{\sigma }_{D_L}}{D_L}\right)}^2+{\left(\frac{{\sigma }_{H_0}}{H_0}\right)}^2$$

$$\tilde{r}(\mathbf{u},\lambda )=\frac{{\sum }_{i=1}^I{\kappa }^{(i)}{\prod }_{j=1}^J\delta \left(u_j-{\bar{u}}_j^{(i)}\right)}{{\sum }_{i=1}^I{\kappa }^{(i)}},$$

$$\tilde{r}(\mathbf{u},\lambda )=\frac{{\sum }_{i=1}^I{\kappa }^{(i)}{\prod }_{j=1}^J\mathcal{N}\left(u_j;{\bar{u}}_j^{(i)},{\sigma }^{(i)}\right)}{{\sum }_{i=1}^I{\kappa }^{(i)}},$$

$$\mathcal{N}\left(u_j;{\bar{u}}_j^{(i)},{\sigma }^{(i)}\right)=\frac{1}{{\sigma }^{(i)}\sqrt{2\pi }}\exp \left(-\frac{u_j-{\bar{u}}_j^{(i)}}{2{\sigma }^{2(i)}}\right).$$

$${\tilde{r}}_{k_1,...,k_J}(\lambda )={\int }_{b_{k_1}}^{B_{k_1}}d{u}_1...{\int }_{b_{k_J}}^{B_{k_J}}d{u}_J\tilde{r}(\mathbf{u},\lambda )$$

$$=\frac{{\sum }_{i=1}^I{\kappa }^{(i)}{\prod }_{j=1}^J{\int }_{b_{k_j}}^{B_{k_j}}\mathcal{N}\left(u_j;{\bar{u}}_j^{(i)},{\sigma }^{(i)}\right)d{u}_j}{{\sum }_{i=1}^I{\kappa }^{(i)}}.$$

$$\sum _{k_1}...\sum _{k_J}{\tilde{r}}_{k_1,...,k_J}=1.$$

$$r_k(\lambda )=N(\lambda ){\tilde{r}}_k(\lambda ).$$

$$p(\mathbf{d}|\lambda )=\prod _k\frac{({r_k(\lambda ))}^{d_k}{e}^{-r_k(\lambda )}}{d_k!}.$$

$$p(\mathbf{d}|\lambda )=\left(\prod _k\frac{({{\tilde{r}}_k(\lambda ))}^{d_k}{e}^{-{\tilde{r}}_k(\lambda )}}{d_k!}\right)\left(N(\lambda {)}^{{\sum }_k{d}_k}{e}^{-N(\lambda )}\right),$$

$$\tilde{p}(\mathbf{d}|\lambda )=\left(\prod _k\frac{({{\tilde{r}}_k(\lambda ))}^{d_k}{e}^{-{\tilde{r}}_k(\lambda )}}{d_k!}\right)\sum _N\left(N^{{\sum }_k{d}_k}{e}^{-N}\right).$$

$$O_{AB}=\frac{p(\mathbf{d}|A)\pi (A)}{p(\mathbf{d}|B)\pi (B)},$$

$$p_A=\frac{O_{AB}}{1+O_{AB}}=\frac{p(\mathbf{d}|A)}{p(\mathbf{d}|A)+p(\mathbf{d}|B)},$$

$$\mathcal{O}=\frac{1}{1-2\mathrm{erf}(\sigma )}.$$

$$p(\lambda |\mathbf{d})=\frac{p(\mathbf{d}|\lambda )\pi (\lambda )}{\int p(\mathbf{d}|\lambda )\pi (\lambda )d\lambda }.$$

$$p(d_k)=\frac{{r_k({\lambda }_{\mathrm{true}})}^{d_k}{e}^{-r_k({\lambda }_{\mathrm{true}})}}{d_k!}.$$

$$\displaystyle\lambda_{\rm true}\equiv\big{\{}f_{\rm 1g+1g},f_{\rm 1g+2g},f_{\rm 2g+2g}\big{\}}=\big{\{}0.6,0.1,0.3\big{\}}\,.$$

Full text

I Introduction

The observation of gravitational waves (GWs) from merging black hole (BH) binaries was a milestone in physics and astronomy Abbott et al. (2016a) (LIGO Scientific Collaboration and Virgo Collaboration); Abbott et al. (2016b) (LIGO Scientific Collaboration and Virgo Collaboration); Abbott et al. (2016c) (LIGO Scientific Collaboration and Virgo Collaboration). During their first observing run (O1), the Advanced LIGO detectors detected two GW events (GW150914 and GW151226) and a marginal candidate LVT151012, which is also likely to be of astrophysical origin. The second observing run (O2) is currently ongoing, and Advanced Virgo is expected to join the detector network soon. Dozens of BH mergers may be detected by the end of O2 or in the third run (O3), allowing for statistical studies of their populations.

These events can further our understanding of the formation channels of binary BHs (Abbott et al. (2016d), LIGO Scientific Collaboration and Virgo Collaboration), because different astrophysical scenarios predict different binary properties. As the number of detections grows, a statistical analysis of the observed binary parameters should eventually allow us to identify or constrain the physical processes responsible for the formation and merger of compact binaries. Currently favored scenarios include stellar evolution of field binaries (Postnov and Yungelson, 2014) and the dynamical capture of BHs in globular clusters (Benacquista and Downing, 2013). Recent work showed that both field formation (Mapelli et al., 2009; Belczynski et al., 2010; Dominik et al., 2012, 2013; Spera et al., 2015; Dominik et al., 2015; Belczynski et al., 2016a, b) and cluster formation (Rodriguez et al., 2015, 2016a; Chatterjee et al., 2017; Rodriguez et al., 2016b) are broadly compatible with current Advanced LIGO observations (Abbott et al. (2016d), LIGO Scientific Collaboration and Virgo Collaboration).

It is quite likely that both field and cluster formation channels are at work in nature. The first event, GW150914, was the most surprising, because the merging BHs are much heavier that those whose masses have been estimated so far in x-ray binaries Remillard and McClintock (2006); Narayan and McClintock (2013), indicating a formation scenario that might differ from “ordinary” stellar evolution. Alternative theoretical scenarios which could explain the unexpected properties of GW150914 include formation via hierarchical triples (Wen, 2003; Antonini et al., 2014, 2016), a Population III origin for the binary members (Kinugawa et al., 2014; Hartwig et al., 2016; Belczynski et al., 2016c), chemically homogeneous evolution in short-period binaries (Mandel and de Mink, 2016; de Mink and Mandel, 2016; Marchant et al., 2016), and a primordial origin for the merging BHs (Bird et al., 2016; Cholis et al., 2016).

One possibility to explain the high mass of the merging BHs in GW150914 is that these BHs did not form following stellar collapse, but rather from previous BH mergers. Field formation scenarios typically predict long delay times between the formation and merger of a BH binary Dominik et al. (2012), so repeated mergers seem unlikely. However gravitational encounters are more common in dense stellar environments, and some scenarios suggest that repeated mergers may be possible Antonini and Rasio (2016); Mapelli (2016); O’Leary et al. (2016). The most likely environment to host multiple mergers are nuclear clusters Antonini and Rasio (2016), which present larger escape speeds compared to globular and open clusters and can, therefore, more easily retains merger remnants with substantial recoils Merritt et al. (2004). Stellar-mass BH binaries may also form in AGN gaseous discs McKernan et al. (2017), where migration traps can be invoked to assemble multiple generations of mergers Bellovary et al. (2016). Primordial BHs are also expected to merge very quickly Bird et al. (2016); Cholis et al. (2016), so the possibility of repeated mergers in this scenario should not be excluded Clesse and García-Bellido (2017).

In this paper, we ask the following question: can GW observations determine whether merging BHs such as those in GW150914 were born directly from the collapse of massive stars (“first-generation” BHs, henceforth 1g) rather than being the end product of previous mergers (“second-generation” BHs, henceforth 2g)?

Roughly speaking, one can expect mergers to leave several statistically observable imprints in 2g BHs, namely,

1. (i)

2g BHs should be more massive than BHs born from stellar collapse; 2. (ii)

quite independently of the distribution of spin magnitudes following core collapse (which is highly uncertain Kushnir et al. (2016)), the spin magnitudes of 2g BHs should cluster (on average) around the dimensionless spin $\sim 0.7$ resulting from the merger of nonspinning BHs Berti and Volonteri (2008); 3. (iii)

statistically, the merger of BH binaries including 2g components should occur later (i.e., at smaller redshift or luminosity distance from GW detectors) because of the delay time between BH formation and merger.

In this paper we make these arguments more quantitative and rigorous by developing a simple but physically motivated model to describe the bulk theoretical properties of 1g and 2g binary BH mergers (Sec. II). Then we consider a set of present and future GW detectors, and we simulate observable distributions by selecting detectable binaries and estimating the expected measurement errors on their parameters (Sec. III). Finally we set up a Bayesian model selection framework (Sec. IV) to address what can be done with current observations, and to quantify the capabilities of future detectors to distinguish between different models (Sec. V). We conclude by summarizing our results and pointing out possible extensions (Sec. VI).

II Theoretical distributions

Our goal in this section is to develop a simple prescription to build populations of binary BHs. Our greatly oversimplified model is not meant to capture the complexity of binary evolution in an astrophysical setting, but just the main features distinguishing 1g and 2g BHs.

As illustrated by the cartoon in Fig. 1, we construct three theoretical distributions, labeled by “1g+1g,” “1g+2g” and “2g+2g”. In this context, “1g” means that one of the binary components is a first-generation BH produced by stellar collapse, whereas “2g” means that it is a second-generation BH produced by a previous merger.

II.1 The 1g+1g population

Following the LIGO-Virgo Scientific Collaboration Abbott et al. (2016c) (LIGO Scientific Collaboration and Virgo Collaboration), for the 1g+1g population, we adopt three different prescriptions for the distribution of source-frame masses:

• (i)

Model “flat”: we assume uniformly distributed source-frame masses $m_{1}$ and $m_{2}$ in the range $m_{i}\in[5M_{\odot},50M_{\odot}]$ ($i=1,\,2$), where hereafter $m_{1}>m_{2}$.

• (ii)

Model “log”: we take the logarithm of the source-frame masses to be uniformly distributed in the same range, so that the probability distribution $p(m_{1},m_{2})\propto 1/m_{1}m_{2}$.

• (iii)

Model “power law”: we adopt a power-law distribution with spectral index $\alpha=-2.5$ for the primary BH (i.e. $p(m_{1})\propto m^{\alpha}$), while the secondary mass is uniformly distributed in $m_{2}\in[5M_{\odot},m_{1}]$.

The upper limit of $50M_{\odot}$ was chosen to be consistent with current LIGO compact binary coalescence searches, and it excludes “by construction” intermediate-mass BH searches, discussed e.g. in Abbott et al. (2017) (LIGO Scientific Collaboration and Virgo Collaboration). Moreover, pair instability and pulsation pair instability in massive helium cores Heger and Woosley (2002); Woosley (2017) may inhibit the formation of 1g BHs with masses larger than $\sim 50M_{\odot}$ Belczynski et al. (2016d). If multiple mergers occur through mass segregation in stellar clusters, the more massive objects will tend to form binaries, thus increasing the component masses of 1g+2g and 2g+2g populations. Our $50M_{\odot}$ upper mass limit is therefore conservative, because physical mechanisms such as pair instabilities and mass segregation would further separate the mass distributions of populations involving multiple mergers and make them more easily distinguishable.

Given the great uncertainties on the spin magnitude and orientation of binary BHs Belczynski et al. (2008); Miller and Miller (2015); Belczynski et al. (2016e); Rodriguez et al. (2016c), in all three cases we assume the dimensionless spin magnitudes $\chi_{1,2}$ to be uniformly distributed in $[0,1]$, and their directions to be isotropically distributed.111Rodriguez et al. Rodriguez et al. (2016c) argued that massive field binaries should typically have aligned spins because “heavy” BHs receive small supernova kicks that are unable to tilt the orbit Kalogera (2000); Gerosa et al. (2013), while the spins of massive binaries produced in dense stellar environments should be isotropically distributed. A more detailed investigation of the correlation between spin alignment and binary BH formation requires astrophysical modeling that is beyond the scope of this paper (see e.g. Belczynski et al. (2008); Gerosa et al. (2013)). We are only interested in the global statistical properties of the population. Since isotropic spin distributions stay isotropic under precession and gravitational radiation reaction Bogdanović et al. (2007); Gerosa et al. (2015), the assumption of isotropy will hold also at the small separations relevant for GW observations. For this reason there is no need to carry out post-Newtonian evolutions of the spin distributions for individual binaries of the kind discussed in Kesden et al. (2015); Gerosa et al. (2015); Gerosa and Kesden (2016).

II.2 The 2g+2g population

In order to construct the 2g+2g population we use the following procedure. We randomly extract two binaries from a given 1g+1g population. For these binaries, we estimate the final mass $M_{f}$ and spin $\chi_{\rm f}$ of the merger remnant using the numerical relativity fitting formulas of Refs. Barausse et al. (2012); Barausse and Rezzolla (2009)222There are several alternative fitting formulas for the final masses and spins Lousto and Zlochower (2014); Healy et al. (2014); Zlochower and Lousto (2015); Husa et al. (2016); Hofmann et al. (2016); Healy and Lousto (2017); Jiménez-Forteza et al. (2017). The difference between different prescriptions is smaller than measurement errors in GW observations, and therefore the choice of a particular fitting formula is of no consequence for our present purpose. as implemented in Gerosa and Kesden (2016). These masses and spins are used as input for the second round of binary mergers.

To perform meaningful comparison with the 1g+1g model described above, we again restrict our population to binaries with component masses in the range $[5,50]M_{\odot}$, because this is the mass range targeted by LIGO compact binary coalescence searches.

II.3 The 1g+2g population

The 1g+2g distribution is the obvious mixture of the two: we draw one binary from the 1g+1g distribution, merge it to obtain a 2g BH, and then consider the merger of this 2g BH with a 1g BH.

II.4 Redshift distribution

The redshift distribution of BH mergers in the three different populations should be different, because on average 2g mergers are expected to happen later than 1g mergers. We can estimate the delay times between the formation and merger of a BH binary using the quadrupole formula

[TABLE]

with the result

[TABLE]

If the binary initial separations $a$ are drawn from a log-flat distribution (i.e., $dn/da\propto 1/a$), the distribution of the merger times is also log-flat (cf. Belczynski et al. (2016c)):

[TABLE]

The “lookback time” $t_{L}$ is given by Hogg (1999)

[TABLE]

where we assume $\Omega_{k}=0$, $\Omega_{M}=0.307$, $\Omega_{\Lambda}=0.693$ and $H_{0}=67.7\,{\rm kms}^{-1}{\rm Mpc}^{-1}$ Ade et al. (2016) (Planck Collaboration). From the lookback time we can compute the time $t_{L}(z_{1})-t_{L}(z_{2})$ necessary for the Universe to evolve from redshift $z_{1}$ to redshift $z_{2}$.

We distribute the 1g+1g sources uniformly in comoving volume with redshifts $z<2$. For the 1g+2g population, we assume that 2g BHs formed at some redshift $\tilde{z}$ drawn from the same distribution used for 1g+1g binaries. We then extract a delay time $t_{D}$ from a flat distribution in $\log(t_{D})$ in the range $t_{D}\in[10^{-4}$ Gyrs$,\,t_{L}(\tilde{z})]$. The lower limit is very conservative, and it roughly corresponds to the merger time for a $10M_{\odot}$ BH binary evolving from an initial orbital separation $a=10R_{\odot}$. The redshift $z$ of a 1g+2g merger is then given by the numerical solution of the equation

[TABLE]

Finally, for the 2g+2g population we extract two values $\tilde{z}_{1}$, $\tilde{z}_{2}$ from the 1g+1g distribution. The redshift $z$ of a 2g+2g merger follows again from a numerical solution of Eq. (5), with the difference that now we set $\tilde{z}=\min(\tilde{z}_{1},\tilde{z}_{2})$.

In Sec. V.4 we will discuss how time delay prescriptions affect our results.

II.5 Measurable parameters

For concreteness and simplicity, we will characterize each binary by its total mass $M=m_{1}+m_{2}$, mass ratio $q=m_{2}/m_{1}\leq 1$, redshift $z$ and “effective spin” Abbott et al. (2016e) (LIGO Scientific Collaboration and Virgo Collaboration)

[TABLE]

The effective spin (a mass-weighted sum of the projection of the spins $\mathbf{S_{i}}=m_{i}^{2}\chi_{i}\mathbf{\hat{S}_{i}}$ along the orbital angular momentum $\mathbf{L}$) is a constant of the motion in post-Newtonian evolutions, at least at 2PN order Racine (2008); Gerosa et al. (2015). It is also the easiest spin parameter to measure Pürrer et al. (2016); Abbott et al. (2016e) (LIGO Scientific Collaboration and Virgo Collaboration).

Let us introduce a vector ${\bf u}$ whose components are the observable variables to use in our statistical analysis, i.e.

[TABLE]

The components of this vector will be labeled by an index $j=1,...,J$ such that $u_{1}=M$, $u_{2}=q$, etcetera; a capital Latin index $J$ will denote the dimensionality of the vector $\bf u$, i.e. the number of observables considered in the analysis. Each binary in our catalog is characterized by a specific set of observable properties $\bar{u}^{(i)}$, where the superscript index $(i=1,...,I)$ labels entries in our synthetic catalog.

The theoretical distributions of measurable source parameters ${\bf u}=\{M,\,q,\,z,\,\chi_{\rm eff}\}$ for 1g+1g, 1g+2g and 2g+2g events are compared in Fig. 2. Each row corresponds to one of the three mass distributions described in Sec. II.1.

The mass distributions have some noteworthy features. First of all, and quite obviously, 2g BHs have higher component masses. Therefore the total mass is higher when 2g BHs are present (for any given assumption on the mass distribution), and this effect is most notable for the 2g+2g distributions. Mergers also tend to increase the number of comparable-mass binaries, in part because of the fixed mass range for the component masses ($m_{i}\in[5,50]M_{\odot}$). For the “power law” mass function, the mass ratio of the 1g+2g population peaks at $q=0.5$. This is because the mass distribution of the primary BH is strongly peaked at the low end of the range (i.e., at $\sim 5M_{\odot}$), so many 2g binaries are nearly equal mass, with component masses close to $5M_{\odot}$.

Redshift distributions also follow the expected trend: most 1g+1g events occur at large redshift, whereas mergers involving one or two 2g BHs occur (on average) at smaller redshift, because there is a time delay between the formation of 1g BHs via core collapse and their subsequent merger.

The most striking differences are found in the distributions of individual spins. To better illustrate this point, in Fig. 3 we show the distribution of the individual BH spins $(\chi_{1},\,\chi_{2})$, as well as the distribution of the spin of the remnant $\chi_{\rm f}$. As discussed in Berti and Volonteri (2008), from a statistical point of view the effect of mergers is to “cluster” BH spins around $\chi_{\rm f}\sim 0.7$, quite independently of the progenitor parameters. While the 1g+1g spin magnitudes are uniform in the range $[0,\,1]$ by construction, spin distributions become peaked at $\sim 0.7$ when 2g BHs are involved. This clustering is evident in the distribution of primary spins $\chi_{1}$ for the 1g+2g and 2g+2g cases, and in the distribution of secondary spins $\chi_{2}$ for the 2g+2g case. For the 1g+2g population, the peak at $\chi_{2}\sim 0.7$ is less pronounced. This is because the lower-mass BH is most likely 1g, and the spin distribution of 1g BHs is by construction uniform in $[0,\,1]$.

Unfortunately low-SNR GW observations of merger events are not very sensitive to $\chi_{1}$ and $\chi_{2}$, but rather to the effective spin $\chi_{\rm eff}$ defined in Eq. (6). The right column of Fig. 2 shows that the effect of mergers is considerably smeared out in $\chi_{\rm eff}$, but more binaries with $\chi_{\rm eff}\sim 0$ are expected if all sources are 1g BHs. Measurements of $\chi_{\rm eff}$ may still be sufficient to distinguish between different populations, especially when comparing 1g+1g against either 1g+2g or 2g+2g. Discriminating between BH progenitors should be considerably easier with future detectors, when high-SNR events will allow for more precise measurements of $\chi_{1}$, $\chi_{2}$ and $\chi_{\rm f}$ Berti et al. (2016); Vitale (2016); Vitale and Evans (2017) .

II.6 Single detections

In the rest of this paper we will study how statistical inference from several detections can be used to constrain the underlying BH population. However, it is possible that single detections with specific parameters can already provide smoking gun evidence for the occurrence of multiple mergers.

One possibility, as mentioned in Sec. II.1, is that pair instabilities may prevent the formation of 1g BHs with masses sensibly above $\sim 50M_{\odot}$ Heger and Woosley (2002); Woosley (2017); Belczynski et al. (2016d). If this is indeed the case, a single detection of a merging BH binary where one of the components has mass larger than $50M_{\odot}$ would indicate the occurrence of multiple mergers. This argument, however, relies on two crucial assumptions: (i) that 1g BHs always form from stellar collapse, while more exotic formation channels (e.g. involving primordial BHs) may produce massive BHs without invoking multiple mergers; (ii) that pair instabilities in core collapse do indeed prevent the formation of massive BHs. Pair instabilities, pair instability pulsations and the exact value of the maximum BH mass that can be produced via core collapse are all topics of current research Belczynski et al. (2016d).

Another possibility involves accurate measurements of the component spins through the detection of a single nearby, non face-on binary merger with comparable, low masses and many precession cycles in the LIGO band. Unfortunately, parameter estimation studies suggest that current-generation detectors could allow dimensionless spin measurement errors $\sim 0.3$ in best-case scenarios Vitale et al. (2017a). Errors of this magnitude are comparable to the width of the peaks in the spins distributions shown in Fig. 3 and there is significant uncertainty in the spin magnitude distribution of astrophysical BHs, so it seems unlikely that single spin measurements may allow us to tell apart 1g BHs from 2g BHs, at least in the near future.

III Observable distributions

From the theoretical distributions described in Sec. II, we construct observable distributions by (i) selecting detectable binaries according to a detection statistic, such as a threshold in the signal-to-noise ratio (SNR), and (ii) folding in measurement errors.

III.1 Detection probability

We first assign a detection probability ${\kappa}^{(i)}<1$ to each binary in our catalogs. This number takes into account the detector sensitivity and antenna pattern, as well as the (random) sky position of the source. We compute ${\kappa}^{(i)}$ following the procedure outlined in Ref. Dominik et al. (2015), where an astrophysical catalog of binaries produced using the Startrack population synthesis code was filtered to produce similar catalogs of observable binaries for a specific set of GW detectors. This procedure is briefly reviewed below.

Each binary produces a GW strain $h(t)$ and an expectation value for the SNR

[TABLE]

where $S_{n}(f)$ is the noise power spectral density of the detector and $\tilde{h}(f)$ is the Fourier transform of the strain $h(t)$. The strain is computed using the IMRPhenomC waveform model Santamaría et al. (2010). In this paper we consider noise power spectral densities for the first AdLIGO observing run (O1), the Advanced LIGO design sensitivity Abbott et al. (2016f) (LIGO Scientific Collaboration and Virgo Collaboration), A+ (Advanced LIGO with squeezing) and Voyager (the most advanced instrument that can be hosted in facilities similar to LIGO) Mc Clelland et al. (2015).

For any binary in our catalog we can compute $\rho_{\rm opt}$, i.e. the single-detector SNR for a binary that is optimally located and oriented in the sky. We then select those binaries in the catalog that are above a detection threshold $\rho_{\rm opt}\geq\rho_{\rm thr}=8$. This criterion has often been used as a simple, reasonable proxy for a more realistic calculation of GW detection rates in multi-detector networks Abadie et al. (2010) (LIGO Scientific Collaboration and Virgo Collaboration); Dominik et al. (2015). Then we compute the detection probability as

[TABLE]

where the function $P(w^{(i)})$ is the cumulative distribution function for the projection parameter $w^{(i)}\equiv\rho_{\rm thr}/\rho^{(i)}_{\rm opt}$. This cumulative distribution function takes into account the geometrical “peanut factor” that characterizes the sensitivity of the detector to the source sky location, inclination and polarization (see Dominik et al. (2015) and references therein). Roughly speaking, $w^{(i)}=1$ means that the source is in a “blind spot” of the detector, while $w^{(i)}=0$ in the high-SNR limit. A tabulated version of $P(w^{(i)})$ is publicly available.333www.phy.olemiss.edu/$\sim$berti/research We use standard spline interpolation to compute this function for generic values of $w^{(i)}$.

III.2 Measurement errors

Ideally we should compute measurement errors for each binary in the catalog using Markov-Chain Monte Carlo methods, and use the obtained posteriors to perform model selection. This is computationally expensive, and unnecessary from the point of view of our proof-of-principle analysis. For our present purpose we adopt a much simpler prescription, described below.

We build on a study by Ghosh et al. Ghosh et al. (2016), who computed BH binary measurement errors using the lalinference code Veitch et al. (2015) (see also Ajith and Bose (2009); Littenberg et al. (2013); Cho et al. (2013); Graff et al. (2015) for more work on the subject). In particular, we use their results for aligned-spin BH binaries detected by a network of 3 advanced detectors. Their data set provides 1$\sigma$ errors on several quantities, including the total mass $M$, mass ratio $q$ and redshift $z$. These are shown in blue in Fig. 4.

The data set is too sparse to perform an efficient binning and interpolation in three dimensions ($M,\,q,\,z$). In order to partially account for the expected degeneracies (e.g., close binaries will generally have smaller errors on the masses), we adopt the following procedure. Consider a binary in our catalog with parameters $({\bar{M},\,\bar{q},\,\bar{z}})$. To estimate measurement errors on the parameters of this binary, we consider the 5 “closest” binaries in the data set of Ref. Ghosh et al. (2016), and compute the average and standard deviation of their measurement errors. Here “closest” is defined in the following sense: given the maximum and minimum value of each of the three parameters ($M,\,q,\,z$), we rescale their actual values so that these parameters are distributed in a cube of size one; then we compute the Euclidean distance between binaries in this cube. The average and standard deviation from the 5 closest binaries are then used to extract the measurement errors ${\sigma_{\bar{M}},\sigma_{\bar{q}},\sigma_{\bar{z}}}$ from a normal distribution. The red dots and histograms in Fig. 4 show the measurement errors obtained from this resampling. The obtained distributions look remarkably close to the original data. Errors on the redshift are slightly overestimated, so (if anything) our resampling procedure seems to yield conservative predictions. Estimates for the errors on $\chi_{\rm eff}$ were not computed in Ref. Ghosh et al. (2016), so we assume $\sigma_{\chi_{\rm eff}}=0.1$ for all binaries measured by LIGO at design sensitivity. This rough estimate is quite conservative, and it is consistent with measurement errors in the first GW detections Abbott et al. (2016c) (LIGO Scientific Collaboration and Virgo Collaboration).

Ref. Ghosh et al. (2016) computed parameter estimation errors for the LIGO-Virgo network at design sensitivity. Fisher matrix arguments Poisson and Will (1995) suggest that the capabilities of other detectors can be estimated rescaling the errors on the total mass, mass ratio, luminosity distance and $\chi_{\rm eff}$ by the ratio of SNRs, i.e.

[TABLE]

Luminosity distance and redshift are related by

[TABLE]

(where we use units such that $c=1$), so that

[TABLE]

The error on the redshift $\sigma_{z}$ is related to the error on $D_{L}$ by

[TABLE]

where we assumed $\sigma_{\Omega_{\Lambda}}\simeq 0$ (see e.g. Hughes (2002); Berti et al. (2005)). Given recent discrepancies in the determination of $H_{0}$, we assume ${\sigma_{H_{0}}}/{H_{0}}=0.1$ Di Valentino et al. (2016); Riess et al. (2016); Bernal et al. (2016).

III.3 Binning

Recall that each binary is characterized by a vector of observable parameters ${\bf u}=\{u_{1},\dots,u_{J}\}$. If (for simplicity) we momentarily neglect measurement errors, the observable distribution is just a sum of Dirac deltas centered at $\bar{\bf u}^{(i)}$, and each delta is weighted by the detection probability ${\kappa}^{(i)}$:

[TABLE]

where $\lambda$ labels the model (cf. Sec. IV) and the denominator ensures normalization. Using the procedure described in Sec. III.2 we can obtain estimates of the 1$\sigma$ errors on the measurement of each parameter. The $i$th binary in the catalog now has estimated parameters $\bar{u}^{(i)}$ with errors $\sigma^{(i)}=\sigma({\bar{u}^{(i)}})$. Assuming that errors are normally distributed and neglecting degeneracies, we can substitute the Dirac deltas of Eq. (14) with Gaussian distributions:

[TABLE]

where

[TABLE]

Next, we need to bin the distributions $\tilde{r}({\bf u},\lambda)$. In each direction $j$, we construct bins $k_{j}$ with extrema $b_{k_{j}}$ and $B_{k_{j}}$, i.e. $u_{j}\in(b_{k_{j}},B_{k_{j}})$. The function $\tilde{r}({\bf u},\lambda)$ in each multi-dimensional bin $\{k_{1},...,k_{J}\}$ is given by the integral

[TABLE]

In practice, we spread each source over multiple bins because of measurement errors (see Gair et al. (2010); Sesana et al. (2011) for a similar approach in the LISA context). Eq. (17) is correctly normalized to $1$ only if the bins $k_{j}$ span the entire support of $\tilde{r}(\mathbf{u},\lambda)$. When substituting Dirac deltas with Gaussian distributions we are adding support in the whole range $[-\infty,+\infty]$ for each of the $u_{j}$’s, and inevitably we end up using a finite range. For simplicity, we just renormalize $\tilde{r}_{k_{1},...,k_{J}}(\lambda)$ such that

[TABLE]

From now on we will identify the bins by a multi-index variable $k=\{k_{1},...,k_{J}\}$, so (for example) we can write $\sum_{k}f_{k}\equiv\sum_{k_{1}}...\sum_{k_{J}}f_{k_{1},...,k_{J}}$ for any binned quantity $f$.

III.4 Putting the pieces together

Examples of observable distributions are given in Fig. 5 for Advanced LIGO at design sensitivity (top) and Voyager (bottom) assuming the “flat” mass function. In each panel, dashed lines show the theoretical distribution for the 1g+1g, 1g+2g and 2g+2g populations, as already presented in Fig. 2. The histograms show the observable population, i.e. the distribution of detectable binaries, where the measured parameters take into account also measurement errors. Some trends are visible.

Let us first focus on the top row, which refers to observations with Advanced LIGO at design sensitivity. It is clear that binaries with larger total mass and lower redshift produce stronger signals, and therefore they are more likely to be detected. In particular, Advanced LIGO can hardly detect any binaries at redshift $z\gtrsim 1$. The distribution of $\chi_{\rm eff}$ also shows a mild excess of observable events with $\chi_{\rm eff}\simeq 0$ for the 1g+1g population with respect to the 1g+2g and 2g+2g populations, suggesting that measurements of $\chi_{\rm eff}$ can indeed help to discriminate between populations.

The bottom row of Fig. 5 shows that the increased sensitivity of a Voyager-like detector has two main effects: it makes observable distributions in each of the parameters much closer to the corresponding theoretical distributions, and (quite importantly) it extends the reach of the detector to high-$z$ binaries. We obviously expect that more sensitive detectors will allow better discrimination between the different populations.

The 2g+2g population presents a peak at $M\sim 80M_{\odot}$ and $q\sim 1$. Equal-mass binaries of $\sim 40M_{\odot}+40M_{\odot}$ can only be detected by Advanced LIGO at design sensitivity if they are located at very small redshift (cf. e.g. Abbott et al. (2016d) (LIGO Scientific Collaboration and Virgo Collaboration)). This explains the significant drop in the number of observed events as $q\to 1$. The effect is strongly mitigated in Voyager, because the instrument is more sensitive at low frequency.

IV Statistical tools

In this section we briefly introduce statistical tools to perform Bayesian model selection. We label models by a parameter $\lambda$ that can be either discrete (if we want to distinguish two competing models $A$ and $B$) or continuous (if want to measure the “mixing fraction” between competing models that best describes the data).

IV.1 Number of observations

Our goal is to infer which model $\lambda$ best describes a set of data. As explained above, our binned distributions $\tilde{r}_{k}(\lambda)$ are normalized. To compare our models with the data we need an extra parameter $N(\lambda)$, the total number of observations predicted by model $\lambda$. We write

[TABLE]

As for the individual binary parameters, we introduce an array ${\bf d}$ whose elements are the single observations $d^{(i)}$, which in turn are $J-$dimensional arrays. We bin the array ${\bf d}$ on the same grid used for the catalogs to obtain binned values $d_{k}$.

The likelihood of obtaining a data set $d_{k}$ from model $\lambda$ is given by

[TABLE]

In our analysis the total number of observation does not contain information about the given model (this may not be the case for more realistic scenarios, where different models predict different merging rates: see e.g. Stevenson et al. (2015)). We therefore marginalize the likelihood over $N(\lambda)$. Plugging Eq. (19) into Eq. (20) one obtains Sesana et al. (2011)

[TABLE]

and consequently the marginalized likelihood is

[TABLE]

Note that the term $\sum_{N}\left(N^{\sum_{k}n_{k}}e^{-N}\right)$ is a multiplicative coefficient that only depends on the data ${\bf d}$, and not on the model $\lambda$. This term can be ignored because, as we will see below, we are only interested in likelihood ratios, not in the likelihoods themselves.

From now on, to simplify notation, we will drop the tilde on $p$ and assume that likelihoods are always marginalized over the total number of events.

IV.2 Model selection

Let us first look at model comparison between pure models, so that $\lambda$ is a discrete variable. Given models $\lambda=A$ and $\lambda=B$, their odds ratio is defined as

[TABLE]

where $\pi$ is the prior probability assigned to each of the two models. The simplest assumption on the priors is $\pi(A)=\pi(B)=1/2$, such that the odds ratio reduces to the likelihood ratio. If $O_{AB}\gg 1$ ($O_{AB}\ll 1$) the data favors model $A$ ($B$). The probability of model $A$ is

[TABLE]

and the probability of model $B$ is $p_{B}=1-p_{A}$. Sometimes $\sigma$-levels are used to quantify the significance of a discrete model comparison, in analogy with Gaussian measurements. The expression relating the odds ratio $\mathcal{O}$ and $\sigma$ is

[TABLE]

We can also assume that the data are represented by a mixture of two or more models, and assess whether the data themselves are informative about the underlying model mixing fractions. Each pure model $m$ enters the mixed model with a weight $f_{m}$, such that $\sum f_{m}=1$. Model comparison is equivalent to Bayesian inference on the parameters $\lambda=\{{f_{1},f_{2},\dots}\}$, as described by the posterior distribution

[TABLE]

As before, $\pi(\lambda)$ is the prior assigned to each mixed model. We choose $\pi(\lambda)$ to be uniformly distributed on the surface $\sum f_{m}=1$.444For instance, for a mixture of three models $\lambda=\{f_{1},f_{2},f_{3}\}$ the equation $\sum f_{m}=1$ describes a 2-dimensional surface $S$ of area $\sqrt{3}/2$. The uniform prior on $S$ is given by $\pi(f_{1},f_{2},f_{3})=2/\sqrt{3}$, so that $\iint_{S}\pi\,dS=1$. From a computational point of view, we first draw values of $\lambda$ from the uniform prior, and then we produce a statistical sample distributed according to $p({\bf d}|\lambda)$ using a standard Monte Carlo hit-or-miss algorithm.

V Results

So far we have outlined a procedure to build a set of “synthetic” GW observations of merging BH binaries (along with their associated errors) from simple astrophysical considerations. We now wish to understand whether these observations can be used to distinguish between different populations using Bayesian model selection (see e.g. Gair et al. (2010, 2011); Sesana et al. (2011); Berti et al. (2011); Stevenson et al. (2015) for previous studies of this problem in different contexts).

V.1 LIGO O1 data

We first apply our model comparison tool to the three LIGO O1 observations. The data set $\mathbf{d}$ consists of the maximum likelihood values provided in Ref. Abbott et al. (2016c) (LIGO Scientific Collaboration and Virgo Collaboration):

• •

GW150914:

$M=65.3M_{\odot}$, $q=0.81$, $z=0.090$, $\chi_{\rm eff}=-0.06$.

• •

GW151226:

$M=21.8M_{\odot}$, $q=0.52$, $z=0.094$, $\chi_{\rm eff}=0.21$.

• •

LVT151012:

$M=37M_{\odot}$, $q=0.57$, $z=0.201$, $\chi_{\rm eff}=0.03$.

As stressed above, measurements errors are included in this analysis at the level of the catalogs, by spreading each source over multiple bins. A more in-depth study should make use of the posterior distribution of the observed parameters obtained through dedicated parameter-estimation pipelines.

Performing model selection as described in the previous sections and using the O1 sensitivity curve, we obtain the odds ratio reported in Table 1. We also repeat the same exercise assuming the anticipated noise power spectral density of Advanced LIGO at design sensitivity. This basically answers the question: “what if the O1 observations had been carried out with a better detector?”

As shown in Table 1, most of the odds ratios are in the range $0.3\lesssim\mathcal{O}\lesssim 3$, corresponding to $1\sigma$. This simply indicates that three observations are not enough to perform a meaningful statistical analysis. However some of the comparisons return odds ratios $\mathcal{O}\sim 10$, approaching $2\sigma$ evidence. When this happens (i) the 1g+1g population seems to be preferred, and (ii) the odds become higher for a more sensitive detector like Voyager. In these cases the algorithm seems to capture real statistical differences between the catalogs, that become more pronounced when more binaries are detected and measurement errors get smaller.

As a note of caution, we stress here that such discrete model comparison analyses can only tell us which of two competing models better describes a given data set, not which model is correct. For instance, our results in Table 1 show some dependence on the underlying mass distribution. This could be due either to the low dimensionality of the statistical sample (cf. Sec. V.2 below), or to the fact that none of the three mass distributions faithfully describes the observations. To bracket uncertainties in the time delay prescription (cf. Sec. II.4), Table 1 also lists odds ratios computed omitting all redshift information. This calculation shows that assumptions on the time delays do not significantly affect our conclusions, given the limited statistics currently available. It will be straightforward to update our analysis with higher statistics and better motivated BH binary formation models when more data become available.

V.2 Simulated data: Pure models

The results of Sec. V.1 show, not surprisingly, that more than three observations are needed to discriminate between different models. In order to estimate the capabilities of larger data sets and more sensitive detectors, here we perform model selection on simulated observations. Our main goal is to estimate how many observations are needed to distinguish a pair of models with a given confidence level.

Given a model $\lambda_{\rm true}$, we extract the number of events per bin $d_{k}$ assuming a Poisson distribution

[TABLE]

Here the total number of observation $N_{\rm obs}=N(\lambda_{\rm true})$ is a free parameter that we need to specify. We expect model comparison to be easier/harder if more/less observations are available. This statement is made more quantitative in Fig. 6 and Table 2.

Figure 6 shows the odds ratio distribution obtained from several realization of $N_{\rm obs}$ observations. For each pair of models we plot $\mathcal{O}_{AB}$ (when $A$ is the true model) and $\mathcal{O}_{BA}$ (when $B$ is the true model), thus addressing how easy (or hard) it is to identify any of the models if it is correct. Thick lines mark the median odds, while the shaded areas encompass 90% of the realizations (i.e., they cover the range between the 5th and the 95th percentiles).

The odds ratio $\mathcal{O}$ increases roughly exponentially with the number of observations $N_{\rm obs}$, so our ability to distinguish between different models should rapidly improve in the coming years. Table 2 shows that in $5\%$ of the realizations, as few as $\sim 20$ detections are enough to discriminate the 1g+1g population from the 2g+2g population at $5\sigma$ with Advanced LIGO at design sensitivity, while $N_{\rm obs}\sim 80$ observations are necessary to achieve $5\sigma$ confidence in 95% of the realizations.

Model selection involving the 1g+2g population typically requires a larger number of observations. This is clear when comparing the left panels of Fig. 6 to the middle and right panels. In both the (1g+1g vs. 1g+2g) and (1g+2g vs. 2g+2g) comparisons the odds ratio grows (roughly) exponentially, but with smaller slope compared to the (1g+1g vs. 2g+2g) case. However the slope (and the odds ratio $\mathcal{O}$) is larger when 1g+2g is the true model: it is slightly easier to mistake a 1g+1g (or 2g+2g) population for a 1g+2g population than vice versa.

Model comparison is easier with more sensitive detectors. For example, distinguishing 1g+1g from 2g+2g at $5\sigma$ in 90% of the realizations requires only $\sim 30$ Voyager observations (instead of $\sim 80$ for Advanced LIGO at design sensitivity).

In Sec. V.1, where only three observations were considered, the results were greatly dependent on the assumed mass distribution. Table 2 shows that this dependence becomes much weaker when more observations are available and/or the instrumental sensitivity improves. This is largely due to the discriminating power of the redshift distribution of the events, which becomes more relevant when high-$z$ binaries become detectable (cf. Fig. 5).

V.3 Simulated data: Mixed models

Let us now turn to a more ambitious task. As anticipated in Sec. IV.2, we now consider a population of binaries consisting of a mixture of the three pure models 1g+1g, 1g+2g and 2g+2g. The task is to measure their mixing fraction, i.e. to determine how many binaries belong to each of the three pure populations. This is computationally expensive, as it requires many evaluations of the likelihood defined in Eq. (26) through Monte Carlo methods.

As a proof of principle, in Fig. 7 we show results for a specific choice of the mixing parameters. Simulated observations are drawn from a model555Our injected fraction of 2g BHs was chosen only for illustrative purposes. It is higher than current estimates of merger rates in nuclear clusters, which are favorable environments for multiple merger events Antonini and Rasio (2016). where $60\%$ of the binaries are 1g+1g, $10\%$ are 2g+2g, and $30\%$ are 1g+2g:

[TABLE]

For concreteness we assume the “flat” mass prescription and consider several realizations of $N_{\rm obs}=100$ observations performed with the Advanced LIGO network at design sensitivity. Each of the triangles in Fig. 7 shows the surface $f_{\rm 1g+1g}+f_{\rm 1g+2g}+f_{\rm 2g+2g}=1$. The color coding corresponds to the values of the posterior $p(\lambda|{\bf d})$. Pure models lie on the corners of this “Dalitz plot,” while the star marks the injected fraction.

As expected, measuring mixing fractions is sensibly harder than performing discrete model comparison, and it is going to require many more observations (a similar result was obtained in Ref. Sesana et al. (2011)). The injected fractions are recovered only in some of the realizations, suggesting that these data points are probably not enough to confidently perform the measurement.

In any case, we can note some trends. Most (but not all) of the realizations assign a rather low probability to the region where $f_{2g+2g}\sim 0$. Whenever a few 2g+2g events are present, their properties are sensibly different from those involving 1g BHs, and therefore the 2g+2g population can be identified relatively easily. Although we may be unlucky and estimate mixing fractions which are sensibly different from their true values, our model comparison algorithm does return a statistically consistent result. Out of 1000 realizations, we find that the correct mixing fraction is identified within the 50% (90%) confidence interval in $57\%$ (25%) of the cases. The relatively small number of observations is likely to be one of the main reasons for this relatively pessimistic result: if we assume $N_{\rm obs}=1000$, the correct mixing fraction is identified within the 50% (90%) confidence interval in $90\%$ (77%) of the cases.

In conclusion, this preliminary study shows that measuring mixing fractions will be challenging in the near future. Estimating mixing fractions with high confidence may require several hundreds (if not thousands) of observations.

V.4 Caveats on mass functions and time delays

We have shown that, given a sufficient number of detected events, it is possible to distinguish a given 1g+1g BH population from a variant of the same population where repeated mergers occur. Here we discuss how uncertainties in the assumed 1g+1g mass distribution and in time delay prescriptions may bias our conclusions.

In Table 3 we perform pure model comparisons between BH binary populations that differ in both merger generation (1g+1g, 1g+2g, 2g+2g) scenario and in the assumed mass distribution. As shown in Sec. V.2, the true distribution is correctly identified whenever it is among those tested. When the injected population is not among those being compared, differences in the assumed mass distribution can sometimes dominate over differences induced by the occurrence of subsequent mergers. For instance, if injected 2g+2g observations assuming the “flat” mass distributions are examined assuming a “power law” mass model, one would erroneously conclude that the observed population is 1g+1g, rather than 2g+2g.

However this should not be a problem in practice, because the mass distribution should soon be well constrained by observations. Realistic astrophysical scenarios typically predict a small fraction of multiple mergers, i.e., a small fraction of 2g events. Even remaining theory-agnostic, this anticipation is already (although inconclusively) supported by the data. Our Table 1 suggests that 1g+1g mergers may already be favored over 2g scenarios. So, in practice, there are theoretical and (hopefully soon) experimental reason to assume that the majority of detected events have a 1g+1g origin. In this very plausible scenario, the model selection procedure can be “bootstrapped” as follows:

• (i)

The mass distribution is inferred from a large enough number of detections, assuming that most events are 1g+1g;

• (ii)

This observationally inferred mass distribution can be used to replace our “toy” mass distributions (flat, log or power law) for 1g+1g BHs, and the 2g distributions can be constructed through hierarchical mergers as described earlier;

• (iii)

We can now look at all measurable properties of the population to determine whether some (presumably small) fraction of events has a 2g origin.

Table 3 shows that step (i) above does not present problems. Indeed, according to Table 3, while it is indeed possible to wrongly rule in favor of 1g+1g BHs given 2g+2g injections, the converse is unlikely: if model selection favors 2g+2g BHs, the injected data never belong to a 1g+1g population with a different mass spectrum.

To quantify the importance of time delays, we repeated all the comparisons shown in Table 2 excluding redshift information, i.e. taking ${\bf u}=\{M,\,q,\,\chi_{\rm eff}\}$ as our vector of observable quantities. We find that the correct population is always identified. The odds still grow exponentially with the number of observations, although with somewhat shallower slopes. This is expected, because the statistical analysis is performed using less information. Omitting redshift information does not significantly affect the 1g+1g vs. 2g+2g and 1g+2g vs. 2g+2g comparisons, but it plays a more important role in the 1g+1g vs 1g+2g comparisons. This is because (as illustrated in Fig. 5) the mass distributions are very similar for these populations, which are therefore harder to distinguish if redshifts are ignored. For instance, while Fig. 6 shows that $\sim 40$ observations are enough to distinguish the 1g+1g and 1g+2g “flat” populations at $3\sigma$ with LIGO in 50% of the realizations, up to $\sim 200$ events will be necessary to reach the same conclusion in the absence of redshift information.

VI Conclusions

The main result of this paper is that GW observations of merging stellar-mass BH binaries can be used to gather information about their progenitors. Starting from simple, physically motivated populations of “first generation” (1g) BHs born from stellar collapse, we construct populations where merging binaries include “second generation” (2g) BHs, whose masses and spins are computed using numerical relativity fitting formulas. Then we use Bayesian model selection to determine whether current or future ground-based GW interferometers can distinguish different populations. If 2g BHs occur in nature, it should be possible to recover evidence for their existence from GW data; otherwise, the data can be used to constrain astrophysical models that produce 2g BHs.

As a first application of our Bayesian model selection framework, we perform model selection using the two confirmed detections (GW150914 and GW151226) and the LVT151012 trigger from Advanced LIGO’s first observing run. It is quite remarkable that, even with only three data points, some of the comparisons show odds ratios as high as $\sim 10$ in favor of 1g BHs. As expected, model selection performance improves with more observations and more sensitive detectors. Indeed, as shown in Fig. 6, the Bayesian odds ratio for comparisons between two pure models scales (roughly) exponentially with the number of observations. Depending on the actual realization, $\sim 20$–$200$ Advanced LIGO observations at design sensitivity should allow us to discriminate which of the three populations is favored by the data at $5\sigma$ confidence level in one-to-one comparisons. Instrumental upgrades will bring this number down to $15$–$200$ observations for A+, and $5$–$100$ for Voyager.

More realistically, astrophysical populations of merging binaries will be a mixture of all three populations (1g+1g, 1g+2g, 2g+2g), and the real experimental task is to determine the relative mixing fractions. Using simulated data, we construct synthetic catalogs assuming a mixture of models for the different BH generations, and attempt to measure the mixing fractions using Bayesian inference. Our preliminary results suggest that this is a much more challenging task: recovering the mixing fractions may require several hundreds (if not thousands) of observations.

This work should be regarded as a proof-of-principle study that can (and should) be extended in several directions. Our simple models are not supposed to be astrophysically realistic: they were developed solely to show that, at least in principle, GW observations could provide information on the occurrence of multiple stellar-mass BH mergers. The inclusion of detailed spin alignment models and more realistic mass distributions (see e.g. Kovetz et al. (2017)), preferably with input from population synthesis codes, is an important topic for future investigation.

As illustrated in Sec. II.5 (see, in particular, Fig. 3), the spin magnitudes of the merging BHs are very sensitive to their merger history. This is also true for the massive BH binaries observable by LISA: see e.g. Berti and Volonteri (2008); Sesana et al. (2014). Unlike BHs born from stellar collapse, the spin distribution of postmerger BHs should be strongly peaked at $\chi_{\rm f}\sim 0.7$. In this paper we only considered measurements of the “effective spin” $\chi_{\rm eff}$, because this is the spin parameter that enters at lowest PN order in the gravitational waveform. This is a very conservative approach. As shown in Fig. 2, the “memory effect” encoded in the spin magnitudes is largely washed out in this variable. Measurements of the individual spin magnitudes should be possible by considering better waveform models or higher SNR signals: in this sense, our predictions should be regarded as conservative. Moreover, high-SNR ringdown observations will allow measurements of the final (postmerger) spin $\chi_{\rm f}$ within a few percent Berti et al. (2016). These measurements could also be used to identify the progenitors of merging BHs.

The model selection framework developed in this paper is complementary to other studies, which usually focus on discriminating specific astrophysical formation channels (e.g., field binaries vs. dynamical formation scenarios Stevenson et al. (2015); Vitale et al. (2017b); Rodriguez et al. (2016c); Nishizawa et al. (2016); Breivik et al. (2016); Nishizawa et al. (2017), but see also Haster et al. (2016) for work on intermediate-mass BHs). We focused on using statistical distributions consisting of several observations, but it is possible that single events may be smoking guns for (or against) multiple merger scenarios, at the price of making stronger assumptions on the formation mechanism of 1g BHs. For example, binaries with component masses above the pair-instability gap Belczynski et al. (2016d) can point to the occurrence of multiple mergers is we assume that 1g BHs are formed via core collapse, and if we are confident about the upper mass limit on 1g BHs set by pair instabilities. We hope that our approach will spark more studies of the astrophysical information encoded in present and future GW data sets.

While completing our study we learned that Maya Fishbach, Daniel Holz and Ben Farr have been pursuing a similar investigation Fishbach et al. (2017). Their work nicely complements our own, as they focus on the spin distributions and address the detectability of more than two generations of mergers.

Acknowledgements.

We thank Archisman Ghosh, Walter Del Pozzo and Parameswaran Ajith for sharing parameter estimation data from Ref. Ghosh et al. (2016). We also thank Daniel Holz, Maya Fishbach, Ben Farr, Leo Stein and Chris Moore for discussions. D.G. is supported by NASA through Einstein Postdoctoral Fellowship Grant No. PF6-170152 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under Contract NAS8-03060. E.B. was supported by NSF Grant No. PHY-1607130 and by FCT contract IF/00797/2014/CP1214/CT0012 under the IF2014 Program. This work was supported by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP-690904. Computations were performed on the Caltech computer cluster “Wheeler,” supported by the Sherman Fairchild Foundation and Caltech. Partial support is acknowledged by NSF CAREER Award PHY-1151197. This research made use of python packages astropy Robitaille et al. (2013) (Astropy Collaboration) and matplotlib Hunter (2007).