Inferring black hole formation channels in GWTC-4.0 via parametric mass-spin correlations derived from first principles

CERN-TH-2025-240 Inferring black hole formation channels in GWTC-4.0 via parametric mass-spin correlations derived from first principles Emanuele Berti , 1, ∗ Francesco Crescimbeni , 2, 3, † Gabriele Franciolini , 4, 5, 6, ‡ Simone Mastrogiovanni , 3, § Paolo Pani , 2, 3, ¶ and Grégoire Pierra 3, ∗∗ 1 Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA 2 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy 3 INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185, Roma, Italy 4 Dipartimento di Fisica e Astronomia “G. Galilei”, Università degli Studi di Padova, via Marzolo 8, I-35131 Padova, Italy 5 INFN, Sezione di Padova, via Marzolo 8, I-35131 Padova, Italy 6 Department of Theoretical Physics, CERN, Esplanade des Particules 1, P.O. Box 1211, Geneva 23, Switzerland (Dated: December 4, 2025) We investigate the differences between several proposed formation scenarios for binary black holes (BBHs), including isolated stellar evolution, dynamical assembly in dense clusters and AGN disks, and primordial BHs. Our approach exploits the predicted spin features of each formation channel, and adopts parameterized models of the predicted correlations between the spin magnitudes (and orientations) and mass, inspired by first principles. Using hierarchical Bayesian inference on the recent GWTC-4.0 dataset, we compare these features across all models and assess how well each scenario explains the data. We find that the data strongly favor the presence of a positive correlation between mass and spin magnitude, in agreement with previous studies. Furthermore, the hierarchical scenario provides a better fit to the observations, due to the inclusion of second-generation mergers leading to higher spins at larger masses. The current dataset is not informative enough about spin orientation: the cluster (random orientations) and AGN (aligned orientations) scenarios have comparable Bayesian evidence. Finally, the mass-spin correlation predicted by the primordial scenario gives a poor fit to the data, and this scenario can only account for a subset of the observed events. CONTENTS I. Introduction 1 A. Executive summary 2 II. Modeling different BBH populations 2 A. Spin models 4 1. Isolated channel 4 2. Hierarchical channel 4 3. Active galactic nuclei channel 6 4. Primordial channel 6 B. Redshift model 6 III. Hierarchical Bayesian inference set-up 7 IV. Results 8 A. Model comparison 9 B. Posterior predictive distributions 10 1. Spin parameters 12 V. Conclusions 14 Acknowledgments 15 A. The LVK Gaussian component spin model 16 ∗[email protected] †francesco.crescimb[email protected] ‡[email protected] §simone.mastrogiov[email protected] ¶[email protected] ∗∗ [email protected] B. Impact of effective likelihood variance cuts 16 C. Mass distribution parameters 17 References 18 I. INTRODUCTION Since the first gravitational-wave (GW) detection by the LIGO-Virgo Collaboration [ 1 ], the growing catalog of GW candidates has provided increasingly precise constraints on the astrophysical population of merging binary black holes (BBHs) and binary neutron stars [ 2 ]. With the expansion of the search volume to cosmological distances, it is now possible to investigate the statistical properties of these sources and compare them against different models of stellar and binary evolution, as well as more exotic scenarios. A central goal of population studies is to determine how compact binaries form and evolve before merger [3– 5 ]. Proposed astrophysical formation channels include for example isolated binary evolution in galactic fields (see e.g. [ 6 – 15 ]), where processes such as mass transfer and common-envelope evolution [ 16 ] shape the final system, and dynamical assembly in dense stellar environments, like globular clusters or galactic nuclei [ 17 – 27 ] (see [ 28 ] for a review). These different pathways are expected to leave characteristic imprints on the observed mass, spin, eccentricity and redshift distributions of the merging binaries, which can be probed through hierarchical Bayesian analysis [ 29 – 31 ]. The growing number of detections enables arXiv:2512.03152v1 [gr-qc] 2 Dec 2025 2 more detailed tests of stellar evolution predictions, such as the existence of features in the BBH mass distribution associated with pair-instability supernovae [ 32 – 35 ], and the effects of tidal interactions and of the binary’s accretion history on the spin distributions [ 36 ]. Measurements of the merger rate evolution with redshift further connect compact binaries to the cosmic star-formation history. In addition, GW observations offer a unique window onto more exotic possibilities, such as primordial BHs (PBHs) formed in the early Universe [ 37 – 40 ], which may contribute to the observed merger rate and provide clues about physics beyond the standard model of cosmology, including their potential role as dark matter candidates or tracers of high-redshift phenomena. One interesting possibility to study and characterize different formation channels is through the spin distribution of merging BBHs [ 41 – 44 ]. The third observing run of the LIGO-Virgo-KAGRA (LVK) Collaboration [ 45 ] found that most BBH systems are produced with spin magnitudes that have strong support for χ≲ 0 . 4, while the distribution of tilt angles θ seems to prefer systems with spins above the orbital plane. The GW candidates in the GWTC-4.0 catalog [ 2 ] support these conclusions, and also hint at a more detailed structure in the distribution of the effective spin, χeff . In particular, the spin magnitude distribution is concentrated at χ≲ 0 . 4, the spin-tilt distribution may peak away from perfect alignment with the orbital angular momentum, and the χeff distribution is asymmetric around its peak. There have been several attempts to constrain the population distributions of GW events by leveraging spin measurements to infer the underlying BH merger channels, both with parametric models [ 46 – 60 ] and non-parametric ones [ 61 – 69 ] (see e.g. [ 70 ] for an in-depth review of the topic). Having an accurate inference of BH spins, based on a physically motivated parametric models, can potentially result into more robust constraints on the axion parameter space through the phenomenon of BH superradiance (see e.g. [71–75]). In this work, we build upon this approach by modeling the spin distribution and its correlation with mass with physically motivated, simplified models derived from first principles for four of the main formation scenarios: BHs formed in isolation (IBHs), BHs formed hierarchically in clusters (HBHs), BHs formed hierarchically in the disks of active galactic nuclei (AGNs), and primordial BHs (PBHs). These simplified yet insightful models are designed to capture the key features of distinct BBH formation channels, and to provide a more direct physical connection between spin measurements and their astrophysical origins. A. Executive summary Our main findings are based on the modeling of physical correlations between source masses, spin magnitudes, and tilt angles inspired by the four different formation channels listed before, and detailed in Sec. II A below. Here, for the reader’s convenience, we summarize the main results of the paper: • We find strong support for a spin magnitude distribution which broadens at high masses. Within our models, this is closer to hierarchical scenarios (namely, HBHs and AGNs), which include both first and second-generation mergers. However, there is no support either in favor of or against a flat spin direction distribution, compared to spin-angular momentum alignment at small masses. Therefore, we cannot distinguish, among the environments we have considered (clusters or AGN disks), the one in which hierarchical mergers are most likely to occur. • Our analysis shows a weak preference for multiple spin populations, although we observe that even a single hierarchical scenario—either HBHs or AGNs— could be able to fully capture the spin distribution and its correlation with the mass observed in the GWTC-4.0 data. • The sharp mass-spin correlation predicted by the primordial BH scenario, with efficient cosmological mass-spin evolution, is strongly disfavored as the sole explanation of the GWTC-4.0 dataset. • The inferred mass-distribution parameters do not change significantly (i.e., beyond the O (1) σ level) when performing the inference across the different models, or combinations of models, considered in this work, despite their different predictions for the spin distribution. • The information on the merger-rate redshift evolution remains subdominant compared to spin information in the current catalog, and/or the data do not require each subpopulation to follow radically different redshift distributions. • For scenarios that provide a poor fit to the data, such as the PBH-only case, the hierarchical likelihood is evaluated in regions of parameter space where the Monte Carlo integrals used for the event posteriors or for the selection effects may be insufficiently sampled. Accurate stability estimators should be used to ensure proper convergence. II. MODELING DIFFERENT BBH POPULATIONS In this work, we explore different formation scenarios for BBHs that could explain the physical properties of the population of observed GW sources. We focus on four distinct classes of models: • IBHs: binaries formed through isolated stellar evolution in galactic fields. 3 101102 m1[M] 0.0 0.2 0.4 0.6 0.8 1.0 χ 101102 m2[M] 101102 m1[M] −1.00 −0.75 −0.50 −0.25 0.00 0.25 0.50 0.75 1.00 cos θ 101102 m2[M] FIG. 1: GWTC-4.0 catalog: Mass-spin scatter plot of the 153 GW candidates used in this analysis, chosen to have at least IFAR=1yr −1 . The stars represent the median values of the BBH event parameters from the GWTC-4.0 catalog. The x-axis indicates the source-frame masses of the primary (left columns) and secondary (right columns), while the y-axis shows the dimensionless spin magnitudes χ (first row) or the cosine of the polar angle cos θ (second row). Error bars correspond to the 1 σ uncertainties from the official LVK parameter estimation samples for each event. We do not show the non-trivial correlation between parameters in the posterior for simplicity. mi indicates source frame mass, as in the rest of the text. • HBHs: binaries assembled dynamically in dense stellar environments, such as globular clusters, where previous merger remnants can participate in subsequent (hierarchical) mergers. • AGNs: binaries assembled dynamically in AGN disks. This channel is structurally similar to the HBH one, with the only exception that BHs form in a disk environment, which imprints a preferred aligned-spin direction. • PBHs: binaries formed in the early Universe from the collapse of primordial density fluctuations, independent of stellar processes. Each channel is characterized by its own mass, spin, and redshift distributions, which can be derived from first principles—that is, obtained directly from the physical properties predicted by the channel itself, rather than relying on phenomenological or data-driven assumptions. While population studies based on distinct mass and redshift distributions in each formation channel remain the primary approach for identifying the nature of the merger populations [ 76 – 80 ], as these properties are more tightly constrained by the GW data, the spin distribution can offer valuable complementary information. Indeed, the most salient spin properties can be modeled agnostically and the prediction of every channel may be more robust, despite the details of the formation mechanism not being fully understood. For example, isolated binary formation generically leads to preferentially aligned spins [ 36 , 81 – 86 ], whereas hierarchical mergers predict a subpopulation of spinning BBHs clustering around χ∼ 0 . 7, which may dominate at large masses [ 20 , 21 , 31 , 87 ]. Similar arguments also apply to the primordial scenario. For PBHs, the mass distribution, and even its overall range, is essentially unknown [ 88 ] and strongly model dependent, whereas the main feature of their spin, namely that PBHs are formed with nearly zero 4 spin in the standard formation scenario [ 89 – 91 ] and can acquire spin only through mass-dependent accretion [ 92 – 94 ], is robust and generic. In our analysis, our goal is to disentangle the different possible formation channels using only information from the mass-spin correlations, specific to each scenario. In Fig. 1we show the spin magnitudes and orientations of the observed (i.e., after imposing the selection effect) primary and secondary source-frame masses of all the GWTC-4.0 BBH candidates with IFAR=1 yr −1 , which selects NBBH =153 events. Some visible trends can be identified from the figure, such as a positive correlation between mass and spin for the primary BH, a sparsely populated region at negative orientations for low masses, and the outstanding event GW231123 [ 95 ], which is currently the most massive and most rapidly spinning system detected to this day (but see Ref. [ 96 ] for a discussion, as the origin of this event is still debated [97–106]). In the following sections, we introduce the phenomenological parametric models for the spin and redshift distributions used in this work. We stress that we adopt an agnostic, flexible parametric model for the merger rate, without imposing astrophysical priors. Leveraging the existing uncertainties on the astrophysical and primordial rates, we assume a logarithmic flat prior on the rate of each population. As for the mass spectrum, we assume a single common distribution, i.e., we do not model the mass distribution of each channel separately. Although this is admittedly a strong assumption, it remains compatible with current population-synthesis predictions, which show that isolated, multigenerational, and/or primordial BBHs—either individually or in combination—can reproduce the observed phenomenological mass function. Also, this assumption is analogous to some LVK studies [ 2 ], where multi-population analyses are performed with a single mass distribution. The overall mass distribution we adopt, following the LVK analysis, is discussed in Appendix C. In some cases we will also allow for different merger-rate evolutions for each subpopulation, and find that redshift information remains largely subdominant. A. Spin models In this subsection, we describe in detail the physicallyinformed models adopted for the spin distributions as functions of the masses. For each population model, we provide an illustrative example of the mass-conditional spin and tilt-angle distributions in Fig. 2. We note that the spin distributions we employ, depending on the specific formation scenario, involve mass-dependent parameters, as detailed below. This introduces characteristic massspin correlations in the BBH population. 1. Isolated channel In isolated binaries, spin orientations retain memory of the stellar progenitors’ evolutionary history [ 36 , 81 – 86 ]. While stellar spins are often assumed to be initially aligned with the orbital angular momentum, misalignment can be introduced by supernova kicks imparted at BH formation, which tilt the orbital plane [ 36 , 83 , 107 – 109 ]. Tidal interactions during the binary’s evolution can partially counteract this effect by realigning spins [ 82 , 85 , 110 ]. We will remain agnostic about the detailed modeling of these processes and capture the main properties of this scenario as follows. We assume that BH spins are preferentially aligned with the orbital angular momentum. The tilt angles θi (for i = 1 , 2) are modeled such that cos θi follows a Gaussian distribution centered at unity, representing a preference for alignment, with a standard deviation parameter δ controlling the degree of misalignment. The spin magnitudes χi are drawn independently from a Gaussian distribution N[0,1] ( χi| 0 , χmax )centered at zero, indicating a preference for small spin values (see e.g. [ 31 ]), and truncated to the physical range [0 , 1]. The resulting distributions are then given by p(cos θi)=N[−1,1](cos θi|1, δ),(1) p(χi)=N[0,1](χi|0, χmax).(2) The correlation between source mass and spins is encoded in χmax and δ . For simplicity, we model the δ ( m )and χmax ( m )dependence as linear functions within the relevant range (an approximation that should be sufficient at the present level of measurement precision), and take (e.g. [111]) δ=δIBH 0+˙ δIBH m 30M⊙[0.1,2] ,(3) χmax =χIBH 0+ ˙χIBH m 30M⊙[0.1,1] ,(4) where m is the individual source-frame BH mass. We constrain the parameters to be within the physical range δ∈ [0 . 1 , 2] and χIBH max ∈ [0 . 1 , 1], respectively. Hence, in the IBH scenario, the spin model hyperparameters are ΛIBH ={δIBH 0,˙ δIBH, χIBH 0,˙χIBH}.(5) 2. Hierarchical channel For the dynamical (hierarchical) formation channel, we consider a population which can include contributions from hierarchical mergers. In this model, the BBH population is composed of three distinct sub-populations [ 20 , 31 ]: • 1g + 1g (first-generation mergers): both BHs originate from stellar collapse, and later they form a binary by gravitational dynamics (i.e., capture and multi-body interactions). 5 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 4 5 6 -1.0 -0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 FIG. 2: Examples of probability distributions of χ (first row) and cos θ (second row), for the different formation channels. The distributions span from 10 M⊙ (blue) to 100 M⊙ (red). The model parameters of each scenario are fixed to the maximum likelihood (ML) point obtained when fitting the GWTC-4.0 catalog with our models. • 2g + 1g or 1g + 2g (mixed-generation mergers): one component is a remnant from a previous merger. • 2g + 2g (second-generation mergers): both components are merger remnants. We neglect N > 2generation mergers, as their current detection rate is expected to be subdominant with respect to the former generations. Indeed, since 2g BHs inevitably have large spins, they receive larger merger recoils, and only clusters with very high escape velocity can successfully retain a meaningful fraction of 3g or higher-generation mergers (see e.g. [112–115]). For all the spin orientations in this scenario, we assume that the angles θi are isotropic, so cos θi is drawn uniformly in the physically allowed range [−1,1], and p(cos θi) = 1 2,cos θi∈[−1,1].(6) For first-generation binary components, we assume the spin magnitudes χi to be drawn from a Gaussian distribution centered at χ = 0 with a variable width χmax . Therefore, similarly to the IBH channel, we assume p1g(χ)=N[0,1](χi|0, χmax).(7) For binary components of second generation, we assume a spin magnitude drawn from a truncated Gaussian centered at ¯χf (which is ¯χf≈ 0 . 69 for spinless binaries [ 116 , 117 ]) p2g(χ)=N[0,1](χ|¯χf, σf),(8) where σf is a hyperparameter describing the width of the distribution of the remnant spin (which is a relatively narrow distribution around the spinless results, see e.g. Fig. 8 of Ref. [ 118 ]), and the Gaussian is truncated to lie within the physical range [0,1]. To account for all the sub-populations with relative mixing fractions, we model our total spin distribution as p(χ|m) = X x,y=1,2 πx(m1)πy(m2)pxg(χ1)pyg(χ2),(9) where m = ( m1, m2 )denotes the component masses, χ = ( χ1, χ2 )the spin magnitudes, and  θ = ( θ1, θ2 )the tilt angles. We further define πx(m)as πx(m) = (f(m)if x= 1 1−f(m)if x= 2 ,(10) with f(m) = f1g,high +(f1g,low −f1g,high) 1+e m−mt δmt ,(11) where f ( m )is a sigmoid mixture fraction that switches from different fractions of 1g/2g BHs as a function of mass. As 2g BHs should be typically more massive than 1g ones, we expect the fraction f ( m )to transition between 1 and 0. In the above equation, f1g,high is the fraction of 1g BHs at high masses, f1g,low is the fraction of 1g BHs at low masses, mt is a transition mass, and δmt is a transition window. Note that Eq. (10) allows for the possibility that the more massive black hole is 1g and the secondary 2g. 6 This model of the spin distribution of 1g+2g BHs can reproduce the χeff distribution found in populationsynthesis results for GC and NSC environments, as defined in Ref. [ 76 ] based on [ 119 , 120 ], simply by adjusting χmax and the relative fraction as a function of the primary mass. This was shown explicitly in Ref. [ 46 ] (see their Fig. 1). Therefore, this model provides a well-motivated framework to capture the spin properties, and their correlation with mass, in the dynamical scenario. Having fixed {¯χHBH f = 0 . 69 , fHBH 1g,low = 1 , fHBH 1g,high = 0 } to their values as motivated above, the population hyperparameters in the HBH model are ΛHBH ={χHBH max , σHBH f, mHBH t, δmHBH t}.(12) 3. Active galactic nuclei channel We also include a simple model for the AGN channel (see e.g. [ 19 , 121 ]). For the spin magnitude, this channel is assumed to contain both 1g and 2g mergers, with a spin distribution matching the form adopted for the HBH model in Sec. II A 2. However, for this model, the presence of a disk in the AGN environment generates a preferential direction for both BH pair dynamics and gas accretion [ 122 , 123 ], thus leading to binaries that have nearly aligned spins [ 115 , 121 ]. Therefore, we model the tilt distributions as: p(cos θi)=N[−1,1](cos θi|1, δAGN).(13) We allow for the tilt angles to be correlated with masses, therefore: δAGN =δAGN 0+˙ δAGN m 30M⊙[0.1,2] .(14) Following a similar treatment for the spin magnitude as described above, the AGN model is described by the following model hyperparameters: ΛAGN ={δAGN 0,˙ δAGN, χAGN max , σAGN f, mAGN t, δmAGN t}. (15) Accretion effects in the AGN channel are not explicitly included in this model, and are left for future work. 4. Primordial channel In the PBH scenario, binaries are supposed to form in the early Universe with negligible initial spins, due to the quasi-spherical nature of the collapse of large curvature perturbations during radiation domination [ 89 – 91 ]. However, PBHs may acquire spin through gas accretion before re-ionization [ 92 , 93 ], with the efficiency of accretion depending on the PBH mass. Accretion is negligible for light PBHs below O (10) M⊙ , while it can induce significant spins for more massive PBHs (see the upper right panel of Fig. 2), thus introducing a characteristic correlation between spin and mass. The location of this transition depends on the accretion efficiency, and it is encoded in the model hyperparameter zcut−off [ 93 ]. In the absence of subsolar mass merger detections [ 124 – 128 ] and access to high redshfit events with z≳O (30) [ 129 – 135 ], the mass-spin correlation induced by accretion remains the only predictive testable imprint of the primordial scenario [136]. The spin directions are assumed to be isotropically distributed, as expected for independent PBHs with random relative orientations [92]: p(cos θi) = 1 2,cos θi∈[−1,1].(16) The spin magnitudes χi are modeled via an analytical fit describing accretion-driven spin growth as a function of the primary mass m1 , mass ratio q , and a cutoff redshift zcut-off that parametrizes the end of efficient accretion. A detailed description of these functions can be found in Ref. [ 136 ]. Given the uncertainties in the accretion efficiencies, as well as expected scattering in the environmental properties around PBH binaries, we allow both spins to be distributed as a Gaussian centered around the value χ = χPBH ( m|zcut-off )predicted in the model, so that p(χ|m)=N[0,1](χ|χPBH(m|zcut−off ), σχ).(17) In this case, the hyperparameters of the PBH model are ΛPBH ={zcut-off, σχ},(18) as the full spin distribution is determined by these parameters via the analytical fit described above. We summarize the spin model parameters in Table I. B. Redshift model For all channels, we model the merger rate density evolution with redshift through a smooth Madau-Dickinson broken power-law function, ψ(z|Λz) = "1 + 1 (1+zp)γ+k#(1+z)γ 1 + h(1+z) (1+zp)iγ+k,(19) governed by three parameters Λz≡ {γ, k, zp}.(20) This shape is motivated by the star formation rate evolution [ 137 ], although it remains sufficiently flexible, so that it can fit the low-redshift merger rate at z≲O (3)—of relevance for current LVK sensitivity—for all the formation channels considered here. Interestingly, the functional form in Eq. (19) can reproduce specific population-synthesis predictions for the scenarios described above. In particular, as shown in 7 TABLE I: Priors on the BBH spin and redshift-evolution model hyperparameters. U indicates uniform distribution, while LU indicates log-uniform distribution, within the indicated range. Spin models Parameter Prior IBH δIBH 0U[0.1,1] ˙ δIBH U[−1,1] χIBH 0U[0.05,1] ˙χIBH U[−1,1] HBH χHBH max U[0.1,1] σHBH fU[0.1,0.5] mHBH tU[10,100] M⊙ δmHBH tU[1,100] M⊙ AGN δAGN 0U[0.1,1] ˙ δAGN U[−1,1] χAGN max U[0.1,1] σAGN fU[0.1,0.5] mAGN tU[10,100] M⊙ δmAGN tU[1,100] M⊙ PBH zcut−off U[10,30] σχU[0.05,0.2] Redshift evolution Parameter Prior range IBH, HBH, AGN γU[0,10] kU[0,10] zpU[0,10] PBH γPBH 1.17 kPBH −γPBH All models Rc 0LU[10−4,80] Ref. [ 138 ], the isolated channel predicts a merger rate with reference parameters [139] ΛIBH z={γ= 2.57 , k = 3.26 , zp= 2.36}.(21) For the case of HBHs, the merger rate evolution shown in population synthesis studies (see e.g. [ 140 ]) is fitted using ΛHBH z={γ= 1.56 , k = 1.94 , zp= 2.12}.(22) The PBH merger rate is dominated by the binaries formed at high redshift, before matter-radiation equality, and the merger rate evolution is predicted to be of the form ψ∝t−34/37 , where t is the age of the Universe at redshift z (see Ref. [ 141 ] for a recent review). This robust prediction results from the properties of binaries at high redshift and GW-driven evolution through Peter’s formula [ 142 , 143 ]. We can fit this relation at low redshift with better than a few percent accuracy by choosing ΛPBH z={γ=−k= 1.17};(23) this yields a simplified single power-law expression, where the parameter zpdrops out. III. HIERARCHICAL BAYESIAN INFERENCE SET-UP In the following, we describe the setup of the analysis that we perform to identify and disentangle the various channels. We perform a hierarchical Bayesian analysis of the GWTC-4.0 catalog [ 144 , 145 ] with icarogw [ 146 ], a Python code developed to infer astrophysical and cosmological population properties of noisy, heterogeneous, and incomplete observations. We assume standard cosmological parameters for the ΛCDM model [ 147 ]. The core of the hierarchical Bayesian analysis is the construction of the merger rate, for which we adopt two models: •Model I (common mass and redshift distributions). In this minimal model, both the mass distribution and the redshift evolution are shared across all channels. Hence, only the spin distributions are allowed to vary depending on the formation channel: dNBBH(Λ) dλ dz dt =dVc dz ψ(z|Λz) 1+zp(m|Λm) ×X c Rc 0pc(χ, cos  θ|m, Λc),(24) where c takes values in the subset of models {IBH,HBH,AGN,PBH}. •Model II (common mass distribution). All the aforementioned channels share a common mass distribution, while both the redshift evolution and spin distributions remain channel-dependent. The overall merger rate is then defined as dNBBH(Λ) dλ dz dt =dVc dz 1 1+zp(m|Λm) ×X c Rc 0ψ(z|Λc z)pc(χ, cos  θ|m, Λc). Here, λ = ( m, χ,  θ )denotes the set of intrinsic binary parameters; z is the redshift; t is the source-frame time; and Λare the model hyperparameters. The term dVc/dz is the differential comoving volume element, and R ( z| Λ) describes the redshift-dependent merger rate density. The function p(m|Λ) models the distribution of source-frame component masses, while p ( χ, cos  θ|m, Λ) specifies the joint distribution between spin magnitudes and spin orientations. Next, the hierarchical likelihood for Nobs GW observations {x}can be written as [148] L({x}|Λ) ∝e−Nexp(Λ) Nobs Y i Tobs ZdλdzLobs (xi|λ, z) ×dNBBH(Λ) dλdzdt , (25) 8 where Lobs (xi|λ, z) is the likelihood of the single GW event xi . The factor Nexp , corresponding to the expected number of detectable events for the model with hyperparameters Λ, encodes the selection effects, and is defined as Nexp(Λ) = Tobs Zdλ dz Pdet(λ, z)dNBBH(Λ) dz dλ dt ,(26) where Tobs is the observation time, and Pdet ( λ, z )denotes the probability of detecting an event with intrinsic parameters λ at redshift z . An event is considered detected if it exceeds the threshold defined by the search pipeline, such as a minimum signal-to-noise ratio (SNR) or a falsealarm rate (FAR) limit [ 144 , 145 ]. Following the LVK state-of-the-art analyses, we impose a cut on the FAR for each event being smaller than (FARmin)−1≥ 1yr, where FARmin is the minimum FAR computed among the search pipelines active at the detection of the specific event. The integrals in the hierarchical likelihood are computed numerically using Monte Carlo (MC) integration of a set of finite samples. To ensure numerical stability and good estimates of these integrals, we use a variance cut [ 149 ]. For more details on the numerical likelihood evaluation, see Appendix B. It has been shown that GW population inference can be biased unless the variance of the log-likelihood estimator is below unity, especially when including spin information in the models [ 150 ]. Therefore, following the recent LVK population analysis [ 2 ], we adopt a threshold of σ2 ln ˆ L = 1 to mitigate potential biases in the posterior. Above this threshold, the likelihood estimate may not be sufficiently converged, and posterior samples with larger variances are therefore discarded. In practice, this can exclude substantial regions of the hyperparameter space for some models, limiting the range of populations that can be robustly explored [ 151 – 158 ]. As shown in Appendix B, the likelihood-variance cut can indeed impact the inference, in particular for models which include spin and/or spiky distributions. When a model fits the data poorly, evaluating the likelihood requires sampling points in the far tail of the posterior of some events, which suffer from poor coverage due to the finite sampling of the GW event posterior, or evaluating the model distribution on regions of parameter space which have few injections needed to compute the selection function (see, e.g., Appendix D3 of [2] and Appendix A of [158]). The population models introduced above depend on a set of hyperparameters that describe the underlying distributions of BBH properties. To carry out hierarchical inference, we impose prior distributions on these hyperparameters, as shown in Table I. The priors adopted for the mass distribution are shown in Appendix C. IV. RESULTS This section presents the results of our inference analysis with the GWTC-4.0 catalog [ 144 , 145 ]. We begin Model Mlog10(BM ⋆)RL max RL av 1 population IBH 5.3 4.7 4.9 HBH 6.2 5.7 5.8 AGN 7.0 6.6 6.9 PBH – – – 2 populations IBH + HBH 7.0 6.2 6.2 IBH + AGN 7.2 6.5 6.7 IBH + PBH 6.5 5.7 5.9 HBH + PBH 7.0 5.6 5.7 3 populations Model I: IBH+HBH+PBH 7.6 6.4 6.4 Model II: IBH+HBH+PBH 7.4 6.4 6.4 TABLE II: Log 10 Bayes factors (second column), ratio of maximum likelihood (third column), and ratio of average likelihood (fourth column), for each combination of the models considered in this work, relative to the Gaussian Component Spins spin model for GWTC-4.0, identified as ⋆ as in [ 2 ] (see Eq. (28) for the definitions of RL max and RL av ). Results shown are derived including the σ2 ln ˆ L cut in each Bayesian inference. by assessing which of the formation channels considered (IBH, HBH, AGN, or PBH) are required to reproduce the observed BBH population, first examining each channel individually, then mixtures of two, and finally three-channel configurations. We then discuss, within these combined IBH+HBH+PBH scenarios, how the different channels can jointly contribute to the observed mass and spin distributions. 1 In the 1and 2-population cases, we use Model I only, in which the merger rate evolution is shared across all channels, while for the 3-population case we also explore Model II (see Sec. III). For all the cases listed above, the reconstructed mass distribution is in agreement with the one reconstructed by the LVK Collaboration in Ref. [ 2 ]. In other words, varying the assumptions on the spin model does not significantly affect the inference of the mass distribution. We discuss this in more detail in Appendix C. 1 The HBH and AGN channels, which differ in our modeling only through the spin orientation, are highly similar and degenerate, as the current dataset does not contain enough information to distinguish between different spin-orientation distributions. For this reason, in the three-channel case we retain the HBH channel as our reference model. 9 A. Model comparison In Table II we report the estimated log 10 Bayes factors for each of the scenarios considered. The reference model is taken to be the Gaussian Component Spins model, denoted by ⋆ and employed in the LVK analysis [ 2 ], which does not include any correlation between the binary component masses and their spins. The Bayes factor values were computed by taking into account of the effective volume explored by the samples [ 159 ]. Thus, they are defined as: BM ⋆=ZM Z⋆·VM eff V⋆ eff (27) being ZM and VM eff respectively the evidence and the effective volume of a given model M. In Table II we also list the ratio between the maximum and posterior-averaged likelihoods found by the model and those of the Gaussian Component Spins model, namely RL max ≡log10 max LM max L⋆ ,RL av ≡log10 ⟨L⟩M ⟨L⟩⋆ .(28) We examine the Bayes factors together with the maximum and average likelihoods to assess whether any model is preferred by the GW data. While Bayes factors are inevitably affected by the choice of priors, reporting information about the likelihood allows us to more robustly interpret how well each model, with a varying number of parameters, improves the inference. As the max-likelihood ratios and the ratios of average likelihoods show similar values to the Bayes factors, we are confident that the effect from the prior volume remains subdominant. When including only one formation channel, the results show that the IBH, HBH, and AGN channels are decisively preferred with respect to the Gaussian Component Spins model, indicating that the catalog strongly favors models featuring mass-spin correlations, in agreement with the conclusions of Ref. [ 53 ]. Among the various models composed of a single population, the AGN channel seems to provide the best fit, with only a slight (but not statistically conclusive) preference over the HBH channel. This small advantage of the AGN scenario likely arises from its preferentially aligned spin directions, which are only marginally favored by the data [ 2 ]. Both AGN and HBH are favored over the IBH scenario, with a relative difference of respectively log10 ( BHBH IBH )=0 . 9, and log10 ( BAGN IBH )=1 . 6. No Bayes factor is reported for the PBH-only scenario, since the analysis does not converge due to the likelihood variance cut, which prevents an adequate exploration of the parameter space and thus reliable fits. This is because the PBH channel struggles to reproduce mildly spinning BHs at low masses without simultaneously over-predicting large spins at higher masses. Explaining the low-mass events, in fact, requires efficient accretion (i.e., low zcut-off ), which in turn predicts large spins across the entire mass range. In short, a single population of PBHs cannot describe the observed spins from GWTC-4.0. When combining two formation channels, both the IBH+AGN and IBH+HBH models perform comparably to the AGN or HBH channels alone. This already suggests that, in these mixed scenarios, the IBH component is subdominant relative to the dynamical channels. This conclusion is further supported by Bayes factors of log10 B = O (2) with respect to the IBH-only case, indicating a strong preference for including at least one dynamical channel. Interestingly, the data favor a spin magnitude distribution characteristic of hierarchical mergers, whether HBH or AGN, while also showing a mild preference for an aligned-spin component. In the IBH+AGN model this component is supplied by the AGN channel, whereas in the IBH+HBH model it is provided by the IBH subpopulation. The 2-population model IBH+PBH, which does not include hierarchical mergers (i.e., either HBH or AGN), is essentially equivalent to the HBH model alone. We also explore configurations in which the PBH channel is added to either the IBH or HBH populations. In both cases, the inclusion of a PBH component increases the log10 Bayes factor by approximately one. We do not consider the HBH+AGN combination, as these channels are degenerate apart from their spin-orientation distributions, and thus this scenario is expected to yield evidence comparable to that of the individual channels. Finally, we consider the 3-channel combination IBH+HBH+PBH, for Model I and Model II , respectively. For Model I , which assumes a common redshift distribution for all channels, the inclusion of an additional population is preferred, but with a low statistical significance. Indeed, the 3-channel model is preferred by merely log10 B ∼ O (0 . 4) with respect to the best 2-population case, IBH+AGN. When relaxing the assumption of a common redshift distribution with Model II , we find only minor differences between that and the 2-population scenario. This outcome simply reflects the fact that (i) information on the merger-rate redshift evolution remains subdominant compared to spin information in the current catalog, and/or (ii) the data do not require each subpopulation to follow radically different redshift distributions. Considering the 2-population and 3-population models, we are interested in what fraction of each possible formation channels contributes to the overall BBH population. The corner plots in Fig. 3display the local merger rate densities for the 2and 3-population scenarios, respectively. Whenever present, both the HBH and AGN channels provide the dominant contribution to the merger rate, with R0∼ 20 Gpc−3yr−1 , a value similar to (and fully consistent with) the total merger rate density inferred in [ 2 ]. Only in the 2-population IBH+PBH case we find that both channels may contribute comparably to the overall merger rate, with the PBH component being subdominant and slightly anti-correlated with the IBH rate. In the 3-population scenario, we also find that the 16 and International Cooperation Grant No. PGR01167. This work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility ( https: //www.arch.jhu.edu/ ), which is supported by the NSF Grant No. OAC-1920103. F.C. acknowledges the financial support provided under the “Progetti per Avvio alla Ricerca Tipo 1,” protocol number AR12419073C0A82B. G.F. thanks IFPU and the organizers of the workshop "Primordial BHs in the Multi-Messenger Era" for the stimulating environment where part of this work was carried out and first presented. S.M. and G.P. are supported by the ERC grant GravitySirens 101163912. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. P.P. is supported by the MUR FIS2 Advanced Grant ET-NOW (CUP: B53C25001080001) and by the INFN TEONGRAV initiative. Some numerical computations were performed at the Vera cluster, supported by MUR and Sapienza University of Rome. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation. Appendix A: The LVK Gaussian component spin model In this section, we present the reference spin models adopted in the population analysis carried out by the LVK Collaboration in the GWTC-4.0 release [2]. Following the most recent analysis by the LVK Collaboration, we model the spin magnitudes χi as a truncated Gaussian distribution between 0 and 1, assuming they are identically and independently distributed. Therefore, p(χ1, χ2|µχ, σχ)=N[0,1](χ1|µχ, σχ)N[0,1](χ2|µχ, σχ). (A1) We model the spin inclination distribution as a mixture between a Gaussian distribution truncated on [ − 1 , 1] and an isotropic distribution, assuming they are identically but not independently distributed [45]: p(cos θ1,cos θ2) = (1 −ζ)×1 4 +ζN[−1,1](cos θ1|µt, σt)N[−1,1](cos θ2|µt, σt). (A2) We report the priors on the Gaussian component spins we adopt in Table III, and the corner plot of the hyperparameters in Fig. 8. TABLE III: Priors on the BBH spin and redshiftevolution hyperparameters for the Gaussian component spin model. Spin model Parameter Prior Default spin model µχU[0,1] σχU[0.005,1] µtU[−1,1] σtU[0.001,4] ζU[0,1] 0.2 0.3 0.4 0.5 σχ −0.8 −0.4 0.0 0.4 0.8 µt 0.8 1.6 2.4 3.2 σt 0.08 0.16 0.24 0.32 µχ 0.2 0.4 0.6 0.8 ζ 0.2 0.3 0.4 0.5 σχ −0.8 −0.4 0.0 0.4 0.8 µt 0.8 1.6 2.4 3.2 σt 0.2 0.4 0.6 0.8 ζ GWTC-4.0 LVK: Gaussian Component Spins This work: Gaussian Component Spins FIG. 8: Posterior distribution for the reference spin models considered in the GWTC-4.0 LVK population analysis. In red, we report the posterior distribution we obtained, while black dashed lines show the LVK posteriors [165]. Appendix B: Impact of effective likelihood variance cuts The hierarchical likelihood is sampled numerically for each population model by combining parameter estimation samples from Nobs GW events with a set of detectable injections used to account for selection effects. Both the injections and the parameter estimation samples provide values of the parameters ( λ, z ), which are then employed to evaluate the BBH merger rate. In addition, the priors πPE and πinj used to generate the parameter estimation samples and the injections, respectively, must be deconvolved to recover the underlying population properties. 17 The overall log-likelihood can be approximated as ln L({x}|Λ) ≈ − Tobs Ngen Nobs X j=1 sj+X i ln   Tobs Ns,i Ns,i X j=1 wi,j , (B1) where sj and wi,j are the weights associated with the injections and the parameter estimation samples, respectively. Here, the index i refers to the i -th GW event, while j labels the Monte Carlo (MC) samples. The weights are given by sj=1 πinj(λj) dNBBH(Λ) dt dz dλ j ,(B2) wi,j =1 πPE(λi,j|Λ) dNBBH(Λ) dt dλ dz i,j .(B3) As we evaluate the MC integrals on a finite number of samples from each event and a finite Ndraw , we must carefully account for the intrinsic variance in the likelihood estimation (see e.g. [ 152 ]). To assess the reliability of our MC estimators for the likelihood, we follow the recent literature and evaluate the variance of the log-likelihood estimator, which varies across parameter space because of the resampling procedures used in Eqs. (B2) and (B3) . By propagating this uncertainty along independent degrees of freedom, the variance of the log-likelihood estimator, σ2 ln ˆ L , for the combined likelihood can be estimated as [152] σ2 ln ˆ L(Λ) = Ndet X i=1 σ2 ˆ Li(Λ) ˆ L2 i(Λ) +N2 detσ2 ξ(Λ),(B4) where σ2 ˆ Li(Λ) = 1 NPE   1 NPE −1 NPE X j=1 w2 i,j −ˆ L2 i(Λ) (B5) is the MC variance in the single-event integrals of (B3) , and σ2 ξ(Λ) = 1 Ndraw   1 Ndraw −1 Nfound X j=1 s2 j−ξ(Λ)2 (B6) is the variance in the detection efficiency MC integrals (B2) . In the previous equations, we introduced NPE , Ndraw , and Nfound , denoting respectively the number of posterior samples available for each detected event, the number of population draws used to estimate selection effects, and the subset of those draws that satisfy the detection criteria. As discussed in Appendix D3 of [ 2 ] and Appendix A of [ 158 ], the inference can be sensitive to the likelihood cuts introduced to ensure reliable likelihood evaluations. We therefore also adopt a second cut based on the effective number of samples. We introduce the effective number of posterior samples per event, defined as [149] Neff,i =PNs,i jwi,j2 PNs,i jw2 i,j ,(B7) which quantifies how many samples per event are contributing to the evaluation of the integral. In our case, we require to have at least an effective number of posterior samples equal to 20 for each event and population model supported by the analysis. In case this requirement is not satisfied, icarogw will artificially associate a null likelihood to the specific point in the parameter space of the population model, as it cannot be trusted. Also, following [ 166 ], we impose numerical stability on the injection by defining the quantity Neff,inj =hPNdet jsji2 PNdet js2 j−N−1 gen PNdet jsj2.(B8) We impose that Neff,inj > 4 Nobs . However, this cut has been pointed out to be less stringent [150]. In Table IV we compare the Bayes factors obtained using the σ2 ln ˆ L (Λ) and Neff cuts. The values for the former correspond to those already shown in Table II, and are reported again here side by side for convenience. Overall, the trends remain similar, although the Bayes factors are slightly larger in some cases. The analysis involving the PBH-only model converges and again confirms that this model alone is ruled out as an explanation of the data. We also observe that, particularly for models featuring preferentially aligned spins—such as IBH, AGN, IBH+HBH, IBH+AGN, and IBH+PBH—the Bayes factors exhibit a substantial increase. While this result cannot be fully trusted, one may speculate that part of this behavior arises because the variance cut can downweight regions of parameter space associated with tight spin alignment. This effect may also be partly linked to the use of a uniform and isotropic spin prior in the GWTC-4.0 parameter estimation, which yields fewer posterior samples near the strongly aligned configuration (see [2] for further discussion). Appendix C: Mass distribution parameters Table Vsummarizes the priors adopted for the parameters governing the primary and secondary BH mass distributions adopted for all models in this work. These include the slopes of the power-law components, the location and width of the Gaussian peaks, and the lower and upper mass cutoffs. The distribution of masses can be further factorized as ppop ≡p ( m1| Λ m ) p ( m2|m1, Λ m ), where Λ m are the parameters of the mass population. The primary mass m1 follows a Power Law + 2 Peaks distribution [ 167 , 168 ]: p(m1|Λm) = (1 −λg,low )P(m1|mmin, mmax,−α) +λgλg,low Gm1|µlow g, σlow g +λg(1 −λg,low )Gm1|µhigh g, σhigh g, (C1) where P is a power law truncated between mmin and mmax with slope α , and G is a Gaussian distribution 18 Model Mlog10(BM ⋆) Neff cut log10(BM ⋆) σ2 ln ˆ Lcut 1 population IBH 5.9 5.3 HBH 6.8 6.2 AGN 8.0 7.0 PBH -12 – 2 populations IBH + HBH 7.6 7.0 IBH + AGN 8.6 7.2 IBH + PBH 8.5 6.5 HBH + PBH 7.6 7.0 3 populations Model I: IBH+HBH+PBH 8.6 7.6 Model II: IBH+HBH+PBH 8.7 7.4 TABLE IV: Log 10 Bayes factors relative to the Gaussian Component Spins spin model for GWTC-4.0, identified as ⋆ as in [ 2 ]. The right column shows the same values reported in Table II, to be compared with the results obtained using the Neff cuts. centered at µg with standard deviation σ . The mixing fractions λg, λg,low control the relative contribution of the two components. A low-mass tapering is introduced through an exponential cutoff governed by a smoothing parameter δm [ 169 ]. The secondary mass m2 is drawn from a truncated power law conditional on the primary: p(m2|m1,Λm)=P(m2|mmin, m1, β),(C2) with slope β and lower bound mmin . Overall, the mass model is characterized by the following parameters: Λm={α, β, mmin, mmax, δm, µg,low, σg,low, µg,high, σg,high, λg, λg,low}.(C3) The parameters and model priors are reported in Table V. In Fig. 9, we show the posterior distribution for the mass hyperparameters, obtained in all analyses performed in this work. We observe only statistically insignificant modifications in the PPD when varying across all the different spin model assumptions. The only noticeable difference arises when removing GW231123 from the dataset. In this case mmax is shifted towards smaller values, down to around 90 M⊙ , consistent with the primary mass of GW190521. Secondly, we find that the PBH-only scenario forces the heavy bump to be more narrowly localized around 55 M⊙ . This behavior is likely related to the shape of the mass-spin correlation, which features a sharp transition from low-mass, low-spin systems to high-mass, high-spin systems. The presence of moderate spins at low masses requires this transition to occur at relatively intermediate/lower masses within TABLE V: Priors on the binary BH mass model hyperparameters adopted for all the analyses in this work. Power-Law + 2 Peaks mass model Parameter Prior range Mass distribution αU[2,5] βU[−1,5] mmin U[3,8] M⊙ mmax U[70,200] M⊙ δmU[0.1,10] M⊙ µg,low U[25,40]M⊙ σg,low U[0.4,10]M⊙ µg,high U[40,100]M⊙ σg,high U[0.4,10]M⊙ λgU[0,1] λg,low U[0,1] the LVK range, corresponding to zcut−off ≲ 20. However, this may lead to a tension with the moderately spinning events in the high-mass portion of the catalog, effectively forcing the preferred masses to be lighter. Overall, we stress that this model alone provides a poor fit to the full catalog (see also [ 170 ] for an analysis using GWTC-3 data). [1] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016),arXiv:1602.03837 [gr-qc]. [2] A. G. Abac et al. (LIGO Scientific, VIRGO, KAGRA), (2025), arXiv:2508.18083 [astro-ph.HE]. [3] I. Mandel and A. Farmer, Phys. Rep. 955, 1 (2022), arXiv:1806.05820 [astro-ph.HE]. [4] M. Mapelli, Proc. Int. Sch. Phys. Fermi 200, 87 (2020), arXiv:1809.09130 [astro-ph.HE]. [5] M. Mapelli, “Formation Channels of Single and Binary Stellar-Mass Black Holes,” (2021) arXiv:2106.00699 [astro-ph.HE]. [6] K. Belczynski, V. Kalogera, and T. Bulik, Astrophys. J. 572, 407 (2001),arXiv:astro-ph/0111452. [7] J. R. Hurley, C. A. Tout, and O. R. Pols, Mon. Not. Roy. Astron. Soc. 329, 897 (2002),arXiv:astro-ph/0201220. [8] M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012),arXiv:1202.4901 [astro-ph.HE]. [9] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, 19 0.0 1.5 3.0 4.5 β 3.2 4.0 4.8 5.6 6.4 mmin 90 120 150 180 mmax 2 4 6 8 δm 27 30 33 36 µg,low 2 4 6 8 σg,low 45 60 75 90 µg,high 2 4 6 8 σg,high 0.04 0.08 0.12 0.16 λg 3.0 3.5 4.0 4.5 α 0.4 0.6 0.8 λg,low 0.0 1.5 3.0 4.5 β 3.2 4.0 4.8 5.6 6.4 mmin 90 120 150 180 mmax 2 4 6 8 δm 27 30 33 36 µg,low 2 4 6 8 σg,low 45 60 75 90 µg,high 2 4 6 8 σg,high 0.04 0.08 0.12 0.16 λg 0.4 0.6 0.8 λg,low 0 20 40 60 80 100 m1[M] 10−5 10−4 10−3 10−2 10−1 p(m1) Posterior predictive distribution Gaussian Component Spins -σ2 ln ˆ Lcut Gaussian Component Spins -Neff cut IBH - σ2 ln ˆ Lcut IBH - Neff cut HBH - σ2 ln ˆ Lcut AGN - σ2 ln ˆ Lcut HBH - Neff cut PBH - Neff cut IBH+HBH - σ2 ln ˆ Lcut IBH+HBH - Neff cut IBH+PBH - σ2 ln ˆ Lcut IBH+PBH - Neff cut IBH+AGN - σ2 ln ˆ Lcut HBH+PBH - σ2 ln ˆ Lcut HBH+PBH - Neff cut Model I: IBH+HBH+PBH - σ2 ln ˆ Lcut Model I: IBH+HBH+PBH - Neff cut Model I: IBH+HBH+PBH (no GW231123) - σ2 ln ˆ Lcut Model I: IBH+HBH+PBH (no GW231123) - Neff cut Model II: IBH+HBH+PBH - σ2 ln ˆ Lcut Model II: IBH+HBH+PBH - Neff cut FIG. 9: Posterior distributions of the mass model hyperparameters inferred from the GWTC-4.0 catalog under different spin model assumptions. We show results including and excluding GW231123, in order to test the impact of this unique event on the inference. Noticeable differences are: (i) excluding GW231123 shifts the preferred mmax to smaller values (∼90 M⊙); (ii) in the PBH-only scenario, the heavy bump is forced to lie more narrowly around 55 M⊙. 20 Astrophys. J. 779, 72 (2013),arXiv:1308.1546 [astroph.HE]. [10] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Astrophys. J. 806, 263 (2015), arXiv:1405.7016 [astro-ph.HE]. [11] S. Stevenson, A. Vigna-Gómez, I. Mandel, J. W. Barrett, C. J. Neijssel, D. Perkins, and S. E. de Mink, Nature Commun. 8, 14906 (2017),arXiv:1704.01352 [astro-ph.HE]. [12] J. Riley et al. (COMPAS Team), Astrophys. J. Supp. 258, 34 (2022),arXiv:2109.10352 [astro-ph.IM]. [13] I. Mandel et al. (COMPAS Team), Astrophys. J. Suppl. 280, 43 (2025),arXiv:2506.02316 [astro-ph.SR]. [14] T. Fragos et al.,Astrophys. J. Suppl. 264, 45 (2023), arXiv:2202.05892 [astro-ph.SR]. [15] J. J. Andrews et al.,Astrophys. J. Suppl. 281, 3 (2025), arXiv:2411.02376 [astro-ph.GA]. [16] N. Ivanova et al.,Astron. Astrophys. Rev. 21, 59 (2013), arXiv:1209.4302 [astro-ph.HE]. [17] M. C. Miller and D. P. Hamilton, Mon. Not. Roy. Astron. Soc. 330, 232 (2002),arXiv:astro-ph/0106188. [18] C. L. Rodriguez, S. Chatterjee, and F. A. Rasio, Phys. Rev. D 93, 084029 (2016),arXiv:1602.02444 [astroph.HE]. [19] I. Bartos, B. Kocsis, Z. Haiman, and S. Márka, Astrophys. J. 835, 165 (2017),arXiv:1602.03831 [astroph.HE]. [20] D. Gerosa and E. Berti, Phys. Rev. D95, 124046 (2017), arXiv:1703.06223 [gr-qc]. [21] M. Fishbach, D. E. Holz, and B. Farr, Astrophys. J. Lett. 840, L24 (2017),arXiv:1703.06869 [astro-ph.HE]. [22] U. N. Di Carlo, N. Giacobbo, M. Mapelli, M. Pasquato, M. Spera, L. Wang, and F. Haardt, Mon. Not. R. Astron. Soc. 487, 2947 (2019),arXiv:1901.00863 [astro-ph.HE]. [23] C. Kimball, C. Talbot, C. P. L. Berry, M. Carney, M. Zevin, E. Thrane, and V. Kalogera, Astrophys. J. 900, 177 (2020),arXiv:2005.00023 [astro-ph.HE]. [24] C. Kimball et al.,Astrophys. J. Lett. 915, L35 (2021), arXiv:2011.05332 [astro-ph.HE]. [25] Y. Bouffanais, M. Mapelli, F. Santoliquido, N. Giacobbo, U. N. Di Carlo, S. Rastello, M. C. Artale, and G. Iorio, Mon. Not. R. Astron. Soc. 507, 5224 (2021), arXiv:2102.12495 [astro-ph.HE]. [26] S. Afroz and S. Mukherjee, Phys. Rev. D 112, 023531 (2025),arXiv:2411.07304 [astro-ph.HE]. [27] S. Afroz and S. Mukherjee, (2025), arXiv:2505.22739 [astro-ph.HE]. [28] D. Gerosa and M. Fishbach, Nat. Astron. 5, 8 (2021), arXiv:2105.03439 [astro-ph.HE]. [29] J. Samsing, Phys. Rev. D 97, 103014 (2018), arXiv:1711.07452 [astro-ph.HE]. [30] S. Stevenson, C. P. L. Berry, and I. Mandel, Mon. Not. R. Astron. Soc. 471, 2801 (2017),arXiv:1703.06873 [astro-ph.HE]. [31] V. Baibhav, D. Gerosa, E. Berti, K. W. K. Wong, T. Helfer, and M. Mould, Phys. Rev. D 102, 043002 (2020),arXiv:2004.00650 [astro-ph.HE]. [32] R. Farmer, M. Renzo, S. E. de Mink, P. Marchant, and S. Justham, (2019), 10.3847/1538-4357/ab518b, arXiv:1910.12874 [astro-ph.SR]. [33] C. Karathanasis, S. Mukherjee, and S. Mastrogiovanni, Mon. Not. R. Astron. Soc. 523, 4539 (2023), arXiv:2204.13495 [astro-ph.CO]. [34] H. Tong et al., (2025), arXiv:2509.04151 [astro-ph.HE]. [35] S. Afroz and S. Mukherjee, (2025), arXiv:2509.09123 [astro-ph.HE]. [36] D. Gerosa, E. Berti, R. O’Shaughnessy, K. Belczynski, M. Kesden, D. Wysocki, and W. Gladysz, Phys. Rev. D98, 084036 (2018),arXiv:1808.02491 [astro-ph.HE]. [37] S. W. Hawking, Nature 248, 30 (1974). [38] B. J. Carr, Astrophys. J. 201, 1 (1975). [39] E. Bagui et al. (LISA Cosmology Working Group), Living Rev. Rel. 28, 1 (2025),arXiv:2310.19857 [astro-ph.CO]. [40] C. Byrnes, G. Franciolini, T. Harada, P. Pani, and M. Sasaki, eds., Primordial Black Holes, Springer Series in Astrophysics and Cosmology (Springer, 2025). [41] S. Vitale, R. Lynch, R. Sturani, and P. Graff, Class. Quantum Grav. 34, 03LT01 (2017),arXiv:1503.04307 [gr-qc]. [42] C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, Astrophys. J. Lett. 832, L2 (2016), arXiv:1609.05916 [astro-ph.HE]. [43] W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, Nature 548, 426 (2017), arXiv:1706.01385 [astro-ph.HE]. [44] B. Farr, D. E. Holz, and W. M. Farr, Astrophys. J. Lett. 854, L9 (2018),arXiv:1709.07896 [astro-ph.HE]. [45] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 011048 (2023),arXiv:2111.03634 [astroph.HE]. [46] G. Franciolini and P. Pani, Phys. Rev. D 105, 123024 (2022),arXiv:2201.13098 [astro-ph.HE]. [47] S. Biscoveanu, T. A. Callister, C.-J. Haster, K. K. Y. Ng, S. Vitale, and W. M. Farr, Astrophys. J. Lett. 932, L19 (2022),arXiv:2204.01578 [astro-ph.HE]. [48] T. A. Callister, S. J. Miller, K. Chatziioannou, and W. M. Farr, Astrophys. J. Lett. 937, L13 (2022), arXiv:2205.08574 [astro-ph.HE]. [49] V. Baibhav, Z. Doctor, and V. Kalogera, Astrophys. J. 946, 50 (2023),arXiv:2212.12113 [astro-ph.HE]. [50] G.-P. Li, Astron. Astrophys. 666, A194 (2022), arXiv:2208.11894 [astro-ph.HE]. [51] C. Périgois, M. Mapelli, F. Santoliquido, Y. Bouffanais, and R. Rufolo, Universe 9, 507 (2023),arXiv:2301.01312 [astro-ph.HE]. [52] J. Heinzel, S. Biscoveanu, and S. Vitale, Phys. Rev. D 109, 103006 (2024),arXiv:2312.00993 [astro-ph.HE]. [53] G. Pierra, S. Mastrogiovanni, and S. Perriès, Astron. Astrophys. 692, A80 (2024),arXiv:2406.01679 [gr-qc]. [54] S. Alvarez-Lopez, J. Heinzel, M. Mould, and S. Vitale, (2025), arXiv:2506.20731 [astro-ph.HE]. [55] H. Tong, T. A. Callister, M. Fishbach, E. Thrane, F. Antonini, S. Stevenson, I. M. Romero-Shaw, and F. Dosopoulou, (2025), arXiv:2511.05316 [astro-ph.HE]. [56] V. Tiwari, (2025), arXiv:2510.25579 [astro-ph.HE]. [57] Y.-Z. Wang, Y.-J. Li, S.-J. Gao, S.-P. Tang, and Y.-Z. Fan, (2025), arXiv:2510.22698 [astro-ph.HE]. [58] W.-H. Guo, Y.-J. Li, Y.-Z. Wang, Y. Shao, S. Wu, T. Zhu, and Y.-Z. Fan, Astrophys. J. 975, 54 (2024), arXiv:2406.03257 [astro-ph.HE]. [59] L. Szemraj and S. Biscoveanu, (2025), arXiv:2507.23663 [gr-qc]. [60] Y.-J. Li, Y.-Z. Wang, S.-P. Tang, and Y.-Z. Fan, (2025), arXiv:2509.23897 [astro-ph.HE]. [61] T. A. Callister and W. M. Farr, Phys. Rev. X 14, 021005 (2024),arXiv:2302.07289 [astro-ph.HE]. [62] J. Golomb and C. Talbot, Phys. Rev. D 108, 103009 21 (2023),arXiv:2210.12287 [astro-ph.HE]. [63] S. Rinaldi, W. Del Pozzo, M. Mapelli, A. LorenzoMedina, and T. Dent, Astron. Astrophys. 684, A204 (2024),arXiv:2310.03074 [astro-ph.HE]. [64] J. Heinzel, M. Mould, S. Álvarez-López, and S. Vitale, Phys. Rev. D 111, 063043 (2025),arXiv:2406.16813 [astro-ph.HE]. [65] J. Heinzel, M. Mould, and S. Vitale, Phys. Rev. D 111, L061305 (2025),arXiv:2406.16844 [astro-ph.HE]. [66] S. Rinaldi, Y. Liang, G. Demasi, M. Mapelli, and W. Del Pozzo, Astron. Astrophys. 702, A52 (2025), arXiv:2506.05929 [astro-ph.HE]. [67] N. Guttman, E. Payne, P. D. Lasky, and E. Thrane, (2025), arXiv:2509.09876 [astro-ph.HE]. [68] C. Adamcewicz, N. Guttman, P. D. Lasky, and E. Thrane, (2025), arXiv:2509.04706 [astro-ph.HE]. [69] O. Sridhar, A. Ray, and V. Kalogera, (2025), arXiv:2511.22093 [astro-ph.HE]. [70] T. A. Callister, (2024), arXiv:2410.19145 [astro-ph.HE]. [71] A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky, and R. Lasenby, Phys. Rev. D 95, 043001 (2017),arXiv:1604.03958 [hep-ph]. [72] K. K. Y. Ng, O. A. Hannuksela, S. Vitale, and T. G. F. Li, Phys. Rev. D 103, 063010 (2021),arXiv:1908.02312 [gr-qc]. [73] K. K. Y. Ng, S. Vitale, O. A. Hannuksela, and T. G. F. Li, Phys. Rev. Lett. 126, 151102 (2021), arXiv:2011.06010 [gr-qc]. [74] P. S. Aswathi, W. E. East, N. Siemonsen, L. Sun, and D. Jones, (2025), arXiv:2507.20979 [gr-qc]. [75] A. Caputo, G. Franciolini, and S. J. Witte, (2025), arXiv:2507.21788 [hep-ph]. [76] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, Astrophys. J. 910, 152 (2021),arXiv:2011.10057 [astro-ph.HE]. [77] K. W. K. Wong, K. Breivik, K. Kremer, and T. Callister, Phys. Rev. D 103, 083021 (2021),arXiv:2011.03564 [astro-ph.HE]. [78] Y. Bouffanais, M. Mapelli, F. Santoliquido, N. Giacobbo, G. Iorio, and G. Costa, Mon. Not. Roy. Astron. Soc. 505, 3873 (2021),arXiv:2010.11220 [astro-ph.HE]. [79] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale, Phys. Rev. D 105, 083526 (2022),arXiv:2105.03349 [gr-qc]. [80] S. Colloms, C. P. L. Berry, J. Veitch, and M. Zevin, Astrophys. J. 988, 189 (2025),arXiv:2503.03819 [astroph.HE]. [81] K. Belczynski, R. E. Taam, E. Rantsiou, and M. van der Sluys, Astrophys. J. 682, 474 (2008),arXiv:astroph/0703131. [82] D. Gerosa, M. Kesden, E. Berti, R. O’Shaughnessy, and U. Sperhake, Phys. Rev. D 87, 104028 (2013), arXiv:1302.4442 [gr-qc]. [83] K. Belczynski et al.,Astron. Astrophys. 636, A104 (2020),arXiv:1706.07053 [astro-ph.HE]. [84] M. Mapelli, Front. Astron. Space Sci. 7, 38 (2020), arXiv:2105.12455 [astro-ph.HE]. [85] N. Steinle and M. Kesden, Phys. Rev. D 103, 063032 (2021),arXiv:2010.00078 [astro-ph.HE]. [86] D. Gangardt, N. Steinle, M. Kesden, D. Gerosa, and E. Stoikos, Phys. Rev. D 103, 124026 (2021), arXiv:2103.03894 [gr-qc]. [87] E. Berti and M. Volonteri, Astrophys. J. 684, 822 (2008), arXiv:0802.0025 [astro-ph]. [88] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO]. [89] J. M. Bardeen, J. Bond, N. Kaiser, and A. Szalay, Astrophys. J. 304, 15 (1986). [90] V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra, and A. Riotto, J. Cosmology Astropart. Phys. 05, 018 (2019),arXiv:1903.01179 [astro-ph.CO]. [91] M. Mirbabayi, A. Gruzinov, and J. Noreña, J. Cosmology Astropart. Phys. 2003, 017 (2020),arXiv:1901.05963 [astro-ph.CO]. [92] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, J. Cosmology Astropart. Phys. 06, 044 (2020), arXiv:2005.05641 [astro-ph.CO]. [93] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, J. Cosmology Astropart. Phys. 04, 052 (2020), arXiv:2003.02778 [astro-ph.CO]. [94] V. De Luca and N. Bellomo, “The Accretion, Emission, Mass and Spin Evolution of Primordial Black Holes,” (2025) arXiv:2312.14097 [astro-ph.CO]. [95] A. G. Abac et al. (LIGO Scientific, VIRGO, KAGRA), Astrophys. J. Lett. 993, L25 (2025),arXiv:2507.08219 [astro-ph.HE]. [96] A. Ray, S. Banagiri, and E. Thrane, (2025), arXiv:2510.07228 [gr-qc]. [97] C. Yuan, Z.-C. Chen, and L. Liu, Phys. Rev. D 112, L081306 (2025),arXiv:2507.15701 [astro-ph.CO]. [98] I. Cuceu, M. A. Bizouard, N. Christensen, and M. Sakellariadou, (2025), arXiv:2507.20778 [gr-qc]. [99] A. Tanikawa, S. Liu, W. Wu, M. S. Fujii, and L. Wang, (2025), arXiv:2508.01135 [astro-ph.SR]. [100] D. Croon, J. Sakstein, and D. Gerosa, (2025), arXiv:2508.10088 [astro-ph.HE]. [101] V. De Luca, G. Franciolini, and A. Riotto, (2025), arXiv:2508.09965 [astro-ph.CO]. [102] S. A. Popa and S. E. de Mink, (2025), arXiv:2509.00154 [astro-ph.HE]. [103] G.-P. Li and X.-L. Fan, (2025), arXiv:2509.08298 [astroph.HE]. [104] L. Paiella, C. Ugolini, M. Spera, M. Branchesi, and M. A. Sedda, (2025), arXiv:2509.10609 [astro-ph.GA]. [105] G. Fabj, C. Tiede, C. Rowan, M. Pessah, and J. Samsing, (2025), arXiv:2510.07952 [astro-ph.HE]. [106] L. Passenger, S. Banagiri, E. Thrane, P. D. Lasky, A. Borchers, M. Fishbach, and C. S. Ye, (2025), arXiv:2510.14363 [astro-ph.HE]. [107] V. Kalogera, Astrophys. J. 541, 319 (2000),arXiv:astroph/9911417. [108] S. Vitale, R. Lynch, J. Veitch, V. Raymond, and R. Sturani, Phys. Rev. Lett. 112, 251101 (2014), arXiv:1403.0129 [gr-qc]. [109] S. S. Bavera, T. Fragos, Y. Qin, E. Zapartas, C. J. Neijssel, I. Mandel, A. Batta, S. M. Gaebel, C. Kimball, and S. Stevenson, Astron. Astrophys. 635, A97 (2020), arXiv:1906.12257 [astro-ph.HE]. [110] P. Hut, Astron. Astrophys. 99, 126 (1981). [111] M. Safarzadeh, W. M. Farr, and E. Ramirez-Ruiz, Astrophys. J. 894, 129 (2020),arXiv:2001.06490 [gr-qc]. [112] D. Gerosa and E. Berti, Phys. Rev. D100, 041301 (2019), arXiv:1906.05295 [astro-ph.HE]. [113] K. Kritos, V. Strokov, V. Baibhav, and E. Berti, Phys. Rev. D 110, 043023 (2024),arXiv:2210.10055 [astro- 22 ph.HE]. [114] K. Kritos, E. Berti, and J. Silk, Phys. Rev. D 108, 083012 (2023),arXiv:2212.06845 [astro-ph.HE]. [115] A. Santini, D. Gerosa, R. Cotesta, and E. Berti, Phys. Rev. D 108, 083033 (2023),arXiv:2308.12998 [astroph.HE]. [116] A. Buonanno, L. E. Kidder, and L. Lehner, Phys. Rev. D77, 026004 (2008),arXiv:0709.3839 [astro-ph]. [117] F. Hofmann, E. Barausse, and L. Rezzolla, Astrophys. J. Lett. 825, L19 (2016),arXiv:1605.01938 [gr-qc]. [118] A. Borchers, C. S. Ye, and M. Fishbach, (2025), arXiv:2503.21278 [astro-ph.HE]. [119] F. Antonini, M. Gieles, and A. Gualandris, Mon. Not. R. Astron. Soc. 486, 5008 (2019),arXiv:1811.03640 [astroph.HE]. [120] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, and C. S. Ye, Phys. Rev. D100, 043027 (2019),arXiv:1906.10260 [astro-ph.HE]. [121] H. E. Cook, B. McKernan, K. E. S. Ford, V. Delfavero, K. Nathaniel, J. Postiglione, S. Ray, E. J. McPike, and R. O’Shaughnessy, Astrophys. J. 993, 163 (2025), arXiv:2411.10590 [astro-ph.HE]. [122] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, Astrophys. J. 178, 347 (1972). [123] J. M. Bardeen and J. A. Petterson, Astrophys. J. Lett. 195, L65 (1975). [124] R. Abbott et al. (LVK), Mon. Not. Roy. Astron. Soc. 524, 5984 (2023), [Erratum: Mon.Not.Roy.Astron.Soc. 526, 6234 (2023)], arXiv:2212.01477 [astro-ph.HE]. [125] A. H. Nitz and Y.-F. Wang, Phys. Rev. D 106, 023024 (2022),arXiv:2202.11024 [astro-ph.HE]. [126] F. Crescimbeni, G. Franciolini, P. Pani, and A. Riotto, Phys. Rev. D 109, 124063 (2024),arXiv:2402.18656 [astro-ph.HE]. [127] J. Golomb, I. Legred, K. Chatziioannou, A. Abac, and T. Dietrich, Phys. Rev. D 110, 063014 (2024), arXiv:2403.07697 [astro-ph.HE]. [128] F. Crescimbeni, G. Franciolini, P. Pani, and M. Vaglio, Phys. Rev. D 111, 083538 (2025),arXiv:2408.14287 [astro-ph.HE]. [129] T. Nakamura et al.,PTEP 2016, 093E01 (2016), arXiv:1607.00897 [astro-ph.HE]. [130] S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. 119, 221104 (2017),arXiv:1708.07380 [astro-ph.CO]. [131] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP 11, 039 (2021),arXiv:2106.13769 [astro-ph.CO]. [132] O. Pujolas, V. Vaskonen, and H. Veermäe, Phys. Rev. D104, 083521 (2021),arXiv:2107.03379 [astro-ph.CO]. [133] K. K. Y. Ng, G. Franciolini, E. Berti, P. Pani, A. Riotto, and S. Vitale, Astrophys. J. Lett. 933, L41 (2022), arXiv:2204.11864 [astro-ph.CO]. [134] K. K. Y. Ng et al.,Phys. Rev. D 107, 024041 (2023), arXiv:2210.03132 [astro-ph.CO]. [135] G. Franciolini, F. Iacovelli, M. Mancarella, M. Maggiore, P. Pani, and A. Riotto, Phys. Rev. D 108, 043506 (2023),arXiv:2304.03160 [gr-qc]. [136] G. Franciolini, R. Cotesta, N. Loutrel, E. Berti, P. Pani, and A. Riotto, Phys. Rev. D 105, 063510 (2022), arXiv:2112.10660 [astro-ph.CO]. [137] P. Madau and M. Dickinson, Annu. Rev. Astron. Astrophys. 52, 415 (2014),arXiv:1403.0007 [astro-ph.CO]. [138] K. K. Y. Ng, S. Vitale, W. M. Farr, and C. L. Rodriguez, Astrophys. J. Lett. 913, L5 (2021),arXiv:2012.09876 [astro-ph.CO]. [139] K. Belczynski, D. E. Holz, T. Bulik, and R. O’Shaughnessy, Nature 534, 512 (2016), arXiv:1602.04531 [astro-ph.HE]. [140] C. L. Rodriguez and A. Loeb, Astrophys. J. Lett. 866, L5 (2018),arXiv:1809.01152 [astro-ph.HE]. [141] M. Raidal, V. Vaskonen, and H. Veermäe, “Formation of Primordial Black Hole Binaries and Their Merger Rates,” in Primordial Black Holes, edited by C. Byrnes, G. Franciolini, T. Harada, P. Pani, and M. Sasaki (2025) arXiv:2404.08416 [astro-ph.CO]. [142] P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963). [143] P. C. Peters, Phys. Rev. 136, B1224 (1964). [144] A. G. Abac et al. (LIGO Scientific, VIRGO, KAGRA), (2025), arXiv:2508.18080 [gr-qc]. [145] A. G. Abac et al. (LIGO Scientific, VIRGO, KAGRA), (2025), arXiv:2508.18082 [gr-qc]. [146] S. Mastrogiovanni, G. Pierra, S. Perriès, D. Laghi, G. Caneva Santoro, A. Ghosh, R. Gray, C. Karathanasis, and K. Leyde, Astron. Astrophys. 682, A167 (2024), arXiv:2305.17973 [astro-ph.CO]. [147] P. A. R. Ade et al. (Planck), Astron. Astrophys. 594, A13 (2016),arXiv:1502.01589 [astro-ph.CO]. [148] I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. R. Astron. Soc. 486, 1086 (2019),arXiv:1809.02063 [physics.data-an]. [149] C. Talbot and J. Golomb, Mon. Not. R. Astron. Soc. 526, 3495 (2023),arXiv:2304.06138 [astro-ph.IM]. [150] J. Heinzel and S. Vitale, (2025), arXiv:2509.07221 [astroph.HE]. [151] D. Wysocki, J. Lange, and R. O’Shaughnessy, Phys. Rev. D 100, 043012 (2019),arXiv:1805.06442 [gr-qc]. [152] R. Essick and W. Farr, (2022), arXiv:2204.00461 [astroph.IM]. [153] Z. Doctor, D. Wysocki, R. O’Shaughnessy, D. E. Holz, and B. Farr, (2019), 10.3847/1538-4357/ab7fac, arXiv:1911.04424 [astro-ph.HE]. [154] V. Delfavero, R. O’Shaughnessy, D. Wysocki, and A. Yelikar, (2021), arXiv:2107.13082 [gr-qc]. [155] J. Golomb and C. Talbot, Astrophys. J. 926, 79 (2022), arXiv:2106.15745 [astro-ph.HE]. [156] M. Mould, C. J. Moore, and D. Gerosa, Phys. Rev. D 109, 063013 (2024),arXiv:2311.12117 [gr-qc]. [157] A. Hussain, M. Isi, and A. Zimmerman, (2024), arXiv:2411.02252 [astro-ph.HE]. [158] M. Mancarella and D. Gerosa, Phys. Rev. D 111, 103012 (2025),arXiv:2502.12156 [gr-qc]. [159] M. Mould, N. E. Wolfe, and S. Vitale, Phys. Rev. D 111, 123049 (2025),arXiv:2504.07197 [astro-ph.IM]. [160] J. Sadiq, T. Dent, and A. Lorenzo-Medina, (2025), arXiv:2506.02250 [astro-ph.HE]. [161] Y.-J. Li, Y.-Z. Wang, S.-P. Tang, and Y.-Z. Fan, Phys. Rev. Lett. 133, 051401 (2024),2303.02973 [astro-ph.HE]. [162] S. Banagiri, T. A. Callister, C. Adamcewicz, Z. Doctor, and V. Kalogera, Astrophys. J. 990, 147 (2025), arXiv:2501.06712 [astro-ph.HE]. [163] A. G. Abac et al. (LIGO Scientific, KAGRA), Astrophys. J. Lett. 993, L21 (2025),arXiv:2510.26931 [astro-ph.HE]. [164] S. Banagiri, E. Thrane, and P. D. Lasky, (2025), arXiv:2509.15646 [astro-ph.HE]. [165] L. S. Collaboration, T. V. Collaboration, and T. K. Collaboration, “Gwtc-4.0: Population properties of merging compact binaries,” (2025). [166] W. M. Farr, Research Notes of the AAS 3, 66 (2019), arXiv:1904.10879 [astro-ph.IM]. 23 [167] C. Talbot and E. Thrane, Astrophys. J. 856, 173 (2018), arXiv:1801.02699 [astro-ph.HE]. [168] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021),arXiv:2010.14533 [astro-ph.HE]. [169] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 041039 (2023),arXiv:2111.03606 [grqc]. [170] G. Franciolini, I. Musco, P. Pani, and A. Urbano, Phys. Rev. D 106, 123526 (2022),arXiv:2209.05959 [astroph.CO].