Binary Black Hole Spins: Model Selection with GWTC-3

: The origin of the spins of stellar-mass black holes is still controversial, and angular momentum transport inside massive stars is one of the main sources of uncertainty. Here, we apply hierarchical Bayesian inference to derive constraints on spin models from the 59 most conﬁdent binary black hole merger events in the third gravitational-wave transient catalogue (GWTC-3). We consider up to ﬁve parameters: chirp mass, mass ratio, redshift, effective spin, and precessing spin. For the model selection, we use a set of binary population synthesis simulations spanning drastically different assumptions for black hole spins and natal kicks. In particular, our spin models range from the maximal to minimal efﬁciency of angular momentum transport in stars. We ﬁnd that if we include the precessing spin parameter into our analysis, models predicting only vanishingly small spins are in tension with GWTC-3 data. On the other hand, models in which most spins are vanishingly small but that also include a subpopulation of tidally spun-up black holes are a good match to the data. Our results show that the precessing spin parameter has a crucial impact on model selection.


INTRODUCTION
The third observing run (O3) of the Advanced LIGO (Aasi et al. 2015) and Virgo (Acernese et al. 2015) detectors has brought the number of compact binary merger observations up to 90 events with a probability of astrophysical origin > 0.5 (Abbott et al. 2019(Abbott et al. , 2021d,a,b),a,b).In particular, the 63 confident detections of binary black hole (BBH) mergers (with a false alarm rate FAR< 0.25 yr −1 ) lead to more accurate constraints on the mass and spin distribution of these systems (Abbott et al. 2021c).
Several models have been proposed to infer the spin magnitude of the BH from that of the progenitor star.The main open question concerns the efficiency of angular momentum transport within a star (e.g., Maeder & Meynet 2000;Cantiello et al. 2014;Fuller et al. 2019).If angular momentum is efficiently transferred from the core to the outer layers, mass loss by stellar winds can dissipate most of it, leading to a low-spinning stellar core and then to a low-spinning BH.If instead the core retains most of its initial angular momentum until the final collapse, the BH will be fast spinning.
In the shellular model (Zahn 1992;Ekström et al. 2012;Limongi & Chieffi 2018;Costa et al. 2019), angular momentum is mainly transported by meridional currents and shear instabilities, leading to relatively inefficient spin dissipation.In contrast, according to the Tayler-Spruit dynamo mechanism (Spruit 2002), differential rotation induces the formation of an unstable magnetic field configuration, leading to an efficient transport of angular momentum via magnetic torques.Building upon the Tayler-Spruit mechanism, Fuller & Ma (2019) derived a new model with an even more efficient angular momentum dissipation, predicting that the core of a single massive star might end its life with almost no rotation.
Electromagnetic observations yield controversial results.Asteroseismology favours slowly rotating cores in the late evolutionary stages, but the vast majority of stars with an asteroseismic estimate of the spin are low-mass stars (Mosser et al. 2012;Gehan et al. 2018;Aerts et al. 2019).Continuum-fitting derived spins of BHs in high-mass X-ray binaries are extremely high (e.g., Reynolds 2021;Miller-Jones et al. 2021;Fishbach & Kalogera 2022), but such measurements might be affected by substantial observational biases (e.g., Reynolds 2021).Finally, BH spins inferred from quasi periodic oscillations yield notably smaller values than continuum fitting.For example, the estimate of the dimensionless spin of the BH in GRO J1655-40 is  = 0.7 ± 0.1 and 0.290 ± 0.003 from continuum fitting (Shafee et al. 2006) and quasi-periodic oscillations (Motta et al. 2014), respectively.In a binary system, the evolution of the spin is further affected by tidal forces and accretion, which tend to spin up a massive star, whereas non-conservative mass transfer and common-envelope ejection enhance mass loss, leading to more efficient spin dissipation (Kushnir et al. 2016;Hotokezaka & Piran 2017;Zaldarriaga et al. 2018;Qin et al. 2018).For example, the model by Bavera et al. (2020) shows that the second-born BH can be highly spinning if its progenitor was tidally spin up when it was a Wolf-Rayet star orbiting about the first-born BH.
Furthermore, the orientation of the BH spin with respect to the orbital angular momentum of the binary system encodes information about binary evolution processes.In a tight binary system, tides and mass transfer tend to align the stellar spins with the orbital angular momentum (Gerosa et al. 2018, but see Stegmann & Antonini 2021 for a possible spin flip process induced by mass transfer).If the binary system is in the field, the supernova kick is the main mechanism that can misalign the spin of a compact object with respect to the orbital angular momentum, by tilting the orbital plane (e.g., Kalogera 2000).Finally, the spins of BHs in dynamically formed binary systems are expected to be isotropically distributed, because close encounters in a dense stellar cluster reset any previous signature of alignment (e.g., Rodriguez et al. 2016;Mapelli et al. 2021).
Here, we perform a model-selection hierarchical Bayesian analysis on confident LVK BBHs ( astro > 0.9 and FAR < 0.25 yr −1 ).We consider models of field BBHs for three of the most used angularmomentum transport models: (i) the shellular model as implemented in the Geneva stellar evolution code (Ekström et al. 2012), (ii) the Tayler-Spruit dynamo model as implemented in the code (Cantiello et al. 2014), and (iii) the model by Fuller & Ma (2019).Hereafter, we will refer to these three models simply as GENEVA (G), MESA (M) and FULLER (F) models, following the description in Belczynski et al. (2020).
For each of these models, we consider an additional variation accounting for the Wolf-Rayet (WR) star tidal spin-up mechanism described by Bavera et al. (2020).Also, we account for spin tilts induced by core-collapse supernova explosions.
This paper is organized as follows.Section 2 presents our population-synthesis models.Section 3 describes the hierarchical Bayesian framework we used and discusses the LVK events used in our study.We lay down the results in Section 4, and summarize our conclusions in Section 5.

