New insights on binary black hole formation channels after GWTC-2: young star clusters versus isolated binaries
Yann Bouffanais, Michela Mapelli, Filippo Santoliquido, Nicola Giacobbo, Ugo N. Di Carlo, Sara Rastello, M. Celeste Artale, Giuliano Iorio
MMNRAS , 000–000 (0000) Preprint 26 February 2021 Compiled using MNRAS L A TEX style file v3.0
New insights on binary black hole formation channels afterGWTC-2: young star clusters versus isolated binaries
Yann Bouffanais , (cid:63) , Michela Mapelli , , † , Filippo Santoliquido , , Nicola Giacobbo , , , ,Ugo N. Di Carlo , , , , Sara Rastello , , M. Celeste Artale , Giuliano Iorio , Physics and Astronomy Department Galileo Galilei, University of Padova, Vicolo dell’Osservatorio 3, I–35122, Padova, Italy INFN-Padova, Via Marzolo 8, I–35131 Padova, Italy INAF-Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I–35122, Padova, Italy Dipartimento di Scienza e Alta Tecnologia, University of Insubria, Via Valleggio 11, I–22100, Como, Italy Institut f¨ur Astro- und Teilchenphysik, Universit¨at Innsbruck, Technikerstrasse 25/8, A-6020, Innsbruck, ¨Osterreich School of Physics and Astronomy & Institute for Gravitational Wave Astronomy, University of Birmingham, Birmingham, B15 2TT, UK
26 February 2021
ABSTRACT
With the recent release of the second gravitational-wave transient catalogue (GWTC-2), which introduced dozensof new detections, we are at a turning point of gravitational wave astronomy, as we are now able to directly inferconstraints on the astrophysical population of compact objects. Here, we tackle the burning issue of understanding theorigin of binary black hole (BBH) mergers. To this effect, we make use of state-of-the art population synthesis andN-body simulations, to represent two distinct formation channels: BBHs formed in the field (isolated channel) andin young star clusters (dynamical channel). We then use a Bayesian hierarchical approach to infer the distributionof the mixing fraction f , with f = 0 ( f = 1) in the pure dynamical (isolated) channel. We explore the effects ofadditional hyper-parameters of the model, such as the spread in metallicity σ Z and the parameter σ sp , describingthe distribution of spin magnitudes. We find that the dynamical model is slightly favoured with a median value of f = 0 .
26, when σ sp = 0 . σ Z = 0 .
4. Models with higher spin magnitudes tend to strongly favour dynamicallyformed BBHs ( f (cid:54) . σ sp = 0 . σ Z , have a large impact on the inference of the mixing fraction, which rises from 0 .
18 to 0 .
43 when we increase σ Z from 0.2 to 0.6, for a fixed value of σ sp = 0 .
1. Finally, our current set of observations is better described by acombination of both formation channels, as a pure dynamical scenario is excluded at the 99% credible interval, exceptwhen the spin magnitude is high.
Key words: black hole physics – gravitational waves – methods: numerical – methods: statistical analysis
In 2015, the LIGO–Virgo collaboration (LVC) announced thefirst direct detection of gravitational waves (GWs) emitted bythe merger of a binary black hole (BBH, Abbott et al. 2016b;Abbott et al. 2016c). Nowadays, the LVC has published atotal of 50 binary compact object mergers, mostly BBHs, asa result of the first (O1, Abbott et al. 2016a), the second (O2,Abbott et al. 2019a,b) and the first half of the third observingrun (O3a, Abbott et al. 2020a,b). In addition, Zackay et al.(2019), Udall et al. (2019), Venumadhav et al. (2020) and Nitzet al. (2020) claim several additional BBH candidates, basedon an independent analysis of the O1 and O2 LVC data.This extraordinary wealth of data gives an unprecedentedopportunity to understand the formation channels of BBHs (cid:63)
E-mail:bouff[email protected], yann.bouff[email protected] † E-mail:[email protected] (see, e.g., Mandel & Farmer 2018 and Mapelli 2018 for areview). In the isolated formation channel, BBHs originatefrom the evolution of massive binary stars. A BBH that formsin isolation can reach coalescence within a Hubble time onlyif its stellar progenitors evolve via common envelope (e.g.,Tutukov & Yungelson 1973; Bethe & Brown 1998; PortegiesZwart & Yungelson 1998; Belczynski et al. 2002, 2008, 2016;Eldridge & Stanway 2016; Stevenson et al. 2017a; Mapelli et al.2017, 2019; Kruckow et al. 2018; Spera et al. 2019; Tanikawaet al. 2020; Belczynski et al. 2020), stable mass transfer (e.g.,Giacobbo et al. 2018; Neijssel et al. 2019; Bavera et al. 2020;Bouffanais et al. 2020; Andrews et al. 2020), or chemicallyhomogeneous evolution (e.g., Marchant et al. 2016; Mandel& de Mink 2016; de Mink & Mandel 2016; du Buisson et al.2020). These processes tend to pose a limit of ∼ −
100 M (cid:12) to the total maximum mass of a BBH merger (e.g., Bouffanaiset al. 2019) and to align the spins of the final black holes (BHs)with the orbital angular momentum of the binary system (e.g., © a r X i v : . [ a s t r o - ph . H E ] F e b Bouffanais et al.
Mandel & de Mink 2016; Rodriguez et al. 2016b; Gerosa et al.2018). Only supernova explosions can partially misalign thesystems (e.g., Kalogera 2000).Alternatively, BBHs can form dynamically: BHs can pair upwith other BHs in dense stellar systems, such as nuclear starclusters (e.g., Antonini & Rasio 2016; Petrovich & Antonini2017; Arca Sedda 2020; Fragione & Silk 2020), globular clus-ters (e.g., Portegies Zwart & McMillan 2000; Tanikawa 2013;Rodriguez et al. 2016a; Askar et al. 2017; Fragione & Kocsis2018; Choksi et al. 2018, 2019; Hong et al. 2018; Rodriguez& Loeb 2018), and young star clusters (e.g., Banerjee et al.2010; Ziosi et al. 2014; Mapelli 2016; Banerjee 2017, 2020;Di Carlo et al. 2019, 2020b; Kumamoto et al. 2019, 2020).Finally, gas torques in AGN discs (e.g., Bartos et al. 2017;Stone et al. 2017; McKernan et al. 2018; Yang et al. 2019;Tagawa et al. 2019) and hierarchical triples (e.g., Antoniniet al. 2017; Silsbee & Tremaine 2017; Fragione & Kocsis 2020;Vigna-G´omez et al. 2021) can also facilitate the formation ofBBH mergers. Unlike isolated BBHs, dynamical BBHs canhave total masses in excess of ∼
100 M (cid:12) as a result of (run-away) stellar mergers in young star clusters (e.g., PortegiesZwart et al. 2004; Mapelli 2016; Di Carlo et al. 2020a; Rizzutoet al. 2021) or hierarchical BBH mergers (e.g., Miller & Hamil-ton 2002; Giersz et al. 2015; Fishbach et al. 2017; Gerosa &Berti 2017; Rodriguez et al. 2019). Moreover, dynamics resetthe memory of spin orientation: we expect that dynamicalBBHs have isotropically oriented spins (e.g., Rodriguez et al.2016b).Based on these two main differences on the mass and spindistribution of dynamical versus isolated BBHs, several studieshave explored the possibility of using LVC data to constrainthe formation channels of BBHs (e.g., Mandel & de Mink2016; Zevin et al. 2017; Stevenson et al. 2017b; Bouffanaiset al. 2019; Callister et al. 2020; Zevin et al. 2020; Wonget al. 2020). Considering O1, O2 and O3a events, the LVCcollaboration recently showed that 12–44% of BBHs havespins tilted by more than 90 ◦ with respect to their orbitalangular momentum (Abbott et al. 2020b). This suggests thatisolated BBHs can hardly account for the entire sample. Here,we compare the properties of BBHs in the second GW tran-sient catalogue (hereafter, GWTC-2, Abbott et al. 2020a,b)against our population synthesis and dynamical simulations.We consider isolated binary evolution and dynamical forma-tion in young dense star clusters (Santoliquido et al. 2020).We show that GWTC-2 data strongly disfavour the possibilitythat all observed BBHs come from isolated binary evolution,even when we take into account the main uncertainties onmetallicity evolution and spin magnitudes.This paper is organised as follows. Section 2 describes ourastrophysical models. In Section 3, we summarize our Bayesianhierarchical analysis. We present and discuss our results inSections 4 and 5, respectively. A short summary is providedin Section 6. The isolated BBHs were simulated with the population synthe-sis code mobse (Giacobbo et al. 2018; Mapelli & Giacobbo2018). With respect to its progenitor bse (Hurley et al. 2000,2002), mobse contains updated prescriptions for stellar winds.In particular, mass loss by massive hot stars is expressedas ˙ M ∝ Z γ , where γ is a function of the stellar luminositythrough the Eddington ratio (Giacobbo & Mapelli 2018). In mobse , the mass of a compact object depends on the finaltotal mass and on the final carbon-oxygen core mass of itsprogenitor star, as described in Fryer et al. (2012). Here, weadopt the rapid core-collapse supernova model from Fryeret al. (2012), which enforces a compact-object mass gap be-tween 2 and 5 M (cid:12) . Electron-capture supernovae are describedas in Giacobbo & Mapelli (2019). Pulsational pair instabilityand pair-instability supernovae are modelled as in Mapelliet al. (2020). These models yield a maximum BH mass of ≈
65 M (cid:12) .We assign natal kicks to the BHs as v BH = v Max (1 − f fb ), where f fb is the fraction of fallback as defined in Fryeret al. (2012), while v Max is a random number drawn froma Maxwellian distribution with one-dimensional root-meansquare velocity σ (Hobbs et al. 2005). In this manuscript, weadopt σ = 15 km s − . The main binary evolution processes(mass transfer, common envelope and tides) are described asin Hurley et al. (2002). For common envelope, we adopt the α formalism with α = 5, while the parameter λ is calculatedself-consistently with the prescriptions derived by Claeys et al.(2014). Orbital decay by GW emission is implemented as inPeters (1964). All BBHs evolve by GW emission, regardlessof their orbital separation.BH spin magnitudes are drawn from a Maxwellian distri-bution with σ sp = 0 .
1, in the fiducial case. We also discussthe two extreme cases with σ sp = 0 .
01 and 0.3. Theoreticalmodels of BH spin magnitudes are affected by substantialuncertainties, mostly because of our poor understanding ofangular momentum transfer in the stellar interior (e.g., Bel-czynski et al. 2020). The case with σ sp = 0 .
01 corresponds toassuming extremely efficient angular momentum dissipation(Fuller & Ma 2019), while values of σ sp = 0 . − . θ = cos ( ν ) cos ( ν ) + sin ( ν ) sin ( ν ) cos ( φ ) , (1)where ν is the angle between the orbital angular momentumvector after ( (cid:126)L new ) and before a SN explosion ( (cid:126)L old ), so thatcos ( ν ) = (cid:126)L new L new · (cid:126)L old L old . (2)In eq. 1, ν ( ν ) corresponds to the tilt induced by the first https://mobse-webpage.netlify.app/ MNRAS000
01 corresponds toassuming extremely efficient angular momentum dissipation(Fuller & Ma 2019), while values of σ sp = 0 . − . θ = cos ( ν ) cos ( ν ) + sin ( ν ) sin ( ν ) cos ( φ ) , (1)where ν is the angle between the orbital angular momentumvector after ( (cid:126)L new ) and before a SN explosion ( (cid:126)L old ), so thatcos ( ν ) = (cid:126)L new L new · (cid:126)L old L old . (2)In eq. 1, ν ( ν ) corresponds to the tilt induced by the first https://mobse-webpage.netlify.app/ MNRAS000 , 000–000 (0000) ew insights on BBH formation (second) supernova, while φ is the phase of the projectionof (cid:126)L into the orbital plane. Our formalism neglects possiblere-alignments of the spin of the first born BH before theformation of the second BH. We get ν and ν directly from mobse , while we have to generate φ as a uniform randomnumber between 0 and 2 π (Gerosa et al. 2013; Rodriguezet al. 2016b).We have simulated binary stars with 12 different metallici-ties: Z = 0 . bina-ries per each metallicity comprised between Z = 0 . × binaries per each metallicity Z (cid:62) . . × isolated bi-naries. The zero-age main-sequence masses of the primarycomponent of each binary star are distributed according toa Kroupa (Kroupa 2001) initial mass function in the range[5 , (cid:12) . The orbital periods, eccentricities and mass ra-tios of binaries are drawn from Sana et al. (2012). In par-ticular, we derive the mass ratio q as F ( q ) ∝ q − . with q ∈ [0 . − P from F (Π) ∝ Π − . withΠ = log ( P/ day) ∈ [0 . − .
5] and the eccentricity e from F ( e ) ∝ e − . with 0 (cid:54) e (cid:54) . The dynamical simulations have already been presented inprevious work (Di Carlo et al. 2020b; Rastello et al. 2020;Santoliquido et al. 2020). Here, we summarize their mainfeatures, while we refer to the aforementioned papers formore details. We simulated dense young star clusters withtotal masses M SC between 300 and 3 × M (cid:12) , randomlydrawn from a distribution d N/ d M SC ∝ M − , consistent withobservations (Lada & Lada 2003; Portegies Zwart et al. 2010).In particular, we analyze 100002 simulations of star clusterswith mass 300 − (cid:12) from Rastello et al. (2020) and6000 simulations of star clusters with mass 1000 − (cid:12) , corresponding to the union of set A and B presentedin Di Carlo et al. (2020b). The dynamical simulations areequally divided between three metallicities: Z = 0 . , nbody6++gpu (Wang et al. 2015, 2016),coupled with our population synthesis code mobse (Giacobboet al. 2018). This guarantees that the treatment of stellarand binary evolution is the same for isolated and dynamicalBBHs.The global initial binary fraction of the dynamical simu-lations is f bin = 0 .
4, but the binary fraction is assumed tocorrelate with mass (K¨upper et al. 2011). As a result, all starswith masses > (cid:12) are initially members of a binary system(Di Carlo et al. 2019) , consistent with observations (Sanaet al. 2012; Moe & Di Stefano 2017). The orbital periods, ec-centricities and mass ratios of the original binaries are drawnfrom the same distribution as the isolated binaries.Spin magnitudes of dynamical BBHs are drawn from aMaxwellian distribution in the same way as we did for isolatedBBHs. In particular, we adopt σ sp = 0 . σ sp = 0 .
01 and 0.3. Dynamicalprocesses such as exchanges, flybys and captures tend tomisalign the spins. Hence, we assume that all components of the dynamically formed BBHs have isotropic spins over thesphere.
We calculated the merger rate density as R ( z ) = dd t lb ( z ) (cid:90) zz max (cid:20)(cid:90) Z max Z min η ( Z ) F ( z (cid:48) , z, Z ) d Z (cid:21) (3) × ψ ( z (cid:48) ) d t lb ( z (cid:48) )d z (cid:48) d z (cid:48) , where t lb ( z ) is the look-back time at redshift z , Z min and Z max are the minimum and maximum metallicity, ψ ( z (cid:48) ) isthe cosmic star formation rate (SFR) density at redshift z (cid:48) .In eq. 3, we used the fit to the SFR density from Madau &Fragos (2017): ψ ( z ) = 0 .
01 (1 + z ) . z ) / . . M (cid:12) Mpc − yr − . (4)In eq. 3, η ( Z ) is the merger efficiency, namely the ratio be-tween the total number N TOT ( Z ) of compact binaries (formedfrom a coeval population) that merge within an Hubble time( t H (cid:46)
14 Gyr) and the total initial mass M ∗ ( Z ) of the simu-lation with metallicity Z .Finally, F ( z (cid:48) , z, Z ) is the fraction of compact binaries thatform at redshift z (cid:48) from stars with metallicity Z and mergeat redshift z : F ( z (cid:48) , z, Z ) = N ( z (cid:48) , z, Z ) N TOT ( Z ) p ( z (cid:48) , Z ) , (5)where p ( z (cid:48) , Z ) is the distribution of stellar metallicities atredshift z (cid:48) and N ( z (cid:48) , z, Z ) is the total number of compactbinaries that merge at redshift z and form from stars withmetallicity Z at redshift z (cid:48) .The merger rate density is dramatically affected by themetallicity evolution of stars across cosmic time, which isneeded to calculate the term p ( z (cid:48) , Z ) of eq. 5. Here we usethe fit to the average metallicity evolution given by Madau &Fragos (2017):log (cid:104) Z/Z (cid:12) (cid:105) = 0 . − . z . . (6)To describe the spread around the average metallicity, weassume that metallicities are distributed according to a log-normal distribution with average log (cid:104) Z/Z (cid:12) (cid:105) and standarddeviation σ Z : p ( z (cid:48) , Z ) = 1 (cid:112) π σ Z exp (cid:40) − [log ( Z ( z (cid:48) ) / Z (cid:12) ) − log (cid:104) Z ( z (cid:48) ) /Z (cid:12) (cid:105) ] σ Z (cid:41) . (7)The standard deviation σ Z is highly uncertain. Here, weprobe different values of σ Z = 0 . , cosmo R ate . We refer to Santoliquido et al. (2020) and San-toliquido et al. (2021) for further details. Our astrophysical models can be described as a functionof a set of hyper-parameters λ . In this study, the hyper-parameters λ are a combination of the formation channeltype (either isolated or dynamical), the metallicity dispersion MNRAS , 000–000 (0000)
Bouffanais et al. ( σ Z = { . , . , . } ) and the spin magnitude root-meansquare ( σ sp = { . , . , . } ).A given model will have a prediction on the distribution ofmerging BBHs, d N d θ ( λ ) = N λ p ( θ | λ ) , (8)where θ are the parameters of BBH mergers and N λ is thetotal number of mergers predicted by the model computed as N λ = (cid:90) z h R ( z ) d V c d z T obs z d z, (9)where d V c / d z is the comoving volume element and T obs is theobservation time considered in the analysis. Note that theintegral in eq. 9 is done over redshift ranging from 0 to z h ,where z h is the horizon redshift, i.e. the redshift correspondingto the instrumental horizon of LIGO and Virgo. For thisspecific analysis, we selected a value z h = 2, as it sets a safeboundary for the detection limit of the detectors during thefirst three observing runs.To model the population of merging BBHs, we have chosenthe parameterisation θ = {M c , q, z, χ eff } where M c is thechirp mass and χ eff the effective spin, χ eff = ( m (cid:126)χ + m (cid:126)χ ) m + m · (cid:126)LL , (10)where m ( m ) is the primary (secondary) BH mass, χ ( χ )is the primary (secondary) dimensionless spin magnitude and (cid:126)L is the orbital angular momentum of the BBH.To compute the distribution p ( θ | λ ), we first constructedcatalogues of 3 × sources for all possible combinations ofhyper-parameters λ , using the merger rate and metallicitydistribution calculated by cosmo R ate . From these catalogues,we can derive a continuous estimation of p ( θ | λ ) by makinguse of Gaussian kernel density estimation. A value of thebandwidth of 0 .
075 proved to be our optimal choice to describeour distributions.Figure 1 shows the sources in our catalogues, for both for-mation channels, assuming σ Z = { . , . } and σ sp = 0 . ≈
60 M (cid:12) ,while M c caps at ≈
35 M (cid:12) for the isolated model. Moreover,the distribution of mass ratios is more peaked towards q ≈ σ Z on masses’ distribution is relatively small forboth formation channels.Both the metallicity spread σ Z and the type of formationchannels do have an impact on the shape of the redshiftdistribution. In particular, for a given formation channel, thepeak of the redshift distribution is shifted towards lower valuesof z as the value of σ Z increases. This happens because theBBH merger rate strongly depends on progenitor’s metallicity:BBH mergers are ∼ − ∼ −
2) orders of magnitude moreefficient at low metallicity than at high metallicity for isolated(dynamical) binaries. Hence, a larger metallicity spread, whichmeans a larger fraction of metal-poor stars at low redshift,implies more mergers at low redshift (Santoliquido et al. 2021).In addition, for a fixed value of σ Z , the distribution of redshiftcarries information on the formation channel. The peak of the distribution moves depending on formation channels, withfor instance a peak close to z = 1 . z = 1 . σ Z = 0 .
6. Anotherimportant feature is represented by the different values ofcurvature and slope at low redshift ( z ∈ [0 , χ eff depending on formation channels. In particular, thedynamical model has equal support for positive and negativevalues of χ eff , and has a spread close to σ sp ; the isolatedmodel has a very strong support towards positive valuesand is centered around 0 .
15. By construction, our isolatedmodel allows for sources with values of χ eff as low as − . χ eff (cid:54) { . , . , . } for σ Z = { . , . , . } , respectively.These sources do nevertheless have an importance from theanalysis point of view, since they allow the isolated model tohave some support for negative values of χ eff . Hierarchical Bayesian inference has proved to be an invalu-able tool to estimate and constrain features of a populationof merging compact objects. The approach we used has beendescribed in previous studies (e.g., Mandel et al. 2019; Bouf-fanais et al. 2019), so we only present here the main equations.Given an ensemble of N obs observations, { h } k , the poste-rior distribution of the hyper-parameter λ is described as aninhomogeneous Poisson distribution p ( λ, N λ |{ h } k ) ∼ e − µ ( λ ) π ( λ, N λ ) N obs (cid:89) k =1 N λ (cid:90) θ L k ( h k | θ ) p ( θ | λ ) d θ, (11)where π ( λ, N λ ) is the prior distribution on λ and N λ , µ ( λ ) isthe predicted number of detections for the model and L k ( h k | θ )is the likelihood of the k th detection. The predicted numberof detections is given by µ ( λ ) = N λ × β ( λ ) , (12)where β ( λ ) is the detection efficiency of the model with hyper-parameter λ , that can be computed as β ( λ ) = (cid:90) p ( θ | λ ) p det ( θ ) d θ, (13)where p det ( θ ) is the probability of detecting a source with pa-rameters θ . This probability can be inferred by computing theoptimal signal-to-noise ratio (SNR) of the source and compar-ing it to a detection threshold. In our case, we computed theoptimal SNR using LIGO Livingston as a reference, for whichwe approximated the sensitivity using the averaged sensitivityover all detections for all three observing runs separately. Fur-thermore, by putting the detection threshold at ρ thr = 8, itwas shown that this single-detector approximation is a goodrepresentation of more complex analysis with a network ofdetectors (Abadie et al. 2010; Abbott et al. 2016d; Wysockiet al. 2018).The values for the event’s log-likelihood were derived fromthe posterior and prior samples released by the LVC, such MNRAS000
15. By construction, our isolatedmodel allows for sources with values of χ eff as low as − . χ eff (cid:54) { . , . , . } for σ Z = { . , . , . } , respectively.These sources do nevertheless have an importance from theanalysis point of view, since they allow the isolated model tohave some support for negative values of χ eff . Hierarchical Bayesian inference has proved to be an invalu-able tool to estimate and constrain features of a populationof merging compact objects. The approach we used has beendescribed in previous studies (e.g., Mandel et al. 2019; Bouf-fanais et al. 2019), so we only present here the main equations.Given an ensemble of N obs observations, { h } k , the poste-rior distribution of the hyper-parameter λ is described as aninhomogeneous Poisson distribution p ( λ, N λ |{ h } k ) ∼ e − µ ( λ ) π ( λ, N λ ) N obs (cid:89) k =1 N λ (cid:90) θ L k ( h k | θ ) p ( θ | λ ) d θ, (11)where π ( λ, N λ ) is the prior distribution on λ and N λ , µ ( λ ) isthe predicted number of detections for the model and L k ( h k | θ )is the likelihood of the k th detection. The predicted numberof detections is given by µ ( λ ) = N λ × β ( λ ) , (12)where β ( λ ) is the detection efficiency of the model with hyper-parameter λ , that can be computed as β ( λ ) = (cid:90) p ( θ | λ ) p det ( θ ) d θ, (13)where p det ( θ ) is the probability of detecting a source with pa-rameters θ . This probability can be inferred by computing theoptimal signal-to-noise ratio (SNR) of the source and compar-ing it to a detection threshold. In our case, we computed theoptimal SNR using LIGO Livingston as a reference, for whichwe approximated the sensitivity using the averaged sensitivityover all detections for all three observing runs separately. Fur-thermore, by putting the detection threshold at ρ thr = 8, itwas shown that this single-detector approximation is a goodrepresentation of more complex analysis with a network ofdetectors (Abadie et al. 2010; Abbott et al. 2016d; Wysockiet al. 2018).The values for the event’s log-likelihood were derived fromthe posterior and prior samples released by the LVC, such MNRAS000 , 000–000 (0000) ew insights on BBH formation q p ( q ) Isolated Z = 0.2Isolated Z = 0.6Dynamical Z = 0.2Dynamical Z = 0.60.0 0.5 1.0 1.5 2.0 z p ( z ) eff p ( e ff )
10 20 30 40 50 60 c /M p ( c ) Figure 1.
Distribution of the model’s parameters as inferred from a catalogue of N tot = 3 × sources. From left to right and from topto bottom, we plot the distribution of the chirp mass M c , mass ratio q , redshift z and effective spin χ eff for the isolated and dynamicalformation channels. The spread in metallicity is equal to σ Z = 0 . σ Z = 0 .
6) for the solid (dashed) lines. The root-mean square of theMaxwellian for the spins is σ sp = 0 . that the integral in eq. 11 is approximated with a Monte Carloapproach as (cid:90) L k ( h k | θ ) p ( θ | λ ) d θ ≈ N ks N ks (cid:88) i =1 p ( θ ki | λ ) π k ( θ ki ) , (14)where θ ki is the i th posterior sample for the k th detectionand N ks is the total number of posterior samples for the k th detection. To compute the prior term in the denominator, wealso used Gaussian kernel density estimation.Finally, we can also choose to neglect the information com-ing from the number of sources predicted by the model whenestimating the posterior distribution. By doing so, we canhave some insights on the impact of the rate on the analysis.In practice, this can be done by marginalising eq. 11 over N λ using a prior π ( N λ ) ∼ /N λ (Fishbach et al. 2018), whichyields the following expression p ( λ |{ h } k ) ∼ π ( λ ) N obs (cid:89) k =1 (cid:34) (cid:82) L k ( d | θ ) p ( θ | λ ) d θβ ( λ ) (cid:35) , (15)where the integral can be approximated in the same way asin eq. 14 and β ( λ ) is given by eq. 13. Our primary goal is to put constraints on the value of themixing fraction, f ∈ [0 , p ( θ | f, σ Z , σ sp ) = f p ( θ | iso , σ Z , σ sp )+(1 − f ) p ( θ | dyn , σ Z , σ sp )(16)where p ( θ | iso , σ Z , σ sp ) and p ( θ | dyn , σ Z , σ sp ) are the distribu-tions corresponding to the isolated and dynamical models,respectively. A value of f = 0 ( f = 1) indicates that our mixedmodel is composed of BBHs formed only via the dynamical(isolated) channel.We inferred only the distribution of the mixing fractionwhile keeping the other hyper-parameters constant. The anal-ysis was then repeated for all combinations of ( σ Z , σ sp ). Togenerate the posterior distribution for the mixing fraction,we have used a Metropolis-Hastings algorithm with a chainrun for 10 iterations. We discarded the first 10 iterations(burn-in) and estimated the autocorrelation length of ourchains, so that we could trim our original chains to obtainquasi-independent samples of the posterior distribution.From a computational point of view, the form of eq. 16allows us to decompose the number of events N λ , the numberof detections µ ( λ ) and the term in the integral of eq. 11, as twodistinct contributions from the isolated and dynamical models.Our strategy was to compute the values associated with theisolated and dynamical models a priori, for the three terms MNRAS , 000–000 (0000)
Bouffanais et al. p ( f ) p ( f ) Z = 0.2 (no 190814) Z = 0.4 (no 190814) Z = 0.6 (no 190814) Z = 0.2 (with 190814) Z = 0.4 (with 190814) Z = 0.6 (with 190814) isolateddynamical Figure 2.
Posterior distribution of the mixing fraction parameter asinferred from MCMC chains for σ sp = 0 . σ Z = { . , . , . } .The upper panel shows the distributions obtained when marginal-ising over the number of events (eq. 15), while the lower panelincludes the rates (eq. 11). The case in which we include (we do notinclude) GW190814 in the analysis is shown with a solid (dotted)line. Formation channel σ Z N det Isolated 0.2 2Dynamical 0.2 38Isolated 0.4 10Dynamical 0.4 48Isolated 0.6 45Dynamical 0.6 67
Table 1.
Number of detections as a function of formation channeland σ Z for O1, O2 and O3a. For reference, the actual number ofBBH events during O1+O2+O3a is 44. We only consider the eventspresented in Table 1 of Abbott et al. (2020b). mentioned above, so that we were able to easily combine themfor any values of f when running the Monte Carlo MarkovChain (MCMC).Figure 2 shows the mixing fraction posterior distributionin the case σ sp = 0 .
1. The distributions are inferred fromtrimmed MCMCs obtained using the marginalised posteriorfrom eq. 15 (top) and the posterior including the rates ineq. 11 (bottom). The Figure also shows two different analyses,depending if we did or not include the event GW190814 inthe detection set. We will begin by discussing the case wherethe event is not included.First, if we look at the case where rates are not included,we see that the three distributions corresponding to the three σ Z are very similar and Gaussian-like. The medians of thedistributions are equal to 0 .
39, 0 .
36 and 0 .
36 for σ Z = 0 .
2, 0 . .
6, respectively, indicating a slight preference towardsthe dynamical scenario. However, a pure dynamical scenariois outside the 99% credible interval, for which we find lowerbound values equal to 0 .
2, 0 .
18 and 0 .
17 for increasing valuesof σ Z .When taking into account the model rates, we do observea significant change between the distributions, with medianslocated at 0 .
18, 0 .
26 and 0 .
43 for σ Z = 0 . , p ( f ) Z = 0.2 Z = 0.4 Z = 0.6 p ( f ) isolateddynamical Figure 3.
Same as Figure 2, but for σ sp = 0 .
01 (low-spin case).We show only the results we obtain including GW190814 in ouranalysis. the number of expected detections as a function of σ Z . Theisolated model only predicts 2 detections at σ Z = 0 .
2, whichis quite far away from the actual 44 BBH mergers detectedduring O1, O2 and O3a. In comparison, the dynamical modelhas a better prediction of 38 detections, explaining why theposterior distribution of the mixing fraction shifts towardslower values when taking into account rates. In contrast,at σ Z = 0 .
6, the isolated model performs better than thedynamical one, with 38 predicted detections compared to 67.As a result, the posterior distribution is shifted towards highervalues of the mixing fraction.Figure 2 also shows the results obtained when includingthe event GW190814. This GW event is peculiar due tothe low mass of its secondary component, with a medianvalue of m = 2 .
59 M (cid:12) , and is an outlier in the observedmass distribution. Abbott et al. (2020b) show that includingthe event in their analysis has a large impact on the massdistribution and rate inference. In fact, the median value ofthe BBH merger rate density rises from 23 Gpc − yr − whenignoring GW190814, up to 58 Gpc − yr − when the outlier isincluded. In our analysis, the inclusion of the event has onlya little impact on the inference of the mixing fraction. Thishappens because our model distributions are fixed and do notdepend on the set of detected events, unlike in Abbott et al.(2020b), where model distributions do depend on considereddetections. In the analysis presented by Abbott et al. (2020b),an outlier like GW190814 has a large impact on the inferenceof the mass distribution, which then impacts the detectablevolume and the rate inference.We now vary the parameter of the spin magnitude σ sp .Figure 3 shows the mixing fraction for σ sp = 0 .
01. Giventhe very low magnitude of the spin, most of the values for χ eff are very close to 0 regardless of the formation channel,suggesting that the relevant parameters are restricted to themasses’ parameters and redshift. In this case, the posteriordistribution of the mixing fraction obtained marginalising overthe rates shifts towards higher values. The resulting medianvalues are equal to 0 .
49, 0 .
49 and 0 .
57 for σ Z = 0 .
2, 0.4 and0.6. When taking the rates into account, we do observe the
MNRAS000
MNRAS000 , 000–000 (0000) ew insights on BBH formation p ( f ) Z = 0.2 Z = 0.4 Z = 0.6 p ( f ) isolateddynamical Figure 4.
Same as Figure 3, but for σ sp = 0 . same pattern as for σ sp , with a clear differentiation betweenthe three σ Z cases. This is dictated once more by the matchbetween the expected number of detections and the actualnumber of detected events.Finally, Figure 4 shows the results for the high-spin case, inwhich σ sp = 0 .
3. In this case, regardless of whether we includethe rates or not in the analysis, all posterior distributionsgive strong support towards the dynamical formation channel.This springs from the fact that the distribution of χ eff in theisolated case peaks at positive values close to 0 .
3, making itvery difficult for this model to explain events with negativevalues of χ eff .In Figure 5, we repeat the analysis for ( σ sp = 0 . , σ Z = 0 . .
19 for O1, to 0 .
44 for O2 and0 .
16 for O3a. The important shift of the median is under-standable given the relatively low number of events, especiallyfor O1 and O2. For instance, the values of the integral fromeq. 14 that we found for GW150914 is an order of magni-tude larger for the dynamical model compared to the isolatedmodel. As only three GW events were detected during O1,this event alone is driving the posterior distribution towardsa dynamically-dominated mixing model.
We have considered a mixing model between the isolatedformation channel and the dynamical formation channel indense young star clusters. Our results exclude a purely isolatedformation channel for the 44 BBH mergers in GWTC-2, at the90% credible level. According to our fiducial model ( σ sp = 0 . σ sp = 0 . σ sp = 0 .
3) stronglyfavours the dynamical channel with respect to the isolatedone, because most spins are aligned with the orbital angularmomentum in the latter model, while some events in GWTC-2have support for negative values of χ eff . p ( f ) Only O1Only O2Only O3aO1+O2+O3a
Figure 5.
Posterior distribution of the mixing fraction parameterwhen considering each of the observation runs separately. Blue line:O1; orange line: O2; green line: O3a. For reference, the dashedblack line shows the case where all three runs are considered simul-taneously. This result was obtained with σ Z = 0 . σ sp = 0 . One of the key uncertainties of our analysis is the metallicityevolution of the stellar progenitors across cosmic time. Westudy its impact on our results by varying the parameter σ Z ,corresponding to the metallicity spread. Larger values of σ Z strengthen the contribution of the isolated channel (medianvalue of f ≈ . σ Z give more support to the dynamical scenario (median valueof f ≈ . − . m , m , z ), marginalising over the number of events (i.e.computing the posterior distribution as in eq. 15) and adopt-ing a constant value for the spread of metallicity, σ Z = 0 . .
2, but both formation channels are required toproperly explain the observations. Zevin et al. (2020) inves-tigate an ensemble of formation channels that cover isolatedand dynamical formation in globular clusters and nuclearstar clusters. They marginalise over the number of events,use a fixed value for the metallicity spread ( σ Z = 0 .
5) andadopt the same parametrisation as we do: θ = {M c , q, z, χ eff } .When they restrict their analysis to just two channels (i.e.binaries evolved through common envelope and dynamicalbinaries formed in globular clusters), they find that the iso-lated scenario is favoured for most values of their spin-modelparameters, except for the highest spin scenario. Their analy-sis also suggests that a combination of formation channels isnecessary to have the better representation of the observed MNRAS , 000–000 (0000)
Bouffanais et al. events and that the higher spin case favours the dynamicalscenario.To summarize, our main results and these of both Wong et al.(2020) and Zevin et al. (2020) point in the same direction: bothisolated and dynamical scenarios are necessary to properlydescribe current observations from GWTC-2.Wong et al. (2020) and Zevin et al. (2020) focus on oldmassive globular clusters and nuclear star clusters, while wetarget dense young star clusters. This is the first time thatthe young star cluster scenario is compared against the newGWTC-2 data. Young star clusters are generally less massive( ∼ − M (cid:12) ) and shorter lived ( (cid:46) t rlx ∼
10 Myr (cid:18) M SC (cid:12) (cid:19) . (cid:18) r vir (cid:19) / , (17)where r vir is the virial radius. For young star clusters t rlx is afew tens of Myr, while it is several hundreds of Myr (or even afew Gyr) for the typical mass and size of globular clusters andnuclear star clusters. Hence, the stellar progenitors of BHshave enough time to sink to the core of a young star clusterand to dynamically interact with each other, even before theycollapse to BHs (Di Carlo et al. 2019; Kumamoto et al. 2019).In contrast, massive stars die before they sink to the core inboth globular clusters and nuclear star clusters. This explainswhy stellar collisions are more important in young star clustersthan in other star clusters, possibly leading to the formationof BHs with masses >
65 M (cid:12) , such as the ones we consideredhere (Portegies Zwart et al. 2004; Di Carlo et al. 2020a).Secondly, the majority of massive binary stars are hard inyoung star clusters (i.e., they have a binding energy higherthan the average kinetic energy of a star in the cluster, Heggie1975). In contrast, most binary stars are soft in globularclusters and nuclear star clusters. This has a crucial impacton the formation channel of BBHs: most original binary starsbreak in globular/nuclear star clusters because of dynamicalencounters, and most BBHs form from interactions amongthree single BHs (three-body captures, e.g., Morscher et al.2015; Antonini & Rasio 2016). In contrast, most binary starssurvive in young star clusters: they harden by flybys andundergo dynamical exchanges. Hence, most BBHs in youngstar clusters form by dynamical exchanges rather than three-body captures.Thirdly, the escape velocity from young star clusters is ∼ − , significantly smaller than the one of globular ( ∼ − ) and nuclear star clusters ( ∼
100 km s − , Antonini& Rasio 2016). Hence, hierarchical mergers among BBHs areextremely rare in young star clusters, while they are commonin nuclear star clusters. For all of these differences amongyoung, globular and nuclear star clusters, the properties ofBBHs in young star clusters are quite peculiar and deservemore investigation. In a follow-up study, we will compare the populations of BBHs in young star clusters, globular clustersand nuclear star clusters together, in order to remove the biasof having only two formation channels (Zevin et al. 2020).Finally, the mass function of BBHs is certainly one of thekey ingredients of our results. In young star clusters, wehave BBHs with primary masses up to ∼
90 M (cid:12) , while themaximum primary mass in isolated BBH mergers is ∼
40 M (cid:12) .This difference is primarily connected with the collapse ofthe hydrogen envelope of the progenitor star. In tight stellarbinaries, non-conservative mass transfer and common envelopelead to the complete ejection of the hydrogen envelope, beforethe formation of the second BH. Hence, the maximum BHmass in merging isolated binaries is ∼
40 M (cid:12) , correspondingto the maximum helium core mass below the pair instabilitythreshold. In contrast, massive single stars and massive stars inloose binary stars (orbital separation (cid:38) R (cid:12) ) can preservea fraction of their initial hydrogen envelope to the very end,leading to the formation of BHs with masses up to ∼ (cid:12) (e.g., Mapelli et al. 2020; Costa et al. 2021). In isolation,these BHs do not merge, but in star clusters they can pair updynamically and lead to massive BBH mergers. In addition,dynamically triggered collisions between massive stars canlead to even more massive BHs, up to ∼
90 M (cid:12) . While primaryBHs with masses >
60 M (cid:12) represent only ∼ .
5% of all BBHmergers in our simulations (Di Carlo et al. 2020a), they givea crucial contribution to the current analysis.
The LVC has recently published the second GW transientcatalogue (GWTC-2, Abbott et al. 2020a), increasing thenumber of binary compact object mergers from 11 to 50events, most of them BBH mergers. This large number ofevents makes it possible to obtain the first constraints on theformation channels of BBHs.Here, we explore two alternative formation channels: i)isolated binary evolution via stable mass transfer and commonenvelope, ii) dynamical formation in dense young star clusters.Young star clusters are generally less massive and shorterlived than globular clusters, but they are the most commonbirthplace of massive stars (see Portegies Zwart et al. 2010,for a review). Hence, most BHs are likely born in dense youngstar clusters.Comparing our simulations to GWTC-2, we estimate themixing fraction, i.e. the fraction of BBHs formed in youngstar clusters versus isolated binaries. Assuming that the spinmagnitudes follow a Maxwellian distribution, we considerthree different models for the spin, corresponding to a standarddeviation parameter σ sp = 0 . , σ Z =0 . , χ eff in some GWTC-2 events.The metallicity spread σ Z is a key ingredient. Assuming alarge (small) spread tends to favour the isolated (dynamical) MNRAS000
The LVC has recently published the second GW transientcatalogue (GWTC-2, Abbott et al. 2020a), increasing thenumber of binary compact object mergers from 11 to 50events, most of them BBH mergers. This large number ofevents makes it possible to obtain the first constraints on theformation channels of BBHs.Here, we explore two alternative formation channels: i)isolated binary evolution via stable mass transfer and commonenvelope, ii) dynamical formation in dense young star clusters.Young star clusters are generally less massive and shorterlived than globular clusters, but they are the most commonbirthplace of massive stars (see Portegies Zwart et al. 2010,for a review). Hence, most BHs are likely born in dense youngstar clusters.Comparing our simulations to GWTC-2, we estimate themixing fraction, i.e. the fraction of BBHs formed in youngstar clusters versus isolated binaries. Assuming that the spinmagnitudes follow a Maxwellian distribution, we considerthree different models for the spin, corresponding to a standarddeviation parameter σ sp = 0 . , σ Z =0 . , χ eff in some GWTC-2 events.The metallicity spread σ Z is a key ingredient. Assuming alarge (small) spread tends to favour the isolated (dynamical) MNRAS000 , 000–000 (0000) ew insights on BBH formation channel versus the dynamical (isolated). This confirms thatmore observational constraints on the evolution of stellarmetallicities across cosmic time are urgently needed, to narrowdown the uncertainties on BBH merger rates. Despite the largeuncertainties on spin magnitudes and metallicity spread, ourresults point towards an exciting direction: more than oneformation channel is needed to explain the properties of BBHsin the second GW transient catalogue, and the dynamicalpath is essential to account for the largest chirp masses andfor negative values of the effective spin. ACKNOWLEDGEMENTS
MM, YB, FS, UND, NG, SR and GI acknowledge financialsupport from the European Research Council for the ERCConsolidator grant DEMOBLACK, under contract no. 770017.MCA and MM acknowledge financial support from the Aus-trian National Science Foundation through FWF stand-alonegrant P31154-N27. NG acknowledges financial support fromthe Leverhulme Trust Grant No. RPG-2019-350 and RoyalSociety Grant No. RGS-R2-202004.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding authors.
REFERENCES
Abadie J., et al., 2010, Class. Quant. Grav., 27, 173001Abbott B. P., et al., 2016a, Phys. Rev., X6, 041015Abbott B. P., et al., 2016b, Phys. Rev. Lett., 116, 061102Abbott B. P., et al., 2016c, ApJ, 818, L22Abbott B. P., et al., 2016d, ApJ, 833, L1Abbott B. P., et al., 2019a, Physical Review X, 9, 031040Abbott B. P., et al., 2019b, ApJ, 882, L24Abbott R., et al., 2020a, arXiv e-prints, p. arXiv:2010.14527Abbott R., et al., 2020b, arXiv e-prints, p. arXiv:2010.14533Andrews J. J., Cronin J., Kalogera V., Berry C., Zezas A., 2020,arXiv e-prints, p. arXiv:2011.13918Antonini F., Rasio F. A., 2016, ApJ, 831, 187Antonini F., Toonen S., Hamers A. S., 2017, ApJ, 841, 77Arca Sedda M., 2020, The Astrophysical Journal, 891, 47Askar A., Szkudlarek M., Gondek-Rosi´nska D., Giersz M., BulikT., 2017, MNRAS, 464, L36Banerjee S., 2017, MNRAS, 467, 524Banerjee S., 2020, arXiv e-prints, p. arXiv:2004.07382Banerjee S., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 371Bartos I., Kocsis B., Haiman Z., M´arka S., 2017, ApJ, 835, 165Bavera S. S., et al., 2020, arXiv e-prints, p. arXiv:2010.16333Belczynski K., Kalogera V., Bulik T., 2002, ApJ, 572, 407Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A.,Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016,Nature, 534, 512Belczynski K., et al., 2020, A&A, 636, A104Bethe H. A., Brown G. E., 1998, ApJ, 506, 780Bouffanais Y., Mapelli M., Gerosa D., Di Carlo U. N., GiacobboN., Berti E., Baibhav V., 2019, ApJ, 886, 25Bouffanais Y., Mapelli M., Santoliquido F., Giacobbo N., Iorio G.,Costa G., 2020, Constraining accretion efficiency in massivebinary stars with LIGO-Virgo black holes ( arXiv:2010.11220 ) Callister T., Fishbach M., Holz D. E., Farr W. M., 2020, ApJ, 896,L32Choksi N., Gnedin O. Y., Li H., 2018, MNRAS, 480, 2343Choksi N., Volonteri M., Colpi M., Gnedin O. Y., Li H., 2019, ApJ,873, 100Claeys J. S. W., Pols O. R., Izzard R. G., Vink J., Verbunt F. W. M.,2014, A&A, 563, A83Costa G., Bressan A., Mapelli M., Marigo P., Iorio G., Spera M.,2021, MNRAS, 501, 4514Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M.,Wang L., Haardt F., 2019, MNRAS, 487, 2947Di Carlo U. N., Mapelli M., Bouffanais Y., Giacobbo N., Santoliq-uido F., Bressan A., Spera M., Haardt F., 2020a, MNRAS, 497,1043Di Carlo U. N., et al., 2020b, MNRAS, 498, 495Eldridge J. J., Stanway E. R., 2016, MNRAS, 462, 3302Fishbach M., Holz D. E., Farr B., 2017, ApJ, 840, L24Fishbach M., Holz D. E., Farr W. M., 2018, ApJ, 863, L41Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103Fragione G., Kocsis B., 2020, MNRAS, 493, 3920Fragione G., Silk J., 2020, MNRAS, 498, 4591Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., KalogeraV., Holz D. E., 2012, ApJ, 749, 91Fuller J., Ma L., 2019, ApJ, 881, L1Gerosa D., Berti E., 2017, Phys. Rev. D, 95, 124046Gerosa D., Kesden M., Berti E., O’Shaughnessy R., Sperhake U.,2013, Phys. Rev. D, 87, 104028Gerosa D., Berti E., O’Shaughnessy R., Belczynski K., Kesden M.,Wysocki D., Gladysz W., 2018, Phys. Rev. D, 98, 084036Giacobbo N., Mapelli M., 2018, MNRAS, 480, 2011Giacobbo N., Mapelli M., 2019, MNRAS, 482, 2234Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959Giersz M., Leigh N., Hypki A., L¨utzgendorf N., Askar A., 2015,MNRAS, 454, 3150Heggie D. C., 1975, MNRAS, 173, 729Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS,360, 974Hong J., Vesperini E., Askar A., Giersz M., Szkudlarek M., BulikT., 2018, MNRAS, 480, 5645Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897Kalogera V., 2000, ApJ, 541, 319Klencki J., Moe M., Gladysz W., Chruslinska M., Holz D. E.,Belczynski K., 2018, A&A, 619, A77Kroupa P., 2001, MNRAS, 322, 231Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G.,2018, preprint, ( arXiv:1801.05433 )Kumamoto J., Fujii M. S., Tanikawa A., 2019, MNRAS, 486, 3942Kumamoto J., Fujii M. S., Tanikawa A., 2020, arXiv e-prints, p.arXiv:2001.10690K¨upper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011,MNRAS, 417, 2300Lada C. J., Lada E. A., 2003, ARA&A, 41, 57Madau P., Fragos T., 2017, ApJ, 840, 39Mandel I., Farmer A., 2018, arXiv e-prints, p. arXiv:1806.05820Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634Mandel I., Farr W. M., Gair J. R., 2019, MNRAS, 486, 1086Mapelli M., 2016, MNRAS, 459, 3432Mapelli M., 2018, arXiv e-prints, p. arXiv:1809.09130Mapelli M., Giacobbo N., 2018, MNRAS, 479, 4391Mapelli M., Giacobbo N., Ripamonti E., Spera M., 2017, MNRAS,472, 2422Mapelli M., Giacobbo N., Santoliquido F., Artale M. C., 2019,MNRAS, 487, 2Mapelli M., Spera M., Montanari E., Limongi M., Chieffi A., Gia-cobbo N., Bressan A., Bouffanais Y., 2020, ApJ, 888, 76Marchant P., Langer N., Podsiadlowski P., Tauris T. M., MoriyaT. J., 2016, A&A, 588, A50 MNRAS , 000–000 (0000) Bouffanais et al.
McKernan B., et al., 2018, ApJ, 866, 66Miller M. C., Hamilton D. P., 2002, MNRAS, 330, 232Moe M., Di Stefano R., 2017, ApJS, 230, 15Morscher M., Pattabiraman B., Rodriguez C., Rasio F. A., UmbreitS., 2015, ApJ, 800, 9Neijssel C. J., et al., 2019, MNRAS, 490, 3740Nitz A. H., Dent T., Davies G. S., Harry I., 2020, arXiv e-prints, p.arXiv:2004.10015Peters P. C., 1964, Physical Review, 136, 1224Petrovich C., Antonini F., 2017, ApJ, 846, 146Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17Portegies Zwart S. F., Yungelson L. R., 1998, A&A, 332, 173Portegies Zwart S. F., Baumgardt H., Hut P., Makino J., McMillanS. L. W., 2004, Nature, 428, 724Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A,48, 431Rastello S., Mapelli M., Di Carlo U. N., Giacobbo N., San-toliquido F., Spera M., Ballone A., 2020, arXiv e-prints, p.arXiv:2003.02277Rizzuto F. P., et al., 2021, MNRAS, 501, 5257Rodriguez C. L., Loeb A., 2018, ApJ, 866, L5Rodriguez C. L., Chatterjee S., Rasio F. A., 2016a, Phys. Rev. D,93, 084029Rodriguez C. L., Zevin M., Pankow C., Kalogera V., Rasio F. A.,2016b, ApJ, 832, L2Rodriguez C. L., Zevin M., Amaro-Seoane P., Chatterjee S., KremerK., Rasio F. A., Ye C. S., 2019, Phys. Rev. D, 100, 043027Sana H., et al., 2012, Science, 337, 444Santoliquido F., Mapelli M., Bouffanais Y., Giacobbo N., Di CarloU. N., Rastello S., Artale M. C., Ballone A., 2020, ApJ, 898,152Santoliquido F., Mapelli M., Giacobbo N., Bouffanais Y., ArtaleM. C., 2021, MNRAS,Silsbee K., Tremaine S., 2017, ApJ, 836, 39Spera M., Mapelli M., Giacobbo N., Trani A. A., Bressan A., CostaG., 2019, MNRAS, 485, 889Stegmann J., Antonini F., 2020, arXiv e-prints, p. arXiv:2012.06329Stevenson S., Berry C. P. L., Mandel I., 2017a, preprint,( arXiv:1703.06873 )Stevenson S., Berry C. P. L., Mandel I., 2017b, MNRAS, 471, 2801Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946Tagawa H., Haiman Z., Kocsis B., 2019, arXiv e-prints, p.arXiv:1912.08218Tanikawa A., 2013, MNRAS, 435, 1358Tanikawa A., Susa H., Yoshida T., Trani A. A., Kinugawa T., 2020,arXiv e-prints, p. arXiv:2008.01890Tutukov A., Yungelson L., 1973, Nauchnye Informatsii, 27, 70Udall R., Jani K., Lange J., O’Shaughnessy R., Clark J., CadonatiL., Shoemaker D., Holley-Bockelmann K., 2019, arXiv e-prints,p. arXiv:1912.10533Venumadhav T., Zackay B., Roulet J., Dai L., Zaldarriaga M., 2020,Phys. Rev. D, 101, 083030Vigna-G´omez A., Toonen S., Ramirez-Ruiz E., Leigh N. W. C.,Riley J., Haster C.-J., 2021, ApJ, 907, L19Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwen-hoven M. B. N., Naab T., 2015, MNRAS, 450, 4070Wang L., et al., 2016, MNRAS, 458, 1450Wong K. W. K., Breivik K., Kremer K., Callister T., 2020, arXive-prints, p. arXiv:2011.03564Wysocki D., Gerosa D., O’Shaughnessy R., Belczynski K., GladyszW., Berti E., Kesden M., Holz D. E., 2018, Phys. Rev. D, 97,043014Yang Y., Bartos I., Haiman Z., Kocsis B., M´arka Z., Stone N. C.,M´arka S., 2019, ApJ, 876, 122Zackay B., Venumadhav T., Dai L., Roulet J., Zaldarriaga M., 2019,Phys. Rev. D, 100, 023007Zevin M., Pankow C., Rodriguez C. L., Sampson L., Chase E.,Kalogera V., Rasio F. A., 2017, ApJ, 846, 82 Zevin M., et al., 2020, arXiv e-prints, p. arXiv:2011.10057Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS,441, 3703de Mink S. E., Mandel I., 2016, MNRAS, 460, 3545du Buisson L., et al., 2020, arXiv e-prints, p. arXiv:2002.11630MNRAS000