Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors
BBayesian Evidence for Both Astrophysical and Primordial Black Holes:Mapping the GWTC-2 Catalog to Third-Generation Detectors
V. De Luca, G. Franciolini, P. Pani,
2, 3 and A. Riotto Département de Physique Théorique and Centre for Astroparticle Physics (CAP),Université de Genève, 24 quai E. Ansermet, CH-1211 Geneva, Switzerland Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185, Roma, Italy (Dated: February 9, 2021)We perform a hierarchical Bayesian analysis of the GWTC-2 catalog to investigate the mixed scenarioin which the merger events are explained by black holes of both astrophysical and primordial origin.For the astrophysical scenario we adopt the phenomenological model used by the LIGO/Virgocollaboration and we include the correlation between different parameters inferred from data, therole of the spins in both the primordial and astrophysical scenarios, and the impact of accretion inthe primordial scenario. Our best-fit mixed model has a strong statistical evidence relative to thesingle-population astrophysical model, thus supporting the coexistence of populations of black-holemergers of two different origins. In particular, our results indicate that the astrophysical mergersaccount for roughly four times the number of primordial black hole events and predict that third-generation detectors, such as the Einstein Telescope and Cosmic Explorer, should detect up tohundreds of mergers from primordial black hole binaries at redshift z (cid:38) . I. INTRODUCTION
The several detections of the Gravitational Waves (GWs) coming from black hole (BH) mergers during thefirst three observing runs performed by the LIGO/Virgo Collaboration (LVC) [1, 2] has ignited the interest inunderstanding the origin and the properties of the Binary Black Hole (BBH) population [3].Multiple formation scenarios have been proposed to explain the formation of the observed BBHs and theirmerger in the local universe. They can be broadly distinguished in two categories. One is the AstrophysicalBlack Hole (ABH) formation scenario (recently reviewed in Refs. [4, 5]) where the BHs form at the end of thestellar evolution in the late-time universe, forming binaries either in isolation through a common envelope phaseor dynamically through three-body encounters in dense stellar clusters. The other possibility is provided bythe Primordial Black Hole (PBH) scenario, firstly proposed in the late ’70s in Refs. [6–8], where the compactobjects are formed in the radiation dominated early-universe out of the collapse of large overdensities [9, 10], andbinaries are assembled via random gravitational decoupling from the Hubble flow before the matter-radiationequality (see Refs. [11, 12] for some recent reviews). Additionally, PBHs in different mass ranges could contributeto a sizeable fraction f PBH ≡ Ω PBH / Ω DM of the Dark Matter (DM) [13]. This has motivated recent efforts inconfronting the PBH scenario with the GW data [14–27] in search for evidences for PBHs. Also, PBHs may beresponsible for the formation of binaries like GW190521 [28] where at least one component confidently falls inthe so-called upper mass gap [29, 30] in which ABH binaries are not generally expected to be produced.One recent effort of the community has been on inferring the properties of PBHs, e.g. their mass distributionand abundance, assuming that all BBH events detected so far have a primordial origin. While this possibilityis certainly still allowed by current constraints (see, e.g., Ref. [26] for the most updated analysis using thesecond GW transient catalog (GWTC-2)), it is plausible that all the events may not be explained by meansof a single population, either astrophysical [31] or primordial [27]. In order to distinguish the contribution ofpopulations of different origin in the data, one needs to rely on different characteristic features which may bea smoking gun for either PBHs or ABHs. One promising feature of the PBH model is the peculiar evolutionof the merger rate, which is monotonically increasing with redshift. This differs from what is expected from anABH population, whose merger rate should peak around redshift of a few. This prediction can be employed byfuture third-generation (3G) detectors like the Einstein Telescope (ET) [32] and Cosmic Explorer (CE) [33] toconfidently distinguish PBHs from other astrophysical formation channels, in particular from BHs formed fromthe first stars [34–36], which might merge also at high redshift [37–39]. a r X i v : . [ a s t r o - ph . C O ] F e b As a first step towards this future direction, in this paper we aim at inferring the binary population propertiesconsidering both the PBH and ABH formation channels. We restrict our description of the astrophysical sectorto the one used by the LVC by means of a phenomenological model and mix it with the PBH model by performinga hierarchical Bayesian analysis on the GWTC-2 catalog.The purpose of this work is twofold. First, we wish to quantify what fraction of the current observed events maybe ascribed to PBHs. Secondly, we wish to use this mixing fraction to make a forecast for the expected numberof events to be observed by 3G detectors at high redshifts, as well as the corresponding redshift dependence ofthe merger rate in this mixed ABH-PBH formation scenario.Our hierarchical Bayesian analysis improves over the recent analysis of Ref. [27] by including the correlationbetween different parameters inferred from data, the role of the spins in both the PBH and ABH scenarios, andthe impact of accretion in the PBH model.The paper is organised as follows. In Sec. II we summarise the hierarchical Bayesian framework adopted inthis paper to perform the inference on the GWTC-2 catalog. In Secs. III and IV we describe the PBH andABH models considered in this paper, while also reporting both single-population inferences. Sec. V containsour main result showing the mixed population inference and Sec. VI draws the forecasted detectability of theuncontaminated PBH subpopulation at high redshift by 3G detectors. We conclude in Sec. VII with some futureprospects, whereas the Appendices are devoted to some technical aspects of both the binary detectability andthe PBH accretion model, respectively.
II. BAYESIAN INFERENCE
In this work we perform a hierarchical Bayesian inference analysis to select the best fitting models given theGWTC-2 dataset in the full PBH, full ABH and mixed ABH-PBH binary populations scenarios. In doing so, weclosely follow the procedure used in [24, 26] and references therein. In this section we provide a short summaryof the hierarchical Bayesian inference framework, highlighting the key inputs coming from both the modelsconsidered and the dataset. The knowledgeable reader may jump directly to the following sections where themodels considered in this work are presented.On general grounds, a hierarchical Bayesian analysis returns the posterior distribution of a series of modelhyperparameters λ based on the observed dataset d . In the following, we will also define the intrinsic eventparameters, such as masses and redshift of the binary mergers, as θ . In practice, the posterior distribution canbe inferred from the data by sampling over the hyperparameter space λ the function p ( λ | d ) ∝ π ( λ ) (cid:90) p ( d | θ ) p pop ( θ | λ )d θ , (2.1)where p ( d | θ ) is the single-event likelihood, π ( λ ) is the prior on the hyperparameters, and p pop ( θ | λ ) is the popula-tion likelihood , equivalent to a prior on the intrisic parameters θ but parametrized by the model hyperparameters λ . The parameters describing single events (e.g. masses, redshifts) are also referred to as event parameters ,while the hyperparameters describing the entire sample (e.g. the fraction of DM in PBHs) can be also referredto as population parameters . For clarity, in Table I we provide a summary of all the parameters defined in thispaper.In order to speed up the computation, such an integration is written as a weighted average over the posteriorsamples. The likelihood function is then computed according to the following equation p ( λ | d ) = π ( λ ) e − N det ( λ ) [ N ( λ )] N obs N obs (cid:89) i =1 S i S i (cid:88) j =1 p pop ( j θ i | λ ) π ( j θ i ) , (2.2)where N obs is the number of events in the catalog, i.e. the number of detections, N ( λ ) is the number of events inthe model characterised by the set of hyperparameters λ , N det ( λ ) is the expected number of observable events inthe model characterised by the hyperparameters λ computed by accounting for the experimental selection bias,and S i is a normalisation factor which depends on the length of the posterior dataset of each event in the catalog.Finally, π ( θ ) is the prior on the binary parameters used by the LVC when performing the parameter estimation. TABLE I: Relevant parameters in the analysis.
Event parameters θ m Source-frame primary mass m Source-frame secondary mass χ eff Effective spin z Merger redshiftHyperparameters λ of the PBH model M c Peak reference mass of the lognormal distribution σ Variance of the lognormal mass distribution f PBH
Fraction of PBHs in DM at formation z cut-off Accretion cut-off redshiftHyperparameters λ of the ABH model R Integrate merger rate at z = 0 κ Merger rate evolution in redshift, defined as R ( z ) ∝ (1 + z ) κ α Exponent of the total mass factor in the differential rate β Exponent of the symmetric mass ratio factor ζ ABH mass function power law scaling ψ ( m ) ∝ m − ζ m min Minimum ABH mass (step-like cut-off) m max Maximum ABH mass (step-like cut-off) µ χ eff Mean value of χ eff in the “Gaussian” spin model σ χ eff Width of χ eff in the “Gaussian” spin modelInferred parameters in the mixed model f P ≡ N detPBH / ( N detABH + N detPBH ) Fraction of PBH binaries f A ≡ N detABH / ( N detABH + N detPBH ) ≡ − f P Fraction of ABH binaries
This prior is removed to extract the values of the single-event likelihood ensuring only the informative part ofthe event posterior is used.We remark that the last sum in Eq. (2.2) gives the average of the binary parameter distribution over posteriorsamples. The factors proportional to p ( λ | d ) ∝ e − N det ( λ ) N ( λ ) N obs (2.3)characterize an inhomogeneous Poisson process [40–43] and are responsible for introducing the rate informationin the inference.In our analysis we make the following assumptions for the relevant quantities:• the prior π ( λ ) on the hyperparameters λ of the model is considered flat;• for a given model M with hyperparameter λ , p pop ( θ | λ ) is the distribution of the BBH parameters, takento be θ i = { m , m , z, χ eff } , where m i is the source-frame mass of the i -th binary component, z is themerger redshift, and χ eff ≡ χ cos α + qχ cos α q (2.4)is the effective spin parameter, which is a function of the mass ratio q ≡ m /m , of both BH spinmagnitudes χ j ( j = 1 , , with ≤ χ j ≤ ), and of their orientation with respect to the orbital angularmomentum parametrised by the angles α j . We neglect the precession spin component χ p in the inference,since this parameter is still poorly reconstructed for most of the GW events [44];The binary parameter distributions in a given model (either primordial or astrophysical), can be computedfrom the differential merger rate d R/ d m d m d z as p pop ( θ | λ ) ≡ N ( λ ) (cid:20) T obs
11 + z d V d z d R d m d m ( θ | λ ) (cid:21) , (2.5)in terms of the observation time T obs , while the number of expected detections are defined as the integral of thedifferential merger rate as N det ( λ ) ≡ T obs (cid:90) d m d m d z p det ( m , m , z ) 11 + z d V d z d R d m d m ( m , m , z | λ ) , (2.6)accounting for the selection bias by introducing the probability of detection p det as defined in Appendix A. Ascustomary, the prefactor / (1 + z ) accounts for the clock redshift at the source epoch and d V / d z stands for thedifferential comoving volume factor, see for example [45].The analysis is performed by sampling the likelihood function (2.2) in the hyperparameter space by using theMarkov chain Monte Carlo software algorithm emcee [46].In the following we will compare how various models are able to explain the GWTC-2 dataset. For a quan-titative assessment, we shall use the statistical evidence Z . Given a model M , the evidence is defined as themarginal population likelihood computed as the integral of the population posterior p ( λ | d ) , i.e. Z M ≡ (cid:90) d λ p ( λ | d ) . (2.7)In other words, the evidence is a measure of the support for a given model given the data d . One can thencompare different models by computing the so-called Bayes factors B M M ≡ Z M Z M . (2.8)According to Jeffreys’ scale criterion [47], a Bayes factor larger than (10 , . , ) would imply a strong, verystrong, or decisive evidence in favour of model M with respect to model M given the available dataset.Among all the binary events included in the GWTC-2 catalog, we select the same subset events as used in theLVC population property paper [44] for a total of BBH events. Therefore, we discard events with large FAR(false-alarm rate), i.e. GW190426, GW190719, GW190909, and events where the secondary binary componenthas mass smaller than M (cid:12) , i.e. GW170817, GW190425, GW190814, to avoid contamination from putativemerger events involving neutron stars. Similarly to Ref. [26], we adopt the “Overall_posterior” as providedin [48] for the GWTC-1 events, and the “PublicationSamples” in [49] for the GWTC-2 events.Finally, at variance with Ref. [27], the present analysis fully accounts for the correlation between the intrinsicevent parameters θ in the posterior data and it is not sensitive to the particular choice of single-event priorsused by the LVC when performing the parameter estimation. As shown in Refs. [50–52] (and in Refs. [53, 54]in particular for the PBH scenario), the priors may modify the interpretation of some events, also possiblyaffecting the maximum likelihood analysis as performed in Ref. [27]. Furthermore, as far as the PBH model isconcerned, we include the effects of accretion onto PBH binaries. Finally, we also include the spin informationin the inference, which is not included in Ref. [27]. We checked the accuracy of our computation of the Bayes factors by using a nested sampling algorithm, as implemented inDYNESTY [55].
III. THE PRIMORDIAL BLACK HOLE MODEL
For an extensive description of the PBH model we adopt throughout the paper, one can rely on Refs. [26, 56, 57].We consider the differential merger rate of PBH binaries formed in the early universe as given by [57] d R PBH d m i d m i d z = 1 . × Gpc yr f PBH η − i (cid:18) tt (cid:19) − (cid:18) M itot M (cid:12) (cid:19) − S ( M itot , f PBH , ψ i ) ψ i ( m i , z i ) ψ i ( m i , z i ) × exp (cid:34) (cid:90) t cut-off t i d t (cid:32) ˙ M tot M tot + 2 ˙ µµ (cid:33)(cid:35) (cid:18) η ( z cut-off ) η ( z i ) (cid:19) / (cid:18) M tot ( z cut-off ) M tot ( z i ) (cid:19) / , (3.1)with η = m m / ( m + m ) , M tot = m + m , t being the current age of the universe, and S being a suppressionfactor discussed below. In all the expressions the subscript “i” indicates quantities at formation epoch, whichfor a PBH of mass m corresponds to a typical redshift z i (cid:39) · ( m/M (cid:12) ) − / .The mass function considered in this work is a lognormal mass function at high redshift of the form ψ i ( m | M c , σ ) = 1 mσ √ π exp (cid:20) − ln ( m/M c )2 σ (cid:21) , (3.2)which typically arises when PBHs are formed from the collapse of large overdensities in the radiation domi-nated early universe with perturbations characterised by a symmetric peak in the comoving curvature powerspectrum [58, 59]. Other possibilities can be considered, see for example Refs. [24, 27, 60].The cut-off redshift z cut-off is a hyperparameter of the accretion model which accounts for the uncertaintyaffecting the strength of accretion close to the reionization epoch, see discussion in Appendix A. In this work wefollow the approach used in Ref. [26] and consider z cut-off as a free hyperparameter of the model. For each valueof z cut-off there exists a one-to-one correspondence between the initial and final masses, m i → m , which can becomputed according to the accretion model described in details in Refs. [57, 61] and reviewed for completeness inAppendix B. We highlight for clarity that a lower cut-off is associated to stronger accretion and vice-versa. Valuesaround z cut-off (cid:39) corresponds to negligible accretion in the mass range of interest for the LVC observations.As the spin of PBHs formed from the collapse of density perturbations in the early universe is expectedto be negligible at formation [62, 63], a non-zero effective spin parameter χ eff can only be acquired by PBHbinaries through an efficient phase of accretion [57, 61], see Appendix B. As such, a peculiar characteristics ofthe PBH model is the correlation between large values of binary total masses and large values of the spins. Thiscorrelation will be considered fully when computing the event parameter distributions in Eq. (2.5). The PBHspin directions are not correlated and randomly distributed on the two-sphere.We also highlight the enhancing contribution coming from the terms in the second line of Eq. (3.1) whichaccounts for both the hardening of binaries and the reduction of merger timescale through GW emission due tomass accretion of the PBHs participating in the binaries.Let us also stress the presence of the suppression factor S ( M itot , f PBH , ψ ) ≡ S × S , which accounts for boththe suppression of binary formation due to the surrounding DM matter inhomogeneities (not in the form ofPBHs) and the disruption of binaries due to neighbouring PBHs [25, 27, 56, 64–68]. In particular, the secondpiece S specifically accounts for disruption of PBH binaries in early sub-structures and clusters throughout thethe history of the universe as binaries typically form before matter-radiation equality. The two pieces read [27] S ( M tot , f PBH , ψ ) ≈ . (cid:20) (cid:104) m (cid:105) / (cid:104) m (cid:105) ¯ N ( y ) + C + σ M f PBH (cid:21) − / exp (cid:2) − ¯ N ( y ) (cid:3) with ¯ N ( y ) ≡ M tot (cid:104) m (cid:105) f PBH f PBH + σ M (3.3) S ( x ) ≈ min (cid:2) , . · − x − . exp (cid:0) .
03 ln x (cid:1)(cid:3) with x ≡ ( t ( z ) /t ) . f PBH . (3.4)Notice also that, for f PBH (cid:46) . , one always finds S (cid:39) , i.e. the suppression of the merger rate due todisruption inside PBH clusters is negligible. This is also supported by the results obtained through a cosmologicalN-body simulation finding that PBHs are essentially isolated for a small enough abundance [69]. In Eq. (3.3)the constant C is defined as in Eq. (A.5) of Ref. [27].The hyperparameters of the PBH model are { M c , σ, f PBH , z cut-off } and the priors chosen to perform the analysisare reported in Table II. In Fig. 1 we report the result of the Bayesian inference using GWTC-2 dataset and M c = 18 . +1 . − . . . . . . σ σ = 0 . +0 . − . − . − . − . − . l og f P B H log f PBH = − . +0 . − .
12 16 20 24 28 M c z c u t − o ff . . . . . σ − . − . − . − . log f PBH
12 16 20 24 28 z cut − off z cut − off = 21 . +1 . − . GWTC-2PBH no spinPBH spin
FIG. 1:
PBH population inference using the GWTC-2 catalog. Blue and red posteriors correspond to the case when the spininformation is used in the inference or not respectively. This result reproduces the one published in Ref. [26].
TABLE II: Prior ranges for the hyperparameters of the PBH model. We assume uniform distributions.
Parameter z cut-off M c [ M (cid:12) ] σ log f PBH
Prior range [10 ,
30] [5 ,
40] [0 . , .
1] [ − , − . assuming all events have primordial origin. The result is consistent with what found in Ref. [26] by means ofmachine learning empowered hierarchical Bayesian analysis and it is reported here in order to allow for a fullcomparison with the mixed ABH-PBH case.The most striking feature of the posteriors for the full PBH analysis is given by the preference for a large cut-off redshift when no spin information is used. This is the result of the fact that, given accretion is a non-linearprocess which enhances the high mass tail of the binary mass function, lower cut-off are typically associated toa larger fraction of events with masses in the heavy portion of the observable window at odds with the dataset.However, this tendency is inverted when the spin information is used, by the need to explain the O (10) events inthe GWTC-2 catalog with a spin incompatible with zero. Those can only be accommodated in the PBH modelif some accretion is present. IV. THE LVC ASTROPHYSICAL PHENOMENOLOGICAL MODEL
We consider the representative astrophysical phenomenological model dubbed “Truncated Model” adopted bythe LVC in the population analysis of both the GWTC-1 [70] and GWTC-2 [44] catalogs. This model considers
TABLE III: Prior ranges for the hyperparameters of the ABH model. We assume uniform distributions.
Parameter R [ Gpc − yr − ] κ α β ζ m min [ M (cid:12) ] m max [ M (cid:12) ] µ χ eff σ χ eff Prior range [1 , [2 ,
10] [30 , [-1,1] [0,1] the local merger rate R , the mass spectrum slope ζ , both the minimum m min and maximum m max mass, andthe exponent β of the symmetric mass ratio as free parameters. The corresponding differential merger rate canbe parametrised as d R ABH d m d m d z = N R (1 + z ) κ ( m + m ) α (cid:20) m m ( m + m ) (cid:21) β ψ ( m ) ψ ( m ) , (4.1)where N is a normalization factor ensuring the local merger rate is R ( z = 0) ≡ R . Also, the mass functionconsidered is a power law with sharp cut-off at both ends as ψ ( m | ζ, m min , m max ) ∝ m − ζ for m min < m < m max , (4.2)with an overall constant enforcing unitary normalization as (cid:82) ψ ( m )d ln m = 1 .In the recent analysis by the LVC they adopted three different spin models named “Default”, “Gaussian” and“Multi-” spin model, respectively, see Appendix D of Ref. [44] for more details. In this work we adopt the“Gaussian” model and, as mentioned in Sec. II, we neglect the precession spin parameter χ p in the inference,since it is still poorly measured for most of the GW events [2]. In the “Gaussian” model the effective spinparameter χ eff is distributed as a Gaussian function with mean µ χ eff and variance σ χ eff truncated within therange [ − , . Notice that the population reconstructed using either the “Gaussian” and “Default” models are inagreement with each other in the analysis performed by the LVC [44].The choice of the hyperparameters for the phenomenological ABH model is summarized in Table III alongwith the ranges for the priors considered in the analysis. Notice that in the ABH model we fix the evolution ofthe merger rate with redshift to the best fit value found in [44]. We stress however that, being the data onlylimited to redshifts smaller than unity, the result of the inference is still largely insensitive to the merger rateevolution.In this section, we consider the scenario in which all events in the GWTC-2 catalog (as specified in Sec. II)are interpreted as ABH binaries and the resulting posterior distribution is shown in Fig. 2. Being the overallrate consistent with the one inferred by the LVC, R = 15 . +5 . − . / Gpc yr , one relevant feature of the posteriordistribution is its sharp drop of the posterior for values of m max smaller than (cid:46) M (cid:12) . This is due to the lack ofsupport of the model where the posterior of the primary mass of the mass gap event GW190521 falls [28]. Thisalso follows from the choice of a hard cut-off at the maximum mass in the ABH model. An analogous behaviouris observed for m min . It is interesting to notice the negative correlation between R and ζ . This can be explainedas follows. With a steeper mass function the binary masses are expected to be on the lighter portion of theobservable mass range. However, having a smaller SNR, lighter binaries are harder to observe. Therefore, therate must be enhanced in order to predict enough observable events for the model to be compatible with theGWTC-2 catalog. A similar correlation between the slope ζ and the maximum mass is observed since if the massfunction extends to larger masses, the posterior moves to steeper profiles in order to reduce the contribution ofmassive events in the observable range.Furthermore, the preference for a large value of β signals a tendency towards symmetric binaries. This isintuitively understood by noticing that the vast majority of events in the GWTC-2 catalog is compatible with q ∼ . This feature was already observed using the GWTC-1 catalog. Finally, the effective spin parameterdistribution is found to be narrow and slightly biased towards positive χ eff . Even though there is no correlationbetween the spin distribution and the mass distribution in the ABH merger rate model in Eq. (4.1), there existsa well-known parameter degeneracy between ( χ eff , q ) in the GW data [71–76]. We checked that the reconstructedpopulation parameters are only slightly affected by the introduction of spin information in the inference in thiscase. R = 15 . +5 . − . . . . . β β = 10 . +1 . − . − . − . − . − . − . ζ ζ = − . +0 . − . m m i n m min = 8 . +0 . − . m m a x m max = 67 . +6 . − . . . . . µ χ e ff µ χ eff = 0 . +0 . − . R . . . . σ χ e ff . . . . β − . − . − . − . − . ζ m min
64 72 80 88 96 m max .
04 0 .
08 0 .
12 0 . µ χ eff .
04 0 .
08 0 .
12 0 . σ χ eff σ χ eff = 0 . +0 . − . GWTC-2ABH spin
FIG. 2:
ABH population posterior found by running the inference on the GWTC-2 catalog.
V. MIXED PBH + ABH SCENARIO
In this section we investigate the mixed scenario where the GWTC-2 population is supposed to contain both ABHand PBH binaries and perform the inference accounting for possible contributions to the observed binaries fromboth models. Similar analysis accounting for multiple astrophysical formation channels, but without accountingfor a possible PBH contribution, can be found in Refs. [31, 37, 77].Consistently with the previous sections, the set of hyperparameters of the mixed model contains astrophysicalinputs { R , β, ζ, m min , m max , µ χ eff , σ χ eff } as well as primordial parameters { M c , σ, f PBH , z cut-off } . The result ofthe inference is shown in Fig. 3.The first noticeable feature is the reduction of the values of both f PBH and R in the mixed scenario comparedto the corresponding single-population scenarios studied in the previous sections. As both parameters are directlyrelated to the overall merger rate in both the ABH and PBH sectors respectively, their values are reduced asonly a smaller portion of events can be ascribed to either sectors. In particular one can compute the mixing R = 12 . +4 . − . β β = 7 . +3 . − . − − − − ζ ζ = − . +0 . − . m m i n m min = 8 . +1 . − . m m a x m max = 41 . +8 . − . . . . . . µ χ e ff µ χ eff = 0 . +0 . − . . . . . . σ χ e ff σ χ eff = 0 . +0 . − . M c M c = 29 . +6 . − . . . . . . σ σ = 0 . +0 . − . − . − . − . − . − . l og f P B H log f PBH = − . +0 . − . R z c u t − o ff β − − − − ζ m min
45 60 75 90 m max .
00 0 .
05 0 .
10 0 .
15 0 . µ χ eff .
04 0 .
08 0 .
12 0 .
16 0 . σ χ eff
12 18 24 30 36 M c . . . . . σ − . − . − . − . − . log f PBH
12 16 20 24 28 z cut − off z cut − off = 23 . +4 . − . GWTC-2ABH spin + PBH spin
FIG. 3:
Full PBH+ABH population inference using the GWTC-2 catalog. fractions in the GWTC-2 events, defined as f P ≡ N detPBH / ( N detABH + N detPBH ) . (5.1) f A ≡ N detABH / ( N detABH + N detPBH ) ≡ − f P . (5.2)In Fig. 4 we show both the distributions of the mixing fractions. The ABH model confidently outperforms thePBH one, having a fraction of events peaking at around f A (cid:39) . , corresponding to about four times the numberof PBHs in the GWTC-2 catalog. This indicates a preference of the data for the ABH phenomenological modeland a smaller contribution to the GWTC-2 catalog from PBH binaries. It is important to stress though thatthe posterior distribution for the mixing fraction of PBH f P is incompatible with zero. The secondary peakstructure in f P (cid:39) / is due to a subdominant region in the posterior with small m max along with large f PBH .0 FIG. 4:
Distribution of the mixing fraction between ABHs and PBHs as found by a Bayesian inference on the GWTC-2 catalogwith the mixed ABH-PBH model.
TABLE IV: Bayesian evidence ratios.
Model M PBH ABH PBH+ABH log B M PBH + ABH ≡ log ( Z M /Z PBH + ABH ) -4.91 -2.17 0 Overall, there is a statistically significant portion of the posterior for which the mixing fractions are comparable.We also perform a comparison of how well the single-population (either ABH or PBH) scenario is able toexplain the entirety of GW events detected so far as compared to our best-fit ABH+PBH mixed model. InRefs. [24] and [27] using the GWTC-1 and GWTC-2 dataset respectively, it was reported that the LVC phe-nomenological model was strongly favoured when compared to the PBH scenario. This conclusion is confirmedby our analysis, as the Bayes factor (see Table IV) turns out to be log B ABHPBH ≈ . , which implies a decisivepreference for the ABH model. This shows that the inclusion of accretion in the PBH model, as well as spininformation in the inference, is not able to reduce the gap between the ABH and PBH models in the attemptto explain the LVC observation by means of a single population. However, we find that log B ABH + PBHABH ≈ . ,providing decisive evidence (according to Jeffreys’ criterion [47]) for our best-fit mixed ABH-PBH model relativeto a single-population ABH model. This implies that the presence of some PBH events in the GWTC-2 catalogseems to be demanded to better fit the data.Overall, the mass distributions found for both models are generally similar to the one observed in the isolatedscenarios. One important exception is provided by the preferred value of the upper mass cut-off m max in the ABHsector. Indeed, compared to the single-population scenario, the upper mass bound is found to be significantlysmaller being m max = 41 . +8 . − . M (cid:12) . This is due to the strong preference of the inference to interpret massiveevents as primordial binaries, see Figs. 5 and 6 and their discussion below. As a consequence, also the massfunction is found to be slightly less steep, due to the negative correlation between ( m max , ζ ) already discussed inthe previous section. In the PBH sector, one relevant difference is in the reduced value of f PBH , which confirmsthat only a subdominant portion of the DM is allowed to be in the form of PBHs. Also, the preferred PBH massfunction is shifted towards larger masses (with a larger reference mass scale M c ) and a narrower shape, as itonly needs to accomodate a smaller portion of the GWTC-2 events compared to the single-population scenario.Finally, the mixed population inference confirms a preference for weak accretion with a value of z cut-off which iscompatible with the one found in Sec. III. The posterior on z cut-off is broader than in the single-population PBHcase, since in the mixed scenarios PBHs are subdominant and the hyperparameters of their population are lessconstrained by the data.In order to guide the understanding of how the GWTC-2 events are interpreted by the inference, in Fig. 5we show the ratio of the contribution to the likelihood function L ≡ p ( λ | d ) /π ( λ ) in Eq. (2.1) from both theastrophysical and primordial sectors, averaged over the population posterior. The plot clearly shows that events1 G W G W G W G W G W G W G W G W G W G W - - G W G W G W G W * G W G W G W G W G W G W G W G W G W * G W G W G W G W G W G W G W G W G W G W G W G W G W G W G W G W * G W G W G W G W G W FIG. 5:
Plot of the ratio between the ABH and PBH likelihood for each event in the GWTC-2 catalog. The ratio is averaged overthe population posterior in Fig. 3. This plot indicates that the inference confidently identifies massive events in the GWTC-2catalog as PBHs. The labels with an asterisk correspond to GW190413*=GW190413_134308, GW190521*=GW190521_074359and GW190828*=GW190828_065509 respectively. with primary masses which have support in the astrophysical mass gap (i.e. above O (45) M (cid:12) ) are confidentlyassigned to the PBH population.The same information is contained in Fig. 6, where we show the distribution of primary mass m and massratio q for both cases in both sectors separately and combined. The most striking difference between thetwo sectors is the presence of a tail at large masses well within the upper mass gap in the PBH binary m distribution. On the other hand, the distribution of the mass ratio is found to be similar in both sectors, witha slightly enhanced ability of the PBH model to accomodate small mass ratios. We again conclude that thefraction of PBH binaries in a mixed model could be responsible for events in the GWTC-2 catalog with largemasses which do not fit well within the ABH model.It is interesting to compare our findings to the population analysis performed by the LVC collaboration. Asalso found in this work, Ref. [44] concluded that events in the GWTC-2 catalog are suggesting the presenceof an additional distinct population of binaries with primary masses above ≈ M (cid:12) . Within the standardstellar formation scenario, it is difficult to accomodate massive events due to the Pulsational Pair SupernovaInstability (PPSN) preventing the formation of binaries with masses above ≈ M (cid:12) , even though the preciselocation of such a cut-off is still uncertain [78–91]. One interesting possibility within the astrophysical sector isthat massive events in the catalog are coming from second generation mergers in globular clusters or galacticnuclei [92–98]. Our findings show that introducing a PBH population of binaries in the inference naturally leadsto the interpretation of those events as coming from primordial binaries. Interestingly, the distribution of massratio reconstructed from the data is found to be similar to the astrophysical one (i.e. both showing the sametendency for symmetric binaries), in contrast to the expectation of second generation mergers which predicts asmaller mass ratio (see for example Ref. [99]). This is quite an interesting difference, which could be used in thefuture to disentangle primordial and second generation mergers in the high-mass portion of the dataset. VI. IMPLICATIONS FOR 3G DETECTORS
One of the most prominent differences between the PBH scenario and any astrophysical formation scenario isthe predicted redshift evolution of the merger rate. Even though the GWTC-2 catalog is still limited to “local”GW signals coming from sources at redshift up to z (cid:46) , future observations with larger horizons will allow touse the merger rate evolution to distinguish between different binary BH populations at current [100–102] andfuture detectors [37, 103–105].The redshift evolution of the merger rate density for the PBH model is found to be monotonically increasing2 - - - - - FIG. 6:
Top:
Contribution to the primary mass and mass ratio distributions for each sectors separately.
Bottom:
Cumulativepopulation distribution. Notice that the relative contribution of each population to the overall event parameter distribution isweighted by the overall relative number of events expected from each population as shown in Fig. 4. The bands around the dashedlines represent the confidence spread. as [56, 57] d R PBH ( z )d z ≈ t − / ( z ) , (6.1)extending up to redshifts z = O (10 ) . Notice that the evolution of the merger rate with time shown in Eq. (6.1)is entirely dictated by the dynamics of primordial binary formation (i.e. how pairs of PBHs decouple fromthe Hubble flow) before the matter-radiation equality era. As such, it is a robust prediction within the PBHscenario. On the other hand, the merger rate of ABHs is expected to peak at redshift of a few with a possiblesecond peak coming from a Pop III population at redshift z ∼ [37, 106–108]. In particular, astrophysicalmodels predicts that BHs from Pop III stars [34–36, 38] should form at z ∼ [39] and one could conservativelyassume that they form up to z = 30 . Their merger time depends on the formation mechanisms of the binariesand could range from O (Gyr) (in which case they merge at z (cid:46) ) to O (10 Myr) if they form dynamically inPop III clusters, in which case they would merge almost at the redshift of BH formation [39].Therefore, assuming the most conservative scenario (BH mergers formed from Pop III clusters), the detectionof a binary BH at redshift higher than z ∼ would be a smoking-gun signal in favour of PBHs, as no ABHsare expected at those or higher redshifts in a standard cosmology [109]. It is interesting therefore to make aforecast for the observability of possible PBH at high redshifts through 3G detectors by assuming the PBH massfunction providing the subpopulation of the currently available GWTC-2 dataset found in this paper.A description of the sensitivity of 3G detectors is reported in Appendix A. In Fig. 7 we show the predicteddistribution of observable events per year for ET and CE as a function of redshift, assuming the PBH populationidentified by the mixed population inference on the GWTC-2 catalog found in the previous section. The ETtelescope, at design configuration, will be able to observe at least one PBH event per year up to z ∼ . Asimilar figure is found for CE. In the right panel of Fig. 7 we also show the distribution as a function of theprimary mass. As expected by looking at the observable horizon shown in Fig. 8, most of the distant binarieswill have their primary mass around the characteristic PBH mass function scale M c .Finally, we also computed the total number of observable distant events (i.e. with z > ) per year for both3 -
20 40 60 80 100102030405060 - - - - - FIG. 7:
Left:
Distribution of observable events per year as a function of redshift at ET and CE coming from the PBHsubpopulation.
Right:
Distribution of observable events as a function of both primary mass and redshift at ET coming from thesubpopulation of PBHs.
ET and CE detectors, finding N ETdet ( z >
30) = 1315 +305 − / yr , N CEdet ( z >
30) = 12 +22 − / yr . (6.2)Note that the larger event rate expected for ET is due to the superior behaviour of the adopted ET-D sensitivitycurve at low-frequencies compared to the CE phase-1 design (see Fig. 8), which is particularly relevant for heavyand/or high-redshift mergers. We conclude that 3G detectors have the potential of unequivocally confirm in thefuture the presence of a PBH subpopulation in the currently available dataset and to unveil a novel (primordial)family of BBHs. However, this conclusion depends also on the measurement accuracy of the redshift, which canbe low for the most distant events even in the 3G era [37].One might think that a more conservative way of proceeding to estimate the maximum number of eventsfor the future 3G detectors would be to maximise the PBH merger rate adopting the largest value of f PBH currently allowed by the constraints. This would, in principle, allow a forecast independent from the details ofthe astrophysical model used to analyse the GWTC-2 data. However, within the PBH merger rate there existsa degeneracy in the choice of the mass distribution hyperparameters: the constraints on f PBH drastically changewhen moving in the ( M c , σ ) plane. For instance, increasing σ leads to much smaller values of f PBH allowed bythe CMB constraints (see for instance Refs. [26, 57]).
VII. CONCLUSIONS
In this work we have presented a hierarchical Bayesian investigation of the merger events in the GWTC-2catalog with the goal of understanding the nature of the BHs, in particular if they are of mixed astrophysicaland primordial origin. Using the LVC “truncated” phenomenological model for a single ABH population [44, 70],the mixing fraction parameter inferred from the Bayesian analysis indicates that roughly of the BH mergersin the GWTC-2 catalog are of astrophysical origin. Furthermore, our best-fit mixed ABH+PBH model has adecisive statistical evidence relative to the single-population phenomenological ABH model, strongly supportingthe existence of extra features which are not captured by the simple truncated ABH model. This suggests theexistence of at least two distinct populations of BH mergers in LVC data and is compatible with a subdominantpopulation of PBHs. The fraction of putative PBH events could in particular explain the events that are atodds with the standard astrophysical scenarios, such as those with binary components in the mass gap. Basedon our best-fit mixed model, we predicted that the next-generations ground-based GW interferometers ET andCE could detect up to hundreds of BH merger at redshift z (cid:38) , which can therefore be confidently identifiedas being primordial, provided their merger redshift can be measured with sufficient accuracy. The number of4possible detections at z (cid:38) depends significantly on the low-frequency end of the 3G sensitivity curve. Thus,our analysis suggests that the low-frequency response of 3G detectors will be crucial for an important and uniquescience case [110] such as the unambiguous detection of PBHs.Our findings can be improved along several lines. First, on the PBH side, one might consider a mass functiondifferent from the adopted lognormal and improve on the accretion model. Secondly, on the astrophysical side,one might use more sophisticated and motivated ABH models. The simplest extension of our work in thisrespect would be to include a PBH population to the LVC “truncated+peak” model used for the GWTC-2catalogue [44]. However, being the latter a phenomenological model that attempts to fit the data without anyastrophysical input, it would be hard to disentangle its phenomenological features (in particular its Gaussianpeak) from putative physical effects such as those coming from a subdominant PBH population. A much moreinteresting avenue is to consider astrophysically-motivated ABH populations, for instance by accounting forboth isolated and dynamical formation channels [37, 77], as well as including the merger rates from Pop IIIbinaries, which are relevant for forecasts with 3G detectors [37, 39]. Work along this line is in progress and willbe reported elsewhere [112]. VIII. ACKNOWLEDGMENTS
We thank A. Maselli and R. Schneider for insightful discussions and we are indebted to V. Baibhav, E. Berti, K.K. Y. Ng, S. Vitale, and K. W. K. Wong for discussions and comments on the draft. Computations were per-formed at University of Geneva on the Baobab/Yggdrasil cluster. We acknowledge use of the software package pycbc [113]. V.DL., G.F. and A.R. are supported by the Swiss National Science Foundation (SNSF), project
TheNon-Gaussian Universe and Cosmological Symmetries , project number: 200020-178787. P.P. acknowledges finan-cial support provided under the European Union’s H2020 ERC, Starting Grant agreement no. DarkGRA–757480,and under the MIUR PRIN and FARE programmes (GW-NEXT, CUP: B84I20000100001), and support from theAmaldi Research Center funded by the MIUR program ‘Dipartimento di Eccellenza" (CUP: B81I18001170001).
Appendix A: Binary detectability
A compact quasi-circular binary is characterized by the BHs masses m and m , dimensionless spins χ and χ ,and by the merger redshift z . We have neglected the role of spins in computing the BH binary detectability inthe present work. This is a conservative assumption since spinning binaries typically have larger Signal-to-NoiseRatio (SNR) than nonspinning binaries with same masses and redshift and, in any case, the large majority ofthe GWTC-2 events are compatible with zero spin. Each individual binary is also characterised by the positionand orientation with respect to the detectors. Those are customarily defined in terms of right ascension α ,declination δ , orbital-plane inclination ι , and polarization angle ψ . One can thus designate the intrinsic andextrinsic parameters as θ i = { m , m , z } and θ e = { α, δ, ι, ψ } , respectively.One has to compute the detectability by averaging over the extrinsic parameters θ e . For each value of θ i , onecan define the detection probability as p det ( θ i ) = (cid:90) p ( θ e ) Θ[ ρ ( θ i , θ e ) − ρ thr ] dθ e , (A.1)where p ( θ e ) is the probability distribution function of θ e , Θ indicates the Heaviside step function and ρ is theSNR. Notice that, in the case of the GWTC-1 catalog, it was already shown that p det can be computed in theapproximate single-detector semianalytic framework of Refs. [114, 115] and a SNR threshold ρ thr = 8 withoutencountering significant departures from the large-scale injection campaigns by the LVC for the O1 [116] andO2 [70] runs. We will adopt the same procedure to compute the detectability of binaries also for the O3 run. In this context it is also interesting to note that a very recent analysis argues that the LVC phenomenological model could beentirely explained by dynamically-formed mergers in globular clusters [111]. - - - - - - - - - - FIG. 8:
Left:
Noise curves.
Right:
Horizon redshift for equal mass binaries. Both panels shows the result for the LIGOHanford (H) and Livingston (L) during the O3 run [2] and both the 3G detectors Cosmic Explorer (CE) during phase 1 from [33]and ET at design sensitivity from [32].
One can factor out the dependency of θ e on the SNR to obtain ρ ( θ i , θ e ) = ω ( θ e ) ρ opt ( θ i ) , where ρ opt is theSNR of an “optimal” source located overhead the detector with face-on inclination. Computing the marginalizeddistribution p det ( θ i ) is achieved by evaluating the cumulative distribution function P ( ω thr ) = (cid:82) ω thr p ( ω (cid:48) ) dω (cid:48) at ω thr = ρ thr /ρ opt ( θ i ) .In this work we consider the case of isotropic sources where α, cos δ, cos ι , and ψ are uniformly distributed.Then, for the case of a single detector approximation, non-precessing sources, and considering only the dominantquadrupole moment, the function P ( ω thr ) is found as in Ref. [45].The optimal SNR ρ opt of individual GW events for a source with masses m and m at a redshift z is givenin terms of the GW waveform in fourier space ˜ h ( ν ) by ρ opt ( m , m , z ) ≡ (cid:90) ∞ | ˜ h ( ν ) | S n ( ν ) d ν, (A.2)where the adopted strain noise S n for current LIGO/Virgo observation runs and future 3G detectors (suchas Einstein Telescope and Cosmic Explorer) are shown in Fig 8. In particular, when computing the binarydetectabiliy during the O1-O2 (O3a) observation runs, we adopted the aLIGOEarlyHighSensitivityP1200087 ( aLIGOMidHighSensitivityP1200087 ) noise power spectral densities, as implemented in the publicly availablerepository pycbc [113]. Consistently with the computational framework described to compute the binary de-tectability, we also adopt the non-precessing waveform model IMRPhenomD. Notice finally that the optimalSNR scales like the inverse of the luminosity distance, i.e. ρ ∝ /D L ( z ) .In Fig. 8 we also plot the SNR expected at current and future GW experiment for an optimally orientedequal mass binary with total mass M tot . In the plot we refer to the “Horizon” as the maximum distance atwhich a binary can in principle be observed if optimally oriented with respect to the detectors, while , as the redshift at which those fraction of binaries are observable. Those values corresponds to values of SNR = { , , } , respectively. Appendix B: Accretion onto PBHs in binaries
PBHs form deep in the radiation dominated early universe, being their mass distribution described by Eq. (3.2).Throughout their evolution up to the merger epoch, however, such mass function must be properly evolved [61,117, 120, 121] by taking into account the impact of baryonic accretion.PBHs in binaries may experience sizeable mass changes if they are heavier than O (10) M (cid:12) . The binarysystem as a whole attracts baryonic gas from the surrounding background medium depending on the binary’s6Bondi-Hoyle mass accretion rate ˙ M bin = 4 πλm H n gas v − eff M tot . (B.1)The Bondi-Hoyle rate is computed as a function of the hydrogen mass m H and mean density n gas while alsokeeping track of the binary effective velocity v eff . The additional factor λ is introduced to parametrise correctionsto the mass accretion rate coming from the Hubble expansion, the coupling of the gas to the CMB radiationinduced by Compton scattering and, finally, the gas viscosity. The interested reader can find more details in [117]and references therein. The individual mass accretion rates of the PBHs participating in the binary system aremodulated by the orbital motion and can be parametrised in terms of the mass ratio q as [57, 61] ˙ m = ˙ M bin (cid:112) q ) , ˙ m = ˙ M bin √ q (cid:112) q ) . (B.2)We note here that the accretion formula (B.1) is derived in the Newtonian approximation. It was recentlysuggested in Refs. [118, 119] that there may be an increase of the mass accretion rate due to general-relativisticeffects.It is of crucial importance to account also for the catalysing effect of an additional DM component other thanPBHs [120, 122, 123]. The most up to date constraints on the DM fraction which is allowed to be in the formof PBHs in the relevant mass range highlight the necessity for a secondary DM component which would tend tocluster around the isolated PBHs or PBHs binaries. In order to account for the enhanced gravitational potentialwell catalysing the baryonic mass accretion we included this effect in the parameter λ [61, 120].The accretion rate is expected to decrease sharply with the onset of structure formation and reionization [18,21, 124]. We account for this effect by placing an effective cut-off in redshift z cut-off after which we neglectaccretion effects. There exists however large uncertainties in the details of accretion around the reionizationepoch as, for example, the effect of PBH induced X-Ray pre-heating [125], details of the structure formationand both local and global feedbacks [120, 126] as well as mechanical feedbacks [127] which force us to considerthe cut-off redshift z cut-off as an essentially free parameter in the model.To summarise, the physical effects of mass accretion onto PBHs are expected to be the following:• An overall shift of the mass distribution to larger values as well as an enlargement of the high-mass tail[61];• A redshift dependent modification of the PBH abundance with respect to the DM proportional to f PBH ( z ) /f PBH ( z i ) ∼ (cid:104) m ( z ) (cid:105) / (cid:104) m ( z i ) (cid:105) [128];• PBH binaries tend to become more symmetrical (i.e, q closer to unity), depending on the strength ofaccretion, as the secondary PBH in the binary is expected to consistently inherit a larger relative accretionrate;• mass accretion may lead to an efficient spin-up of both PBHs in the binary provided enough angularmomentum is transferred to the PBH [129–136] (as it is the case for a disk-like accretion geometry). Formore details on the spin accretion in the PBH model we refer to Refs. [57, 61]. [1] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. X , no. 3, 031040 (2019) [astro-ph.HE/1811.12907] .[2] R. Abbott et al. [LIGO Scientific and Virgo], [gr-qc/2010.14527] .[3] L. Barack, V. Cardoso, S. Nissanke, T. P. Sotiriou, A. Askar, C. Belczynski, G. Bertone, E. Bon, D. Blas andR. Brito, et al. Class. Quant. Grav. (2019) no.14, 143001 [gr-qc/1806.05195] .[4] I. Mandel and A. Farmer, [astro-ph.HE/1806.05820] .[5] M. Mapelli, Proc. Int. Sch. Phys. Fermi (2020), 87-121 [astro-ph.HE/1809.09130] . [6] Zel’dovich, Y.B. and Novikov, I.D.: 1967, Soviet Astronomy , 602.[7] S. W. Hawking, Nature , 30-31 (1974)[8] G. F. Chapline, Nature , no.5489, 251-252 (1975)[9] P. Ivanov, P. Naselsky and I. Novikov, Phys. Rev. D , 7173 (1994).[10] S. Blinnikov, A. Dolgov, N. K. Porayko and K. Postnov, JCAP , 036 (2016) [astro-ph.HE/1611.00541] .[11] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. (2018) no.6, 063001 [astro-ph.CO/1801.05235] .[12] A. M. Green and B. J. Kavanagh, [astro-ph.CO/2007.10722] .[13] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, [astro-ph.CO/2002.12778] .[14] S. Bird, I. Cholis, J. B. Munoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli and A. G. Riess,Phys. Rev. Lett. , no. 20, 201301 (2016) [astro-ph.CO/1603.00464] .[15] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Phys. Rev. Lett. , no. 6, 061101 (2016) Erratum: [Phys.Rev. Lett. , no. 5, 059901 (2018)] [astro-ph.CO/1603.08338] .[16] S. Clesse and J. García-Bellido, Phys. Dark Univ. (2017), 142-147 [astro-ph.CO/1603.05234] .[17] S. Wang, Y. F. Wang, Q. G. Huang and T. G. F. Li, Phys. Rev. Lett. (2018) no.19, 191102 [astro-ph.CO/1610.08725] .[18] Y. Ali-Haïmoud, E. D. Kovetz and M. Kamionkowski, Phys. Rev. D , no.12, 123523 (2017) [astro-ph.CO/1709.06576] .[19] M. Raidal, C. Spethmann, V. Vaskonen and H. Veermäe, JCAP (2019), 018 [astro-ph.CO/1812.01930] .[20] N. Fernandez and S. Profumo, JCAP (2019), 022 [astro-ph.HE/1905.13019] .[21] G. Hütsi, M. Raidal and H. Veermäe, Phys. Rev. D , no. 8, 083016 (2019) [astro-ph.CO/1907.06533] .[22] V. Vaskonen and H. Veermäe, Phys. Rev. D (2020) no.4, 043015 [astro-ph.CO/1908.09752] .[23] A. D. Gow, C. T. Byrnes, A. Hall and J. A. Peacock, JCAP (2020), 031 [astro-ph.CO/1911.12685] .[24] A. Hall, A. D. Gow and C. T. Byrnes, Phys. Rev. D (2020), 123524 [astro-ph.CO/2008.13704] .[25] V. De Luca, V. Desjacques, G. Franciolini and A. Riotto, JCAP (2020), 028 [astro-ph.CO/2009.04731] .[26] K. W. K. Wong, G. Franciolini, V. De Luca, V. Baibhav, E. Berti, P. Pani and A. Riotto, [gr-qc/2011.01865] .[27] G. Hütsi, M. Raidal, V. Vaskonen and H. Veermäe, [astro-ph.CO/2012.02786] .[28] R. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. (2020), 101102 [gr-qc/2009.01075] .[29] V. De Luca, V. Desjacques, G. Franciolini, P. Pani and A. Riotto, [astro-ph.CO/2009.01728] .[30] K. Kritos, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, [gr-qc/2012.03585] .[31] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holzand C. Pankow, [astro-ph.HE/2011.10057] .[32] S. Hild, M. Abernathy, F. Acernese, P. Amaro-Seoane, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsugliaand M. Beker, et al. Class. Quant. Grav. (2011), 094013 [gr-qc/1012.0908] .[33] Cosmic Explorer, https://cosmicexplorer.org/researchers.html (2020)[34] R. Schneider, A. Ferrara, B. Ciardi, V. Ferrari and S. Matarrese, Mon. Not. Roy. Astron. Soc. (2000), 385 [astro-ph/9909419] .[35] R. Schneider, A. Ferrara, P. Natarajan and K. Omukai, Astrophys. J. (2002), 30-39 [astro-ph/0111341] .[36] Schneider, R., Ferrara, A., Salvaterra, R., et al. 2003, Nature, 422, 869.[37] K. K. Y. Ng, S. Vitale, W. M. Farr and C. L. Rodriguez, [astro-ph.CO/2012.09876] .[38] B. Liu and V. Bromm, Mon. Not. Roy. Astron. Soc. (2020) no.2, 2475-2495 [astro-ph.CO/2003.00065] .[39] Valiante, R., Colpi, M., Schneider, R., et al. 2021, Mon. Not. R. Astron. Soc, 500, 4095.[40] Thrane, E. and Talbot, C.: 2019, Publications of the Astronomical Society of Australia , e010. [astro-ph.IM/1809.02293] [41] T. J. Loredo, AIP Conf. Proc. (2004) no.1, 195-206 [astro-ph/0409387] .[42] S. R. Taylor and D. Gerosa, Phys. Rev. D (2018) no.8, 083017 [astro-ph.HE/1806.08365] .[43] I. Mandel, W. M. Farr and J. R. Gair, Mon. Not. Roy. Astron. Soc. (2019) no.1, 1086-1093 [physics.data-an/1809.02063] .[44] R. Abbott et al. [LIGO Scientific and Virgo], [astro-ph.HE/2010.14533] .[45] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik and F. Pannarale,Astrophys. J. (2015) no.2, 263 [astro-ph.HE/1405.7016] .[46] Foreman-Mackey, D., Hogg, D.W., Lang, D., and Goodman, J.: 2013, Publications of the Astronomical Society ofthe Pacific , 306. [astro-ph.IM/1202.3665] [47] H. Jeffreys, Harold, The Theory of Probability (3rd ed.). Oxford, England (1998)[48] https://dcc.ligo.org/LIGO-P1800370/public[49] https://dcc.ligo.org/LIGO-P2000223/public [50] C. Pankow, L. Sampson, L. Perri, E. Chase, S. Coughlin, M. Zevin and V. Kalogera, Astrophys. J. (2017) no.2,154 [astro-ph.HE/1610.05633] .[51] S. Vitale, D. Gerosa, C. J. Haster, K. Chatziioannou and A. Zimmerman, Phys. Rev. Lett. (2017) no.25, 251103 [gr-qc/1707.04637] .[52] M. Zevin, C. P. L. Berry, S. Coughlin, K. Chatziioannou and S. Vitale, Astrophys. J. Lett. (2020) no.1, L17 [astro-ph.HE/2006.11293] .[53] S. Bhagwat, V. De Luca, G. Franciolini, P. Pani and A. Riotto, [astro-ph.CO/2008.12320] .[54] S. Bhagwat, V. De Luca, G. Franciolini, P. Pani and A. Riotto, in preparation.[55] Speagle, J.S.: 2020, Monthly Notices of the Royal Astronomical Society , 3132.[56] M. Raidal, C. Spethmann, V. Vaskonen and H. Veermäe, JCAP , 018 (2019) [astro-ph.CO/1812.01930] .[57] V. De Luca, G. Franciolini, P. Pani and A. Riotto, JCAP , 044 (2020) [astro-ph.CO/2005.05641] .[58] A. Dolgov and J. Silk, Phys. Rev. D (1993), 4244-4255[59] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen and H. Veermäe, Phys. Rev. D , no. 2, 023514 (2017) [astro-ph.CO/1705.05567] .[60] A. D. Gow, C. T. Byrnes and A. Hall, [astro-ph.CO/2009.03204] .[61] V. De Luca, G. Franciolini, P. Pani and A. Riotto, JCAP , 052 (2020) [astro-ph.CO/2003.02778] .[62] V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra and A. Riotto, JCAP (2019), 018 [astro-ph.CO/1903.01179] .[63] M. Mirbabayi, A. Gruzinov and J. Norena, JCAP (2020), 017 [astro-ph.CO/1901.05963] .[64] V. Vaskonen and H. Veermäe, Phys. Rev. D (2020) no.4, 043015 [astro-ph.CO/1908.09752] .[65] K. Jedamzik, JCAP (2020), 022 [astro-ph.CO/2006.11172] .[66] S. Young and A. S. Hamers, JCAP (2020), 036 [astro-ph.CO/2006.15023] .[67] K. Jedamzik, [astro-ph.CO/2007.03565] .[68] M. Tkachev, S. Pilipenko and G. Yepes, Mon. Not. Roy. Astron. Soc. (2020) no.4, 4854-4862 [astro-ph.CO/2009.07813] .[69] D. Inman and Y. Ali-Haïmoud, Phys. Rev. D , no.8, 083528 (2019) [astro-ph.CO/1907.08129] .[70] B. P. Abbott et al. [LIGO Scientific and Virgo], Astrophys. J. Lett. (2019) no.2, L24 [astro-ph.HE/1811.12940] .[71] C. Cutler and E. E. Flanagan, Phys. Rev. D (1994), 2658-2697 [gr-qc/9402014] .[72] E. Poisson and C. M. Will, Phys. Rev. D (1995), 848-855 [gr-qc/9502040] .[73] E. Baird, S. Fairhurst, M. Hannam and P. Murphy, Phys. Rev. D (2013) no.2, 024035 [gr-qc/1211.0546] .[74] M. Pürrer, M. Hannam, P. Ajith and S. Husa, Phys. Rev. D (2013), 064007 [gr-qc/1306.2320] .[75] K. Chatziioannou, N. Cornish, A. Klein and N. Yunes, Phys. Rev. D (2014) no.10, 104023 [gr-qc/1404.3180] .[76] K. K. Y. Ng, S. Vitale, A. Zimmerman, K. Chatziioannou, D. Gerosa and C. J. Haster, Phys. Rev. D (2018)no.8, 083007 [gr-qc/1805.03046] .[77] K. W. K. Wong, K. Breivik, K. Kremer and T. Callister, [astro-ph.HE/2011.03564] .[78] G. Rakavy, G. Shaviv, Astrophys. J. (1967), 803[79] Z. Barkat, G. Rakavy and N. Sack, Phys. Rev. Lett. (1967), 379-381[80] Gary S. Fraley, Astrophys. & Space Sciences (1968) no.1, 96 – 114[81] A. Heger and S. E. Woosley, Astrophys. J. (2002), 532-543 [astro-ph/0107037] .[82] S. E. Woosley, S. Blinnikov and A. Heger, Nature (2007), 390 [astro-ph/0710.3314] .[83] K. Belczynski, A. Heger, W. Gladysz, A. J. Ruiter, S. Woosley, G. Wiktorowicz, H. Y. Chen, T. Bu-lik, R. O’Shaughnessy, D. E. Holz, C. L. Fryer and E. Berti, Astron. Astrophys. (2016), A97 [astro-ph.HE/1607.03116] .[84] S. Woosley, Astrophys. J. (2017) no.2, 244 [astro-ph.HE/1608.08939] .[85] S. Stevenson, M. Sampson, J. Powell, A. Vigna-Gómez, C. J. Neijssel, D. Szécsi and I. Mandel, [astro-ph.HE/1904.02821] .[86] R. Farmer, M. Renzo, S. E. de Mink, P. Marchant and S. Justham, [astro-ph.SR/1910.12874] .[87] M. Mapelli, M. Spera, E. Montanari, M. Limongi, A. Chieffi, N. Giacobbo, A. Bressan and Y. Bouffanais, Astrophys.J. , 76 (2020) [astro-ph.HE/1909.01371] .[88] M. Renzo, R. J. Farmer, S. Justham, S. E. de Mink, Y. Götberg and P. Marchant, Mon. Not. Roy. Astron. Soc. (2020) no.3, 4333-4341 [astro-ph.SR/2002.08200] .[89] P. Marchant and T. Moriya, Astron. Astrophys. , L18 (2020) [astro-ph.HE/2007.06220] .[90] D. Croon, S. D. McDermott and J. Sakstein, Phys. Rev. D , 115024 (2020) [gr-qc/2007.07889] .[91] J. Ziegler and K. Freese, [astro-ph.HE/2010.00254] .[92] M. Fishbach, D. E. Holz and B. Farr, Astrophys. J. Lett. (2017) no.2, L24 [astro-ph.HE/1703.06869] . [93] D. Gerosa and E. Berti, Phys. Rev. D (2019) no.4, 041301 [astro-ph.HE/1906.05295] .[94] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio and C. S. Ye, Phys. Rev. D (2019) no.4, 043027 [astro-ph.HE/1906.10260] .[95] V. Baibhav, D. Gerosa, E. Berti, K. W. K. Wong, T. Helfer and M. Mould, Phys. Rev. D (2020) no.4, 043002 [astro-ph.HE/2004.00650] .[96] C. Kimball, C. Talbot, C. P. L. Berry, M. Carney, M. Zevin, E. Thrane and V. Kalogera, [astro-ph.HE/2005.00023] .[97] J. Samsing and K. Hotokezaka, [astro-ph.HE/2006.09744] .[98] M. Mapelli, F. Santoliquido, Y. Bouffanais, M. A. Sedda, N. Giacobbo, M. C. Artale and A. Ballone, [astro-ph.HE/2007.15022] .[99] D. Gerosa and E. Berti, Phys. Rev. D , no.12, 124046 (2017) [gr-qc/1703.06223] .[100] M. Fishbach, D. E. Holz and W. M. Farr, Astrophys. J. Lett. (2018) no.2, L41 [astro-ph.HE/1805.10270] .[101] T. Callister, M. Fishbach, D. Holz and W. Farr, Astrophys. J. Lett. (2020) no.2, L32 [astro-ph.HE/2003.12152] .[102] M. Fishbach, Z. Doctor, T. Callister, B. Edelman, J. Ye, R. Essick, W. M. Farr, B. Farr and D. E. Holz, [astro-ph.HE/2101.07699] .[103] S. Vitale, W. M. Farr, K. Ng and C. L. Rodriguez, Astrophys. J. Lett. (2019) no.1, L1 [astro-ph.HE/1808.00901] .[104] Z. C. Chen and Q. G. Huang, JCAP (2020), 039 [astro-ph.CO/1904.02396] .[105] X. Ding, K. Liao, M. Biesiada and Z. H. Zhu, [astro-ph.GA/2002.02981] .[106] T. Kinugawa, K. Inayoshi, K. Hotokezaka, D. Nakauchi and T. Nakamura, Mon. Not. Roy. Astron. Soc. (2014)no.4, 2963-2992 [astro-ph.HE/1402.6672] .[107] T. Kinugawa, A. Miyamoto, N. Kanda and T. Nakamura, Mon. Not. Roy. Astron. Soc. (2016) no.1, 1093-1114 [astro-ph.SR/1505.06962] .[108] T. Hartwig, M. Volonteri, V. Bromm, R. S. Klessen, E. Barausse, M. Magg and A. Stacy, Mon. Not. Roy. Astron.Soc. (2016) no.1, L74-L78 [astro-ph.GA/1603.05655] .[109] S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. (2017) no.22, 221104 [astro-ph.CO/1708.07380] .[110] M. Maggiore, C. Van Den Broeck, N. Bartolo, E. Belgacem, D. Bertacca, M. A. Bizouard, M. Branchesi, S. Clesse,S. Foffa and J. García-Bellido, et al. JCAP (2020), 050 [astro-ph.CO/1912.02622] .[111] C. L. Rodriguez, K. Kremer, S. Chatterjee, G. Fragione, A. Loeb, F. A. Rasio, N. C. Weatherford and C. S. Ye, [astro-ph.HE/2101.07793] .[112] V. Baibhav, E. Berti, V. De Luca, G. Franciolini, K. K. Y. Ng, P. Pani, A. Riotto, S. Vitale and K. W. K. Wong,in preparation.[113] A. Nitz et al., “gwastro/pycbc: Pycbc release v1.15.5", (2020).[114] L. S. Finn and D. F. Chernoff, Phys. Rev. D (1993), 2198-2219 [gr-qc/9301003] .[115] L. S. Finn, Phys. Rev. D (1996), 2878-2894 [gr-qc/9601048] .[116] B. P. Abbott et al. [LIGO Scientific and Virgo], Astrophys. J. Lett. (2016) no.1, L1 [astro-ph.HE/1602.03842] .[117] M. Ricotti, Astrophys. J. , 53 (2007) [astro-ph/0706.0864] .[118] A. Cruz-Osorio and L. Rezzolla, Astrophys. J. (2020) no.2, 147 [gr-qc/2004.13782] .[119] A. Cruz-Osorio, F. D. Lora-Clavijo and C. Herdeiro, [astro-ph.HE/2101.01705] .[120] M. Ricotti, J. P. Ostriker and K. J. Mack, Astrophys. J. , 829 (2008) [astro-ph/0709.0524] .[121] J. R. Rice and B. Zhang, JHEAp , 22 (2017) [astro-ph.HE/1702.08069] .[122] J. Adamek, C. T. Byrnes, M. Gosenca and S. Hotchkiss, Phys. Rev. D (2019) no.2, 023506 [astro-ph.CO/1901.08528] .[123] K. J. Mack, J. P. Ostriker and M. Ricotti, Astrophys. J. , 1277 (2007) [astro-ph/0608642] .[124] G. Hasinger, JCAP (2020), 022 [astro-ph.CO/2003.05150] .[125] S. Oh and Z. Haiman, Mon. Not. Roy. Astron. Soc. (2003), 456 [astro-ph/0307135] .[126] Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D , no. 4, 043534 (2017) [astro-ph.CO/1612.05644] .[127] V. Bosch-Ramon and N. Bellomo, Astron. Astrophys. , A132 (2020) [astro-ph.CO/2004.11224] .[128] V. De Luca, G. Franciolini, P. Pani and A. Riotto, [astro-ph.CO/2003.12589] .[129] E. Berti and M. Volonteri, Astrophys. J. , 822 (2008) [astro-ph/0802.0025] .[130] N. Shakura and R. Sunyaev, Astron. Astrophys. (1973), 337-355[131] Novikov, I.D. and Thorne, K.S., Astrophysics and black holes, 343–550[132] J. M. Bardeen, W. H. Press and S. A. Teukolsky, Astrophys. J. , 347 (1972).[133] R. Brito, V. Cardoso and P. Pani, Class. Quant. Grav. , no. 13, 134001 (2015) [gr-qc/1411.0686] .[134] M. Volonteri, P. Madau, E. Quataert and M. J. Rees, Astrophys. J. , 69 (2005) [astro-ph/0410342] . [135] K. S. Thorne, Astrophys. J. (1974), 507-520[136] C. F. Gammie, S. L. Shapiro and J. C. McKinney, Astrophys. J. , 312 (2004) [astro-ph/0310886][astro-ph/0310886]