and natal kicks
We simulated our binary systems with the code (Mapelli et al. 2017;Giacobbo et al. 2018).
We model natal kicks of neutron stars and BHs according to three different models, as shown in Fig. 1: • A unified kick model, in which both neutron stars and BHs receive a kick  kick ∝  ej / rem , where  ej is the mass of the ejecta and  rem the mass of the compact remnant (Giacobbo & Mapelli 2020, hereafter GM20).This model naturally produces low-kicks for electron-capture, stripped and ultra-stripped supernovae (Tauris et al. 2015(Tauris et al. , 2017)).Hereafter, we call this model GM20.
• A model in which compact-object kicks are drawn from a Maxwellian curve with one-dimensional root-mean-square  = 265 km s −1 , consistent with observations of Galactic pulsars (Hobbs et al. 2005).This realistically represents the upper limit for BH natal kicks.Hereafter, we name this model 265.
• A model in which compact-object kicks are drawn from a Maxwellian curve with  = 150 km s −1 .This value of  is more similar to what suggested from indirect measurements of Galactic BH kicks (e.g., Repetto et al. 2017;Atri et al. 2019).Hereafter, we refer to this model as 150.
For more details about , see Giacobbo & Mapelli (2018). is an open-source code and can be downloaded from this link.

Spin magnitude
We have implemented four models for the spin magnitude in , the first three from Belczynski et al. (2020), and the fourth from Bouffanais et al. (2019).Given the large uncertainties on angular momentum transport, we do not claim that these four models are a complete description of the underlying physics: our models must be regarded as toy models, which bracket the current uncertainties on BH spins.

Geneva (G) model
In the Geneva (hereafer, G) model, the dimensionless natal spin magnitude of a BH () can be approximated as: where  = −0.088for all models,  CO is the final carbon-oxygen mass of the progenitor star, while the values of ,  1 ,  2 , and  low depend on metallicity, as indicated in Table 1.This model springs from a fit by Belczynski et al. (2020) to some evolutionary tracks by the Geneva group (Ekström et al. 2012), in which angular momentum transport is relatively inefficient.

MESA (M) model
In the M model, we use the fits done by Belczynski et al. (2020) to a set of stellar tracks run with the code.models the transport of angular momentum according to the Tayler-Spruit magnetic dynamo (Spruit 2002, see also Cantiello et al. 2014).This yields a dimensionless natal BH spin where  1 ,  1 , and  1 are given in Table 2.

Fuller (F) model
Fuller & Ma (2019) predict that angular momentum transport can be even more efficient than the one predicted by the Tayler-Spruit dynamo.Belczynski et al. (2020) summarize the results of the model by Fuller & Ma (2019) simply as  = 0.01 for all single stars and metallicities.

Maxwellian model (Max)
Finally, we also introduce a toy model in which we represent the spin of a BH as a random number drawn from a Maxwellian curve with one-dimensional root-means square   = 0.1 and truncated to  max = 1.0.This model has been first introduced by Bouffanais et al. ( 2019), because it is a good match to the distribution arising from LVK data (e.g., Abbott et al. 2019Abbott et al. , 2021d,c),c).Hereafter, we will indicate this Maxwellian toy model as Max, for brevity.

Tidal spin up
The progenitor star of the second-born BH can be substantially spun-up by tidal interactions.In the scenario explored by Bavera et al. (2020), a common-envelope or an efficient stable mass transfer episode can lead to the formation of a BH-WR binary system, in which the WR star is the result of mass stripping.The orbital period of this BH-WR binary system can be sufficiently short to lead to efficient tidal synchronisation and spin-orbit coupling.The WR star is then efficiently spun-up.If the WR star then collapses to a BH directly, the final spin of the BH will retain the imprint of the final WR spin.
where  WR is the mass of the WR star, while the coefficients  1 ,  2 and  3 have been determined through non-linear least-square minimization and can be found in Bavera et al. (2021).In , we can use these fits for the spin of the second-born BH, while still adopting one of the models presented in the previous subsections (G, M, F, and Max) for the first-born BH.
Table 3. Description of the runs performed for this work. Model for the spin magnitude (Section 2.2).  Correction of the spin magnitude accounting for tidal spin up, as described in B21 (Section 2.3). Model for the natal kick (Section 2.1).

Spin orientation
We assume that natal kicks are the only source of misalignment between the orbital angular momentum vector of the binary system and the direction of BH spins (Rodriguez et al. 2016;Gerosa et al. 2018).Furthermore, we conservatively assume that accretion onto the first-born BH cannot change the direction of its spin (Maccarone et al. 2007).For simplicity, we also neglect the spin-flip process recently described by (Stegmann & Antonini 2021).Under such assumptions, we can derive the angle between the direction of the spins of the two compact objects and that of the orbital angular momentum of the binary system as (Gerosa et al. 2013;Rodriguez et al. 2016) where   is the angle between the new ( ì  new ) and the old ( ì  old ) orbital angular momentum after a supernova ( = 1, 2 corresponding to the first and second supernova), so that cos () = ì  new • ì  old /( new  old ), while  is the phase of the projection of the orbital angular momentum into the orbital plane.

Setup of runs
Hereafter, we consider eight possible models for the spins (see also Table 3): • the first four models (hereafter, G, M, F, and Max) adopt the Geneva, Mesa, Fuller and Maxwellian models for both the first-and second-born BHs, • the other four models (hereafter, G_B21, M_B21, F_B21, and Max_B21) adopt the fits by Bavera et al. (2021) for the second-born BH and the Geneva, Mesa, Fuller and Maxwellian models for the first-born BH.
For each of these eight spin models we consider three different kick models: the GM20, 265, and 150 models discussed in Section 2.1.
As to the main binary evolution parameters, here we use  = 1 for common envelope, while the parameter  depends on the stellar structure as described in Claeys et al. (2014).The other binary evolution parameters are set-up as described in Santoliquido et al. (2021).

Merger rate density
We estimate the evolution of BBH mergers with redshift by using our semi-analytic code C R (Santoliquido et al. 2020(Santoliquido et al. , 2021)).With C R , we convolve our catalogues (Section 2.5) with an observation-based metallicity-dependent star formation rate (SFR) density evolution of the Universe, SFRD(, ), in order to estimate the merger rate density of BBHs as In the above equation,  0 is the Hubble constant, Ω  and Ω Λ are the matter and energy density, respectively.We adopt the values in Aghanim et al. (2020).The term F ( , , ) is given by: where M TOT () is the total simulated initial stellar mass, and dN ( , , )/d () is the rate of BBHs forming from stars with initial metallicity  at redshift  and merging at , extracted from our catalogues.In C R , SFRD(, ) is given by SFRD( , ) = ( ) ( , ), where ( ) is the cosmic SFR density at formation redshift  , and ( , ) is the log-normal distribution of metallicities  at fixed formation redshift  , with average ( ) and spread   .Here, we take both () and () from Madau & Fragos (2017).Finally, we assume a metallicity spread   = 0.3.

Hyper-parametric model description
For each of our models (Table 3), described by their hyper-parameters , we predict the distributions of BBH mergers d d where  are the merger parameters, and   is the total number of mergers predicted by the model.Assuming an instrumental horizon redshift  max = 1.5,   can be calculated as where d c d is the comoving volume and  obs the observation duration.To model the population of merging BBHs, we have chosen five observable parameters  = {M c , , ,  eff ,  p }, where M c = ( 1  2 ) 3/5 /( 1 +  2 ) 1/5 is the chirp mass in the source frame with  1 ( 2 ) the masses of the primary (secondary) BH of the binary, For detectable distribution we mean the distribution of simulated BBHs with sufficiently high signal-to-noise ratio (Section 3).The shaded gray area is the distribution we obtain by stacking the posterior samples of the 59 confident detections from GWTC-3 (Appendix A).  =  2 / 1 .and  is the redshift of the merger.In addition, we used two spin parameters: the effective spin ( eff ) and the precessing spin ( p ).The effective spin  eff is the mass-weighted projection of the two individual BH spins on the binary orbital angular momentum ì where ì  1,2 = ì  1,2 /(  2 1,2 ) is the dimensionless spin parameter of the two BHs.The precessing spin  p is defined as where  1,⊥ ( 2,⊥ ) is the spin component of the primary (secondary) BH perpendicular to the orbital angular momentum vector ì , and  = (4  + 3) /(4 + 3 ).
To compute the distributions (|), we constructed a catalogue of 10 6 sources for all possible combinations of hyper-parameters , using the merger rate density and the metallicity given by C R .From these catalogues we derived continuous estimations of (|) by making use of a Gaussian kernel density estimation assuming a bandwidth of 0.15.

HIERARCHICAL BAYESIAN INFERENCE
Given a set H = {ℎ  }  obs =1 of  obs GW observations, the posterior distribution of a set of hyper-parameters  associated to an astrophysical model can be described as an in-homogeneous Poisson distribution (e.g., Loredo 2004;Mandel et al. 2019;Thrane & Talbot 2019;Bouffanais et al. 2019Bouffanais et al. , 2021a,b),b): where  obs is the number of events observed by the LVK, with an ensemble of parameters ,   is the number of predicted mergers by the model (as calculated in eq.11),   the number of predicted observations given a model and a detector, (,   ) are the prior distributions on  and   , and L  ({ℎ}  |) is the likelihood of the  th observation.
The predicted number of events   can be written in terms of detection efficiency () for a given model: where  det () is the detection probability for a set of parameters .This probability can be inferred by computing the optimal signal to noise ratio (SNR) of the sources and comparing it to a detection threshold.In our case we chose as reference a threshold  thr = 8 in the LIGO Livingston detector, for which we approximated the sensitivity using the measurements for the three runs separately (Abadie et al. 2010;Abbott et al. 2016b;Wysocki et al. 2018).The values for the event's log-likelihood were derived from the posterior and prior samples released by the LVK.Hence, the integral in eq.14 is approximated with a Monte Carlo approach as where    is the  th posterior sample of the  th detection and    is the total number of posterior samples for the  th detection.To compute the prior term in the denominator, we also used Gaussian kernel density estimation.Finally, we can also choose to neglect the information coming from the number of sources predicted by the model when estimating the posterior distribution.By doing so, we can have some insights on the impact of the rate on the analysis.In practice, this can be done by marginalising eq.14 over   using a prior (  ) ∼ 1/  (Fishbach et al. 2018), which yields to the following expression for a model log-likelihood We adopted the formalism described in eqs.14-17 to perform a hierarchical Bayesian inference to compare the astrophysical models presented Sec. 2 with the third gravitational-wave transient catalogue (GWTC-3), the most updated catalogue of gravitational-wave events from the LVK (Abbott et al. 2021b,c).GWTC-3 contains 90 event candidates with probability of astrophysical origin  astro > 0.5.From GWTC-3, we extract 59 confident detections of BBHs with a false alarm rate FAR < 0.25 yr −1 .In this sub-sample, we do not include binary neutron stars and neutron star -BH systems, and we also exclude the other BBH candidates with an higher FAR.Our chosen FAR threshold ensures a sufficiently pure sample for our analysis (Abbott et al. 2021c).A list of the events used in this study is available in Appendix A. For the observable parameters , we use the choice described in Section 2.7, namely  = {M c , , ,  eff ,  p }.

Chirp mass
The chirp mass distribution (Fig. 2) does not depend on the spin model, by construction.Therefore, we only show different natal kicks.Models 150 and 265 show a similar distribution of chirp masses with two peaks of similar importance, one at M c ≈ 8 M and the other (broader) peak at M c ≈ 15 M .In contrast, model GM20 has a much stronger preference for low-mass BHs, with a dominant peak at M c ≈ 8 M .The reason for this difference is that all BHs in tight binary systems receive slow natal kicks in model GM20 (Fig. 1).This happens because stars in tight binary systems lose their envelope during mass transfer episodes; hence, the mass of supernova ejecta ( ej ) is small, triggering low kicks in model GM20.
Figure 2 also compares the detectable distribution of our models with the stacked posterior samples from the confident BBH detections in GWTC-3.This figure highlights two main differences between the population synthesis models and the posterior samples: the peak at M c ≈ 15 M is stronger in the models than it is in the data, while the data present a more significant excess at M c ≈ 25 − 30 M than the models.Finally, the peak at M c ≈ 9 M in the data approximately matches the peak at M c ≈ 8 M in the models.The main features of our population synthesis models (in particular, the peaks at M c ≈ 8 − 10 M and M c ≈ 15 − 20 M ) are also common to other population-synthesis models (e.g., Belczynski et al. 2020;van Son et al. 2022) and mostly spring from the core-collapse SN prescriptions by Fryer et al. (2012).Alternative core-collapse SN models (e.g., Mapelli et al. 2020;Mandel et al. 2021;Patton et al. 2022;Olejak et al. 2022) produce different features and deserve further investigation (Iorio et al., in prep.).

Spin parameters
Figure 3 shows the detectable distribution of spin parameters  p and  eff for all of our models.By construction, large spins are much more common in models G and G_B21, while models F and F_B21 have a strong predominance of vanishingly small spins.Models M, M_B21, Max and Max_B21 are intermediate between the other two extreme models.
Including or not the correction by B21 has negligible impact on the distribution of  p and  eff for models G, because of the predominance of large spin magnitudes.In contrast, introducing the spin-up correction by B21 has a key impact on models F, because it is the only way to account for mild to large spins in these models.The correction by B21 is important also for models M and Max, being responsible for the large-spin wings.
Finally, our model with slow kicks (GM20) results in a distribution of  p that is more peaked at zero (for models G, M and Max) with respect to the other two kick models (150 and 265).In fact, the supernova kicks in model GM20 are not large enough to appreciably misalign BH spins (see Fig. 1).
A similar effect is visible in the distribution of  eff : model 265 produces a distribution of  eff that is less asymmetric about the zero with respect to models 150 and especially GM20.

Model Selection
Figure 4 and Table 4 report the values of the log-likelihood log L defined in Eq. 17.We can quantify the difference between two models A and B by computing the average absolute difference in percentage on the non-A,B variation  ( would be kick(spin) if A and B are spin(kick) models).For example to compare the two models G and G_B21, A and B become G_B21 and G and  = {GM20, 150, 265}.
The tidal spin-up mechanism (B21) affects the spin of a small part of the population of each model (Fig. 3).However, it improves the likelihood of the F and M models significantly (e.g., ΔlogL (M_B21, M) = 89%, Table 4).This improvement of the loglikelihood can be explained by the presence of higher values of  p and  eff in the distribution of populations M_B21 and F_B21 compared to M and F (Fig. 3).
The F model yields L (F) = −∞ if we do not include the tidal spin-up correction, regardless of the kick model.This indicates that the LVK data do not support vanishingly small BH spins for the entire BBH population.However, it is sufficient to inject a tiny subpopulation of spinning BHs, by switching on the B21 correction, and the F model becomes one of the best considered models.In fact, the F_B21 models only includes 0.4% of BHs with  > 0.01 and achieves log L > 200 (for spin models 150 and 265).
The G and G_B21 spin models exhibit lower log-likelihood values than the others for all kicks models: logL 150 for 150 and 265, and logL < 0 for GM20.This happens because the distribution of  eff has non-negligible support for extreme values  eff < −0.5 and  eff > 0.5 (Fig. 3).
The kick models 150 and 265 show similar results (ΔlogL (150, 265) < 3%) for every spin assumptions.Also, for all spin assumptions, the GM20 kick model scores a significantly lower likelihood than the other models 150 and 265 with ΔlogL (150, GM20) ∼ ΔlogL (265, GM20) ∼ 150%.This result can be explained by the high peak of model GM20 at low chirp masses (M c ∼ 8M , see Sec.4.1 and Fig. 2) and by the low value of  p compared to the other kick models (Fig. 3).Models Max and Max_B21 are possibly the best match to the data, but this is not surprising, because they were built as a toy model to visually match the data.Among the astrophysically-motivated models (i.e., after excluding the Max model), M, M_B21 and F_B21 (with kick models 150 and 265) are the most favoured by the data.This might be interpreted as a support for the Tayler-Spruit instability mechanism (adopted in models M) and for the tidal spin-up model by B21.

Importance of 𝜒 p
The  p parameter encodes information on the spin component in the orbital plane.Its impact on gravitational-wave signals is much lower than that of  eff , and therefore its measurement is less precise.To understand the impact of  p on our results, we re-ran the analysis without this parameter.The results are shown in Table 5 and in Fig. 4 with empty markers.Fig. 4 shows that, if we do not include  p , the models M and M_B21 have almost the same log-likelihood, and even the F model yields a positive log-likelihood.Furthermore, the analysis without  p results in significantly larger values of L for the kick model GM20.Our results demonstrate that the measured  p of GWTC-3 BBHs carries substantial information, despite the large uncertainties.

DISCUSSION
The spin magnitude of BHs is largely uncertain, mostly because we do not fully understand angular momentum transport in massive stars.Here, we have taken a number of spin models bracketing the main uncertainties, we have implemented them into our populationsynthesis code , and compared them against GWTC-3 data within a hierarchical Bayesian framework.
The data do not support models in which the entire BH population has vanishingly small spins (model F).This result is mainly driven by the  p parameter.This is in agreement with, e.g., the complementary analysis presented in Callister et al. (2022).They employed a variety of complementary methods to measure the distribution of spin magnitudes and orientations of BBH mergers, and concluded that the existence of a sub-population of BHs with vanishing spins is not required by current data.Callister et al. (2022) 2022) claimed the existence of a subpopulation of zero-spin BHs.From our analysis, we cannot exclude the existence of such sub-population, as the F model with B21 correction (F_B21) still represents a good match of the data.Similarly to Belczynski et al. (2020) and Gerosa et al. (2018), we find that models with large spins (G, G_B21) are less favoured by the data, but they are still acceptable if we allow for large kicks.
Overall, we find a preference for large natal kicks.This result goes into the same direction as the work by Callister et al. (2021).Actually, this preference for large natal kicks is degenerate with the adopted formation channel.Had we included the dynamical formation channel in dense star clusters, we would have added a sub-population of isotropically oriented spins (see, e.g., Figure 8 of Mapelli et al. 2022).In a forthcoming study, we will extend our analysis to a multi-channel analysis.While it is unlikely that BBH mergers only originate from one single channel, adding more formation channels to a hierarchical Bayesian analysis dramatically increases the number of parameters, making it more difficult to reject some portions of the parameter space.

SUMMARY
The origin of BH spins is still controversial, and angular momentum transport inside massive stars is one of the main sources of uncertainty.Here, we apply hierarchical Bayesian inference to derive constraints on spin models from the 59 most confident BBH merger events in GWTC-3.We consider five parameters: chirp mass, mass ratio, redshift, effective spin, and precessing spin.
For model selection, we use a set of binary population synthesis simulations spanning different assumptions for black hole spins and natal kicks.In particular, our spin models account for relatively inefficient (G), efficient (Max and M), and very efficient angularmomentum transport (F).A higher efficiency of angular momentum transport is associated with lower BH spins.In particular, model F predicts vanishingly small spins for the entire BH population.For each of our models, we also include the possibility that some BHs are tidally spun-up (B21).We considered three different natal kick models: according to models 265 and 150, we randomly draw the kicks from a Maxwellian curve with  = 265 and 150 km s −1 , respectively; in the third model (G20), we also derive the kicks from a Maxwellian curve with  = 265 km s −1 , but the kick magnitude is then modulated by the ratio between the mass of the ejecta and the mass of the BH.
We summarize our main results as follows.
• The data from GWTC-3 do not support models in which the entire BH population has vanishingly small spins (model F).
• In contrast, models in which most spins are vanishingly small, but that also include a sub-population of tidally spun-up BHs (model F_B21) are a good match to the data.
• The models in which angular momentum transport is relatively inefficient (G and G_21) yield log-likelihood values that are much lower than models with efficient angular momentum transport (M, M_B21, Max, and Max_B21).
• Models with large BH kicks (150 and 265 ) are favoured by our analysis with respect to low-kick models (G20).
• Our results show that the precessing spin parameter  p plays a crucial impact to constrain the spin distribution of BBH mergers.

Figure 1 .
Figure1.Probability distribution function (PDF) of the binary kick velocities in the centre of mass ( CM ), for our sample of simulated BBH mergers.The centre-of-mass kick velocity takes into account both the first and the second supernova event in each binary system(Perna et al. 2022).Dashed dark-cyan line: model GM20; solid black line: 150; dotted red line: 265.This figure only shows the kick velocity of the stellar progenitors of BBHs that merge within the lifetime of the Universe.
Based on the simulations byBavera et al. (2020),Bavera et al. (2021) derive a fitting formula to describe the spin-up of the WR star and the final spin of the second-born BH:  =  WR log 2 10 (/[day]) +  WR log 10 (/day) if ≤ 1 d 0 otherwise, (3) where  is the orbital period of the BH-WR sytem,  WR =   WR ,   1 ,   2 ,   3 and  WR =     ,

Figure 2 .
Figure2.Predicted detectable distribution of chirp mass, for each kick model: GM20 (solid dark-cyan line), 150 (dotted black line) and 265 (dashed red line).For detectable distribution we mean the distribution of simulated BBHs with sufficiently high signal-to-noise ratio (Section 3).The shaded gray area is the distribution we obtain by stacking the posterior samples of the 59 confident detections from GWTC-3 (Appendix A).

Figure 3 .
Figure 3. Predicted detectable distribution of  p (left) and  eff (right) for all of our models.Different colours refer to the spin model: G, M, F and Max.Solid (dashed) lines include (do not include) the tidal spin-up model by B21.From top to bottom: GM20, 150, and 265.The shaded gray areas are the distributions we obtain by stacking the posterior samples of the 59 confident detections from GWTC-3 (Appendix A).

Table 1 .
Parameters adopted in model G. See Eq. 1 for details.

Table 2 .
Parameters adopted in model M. See Eq. 2 for details.