Potential Subpopulations and Assembling Tendency of the Merging Black Holes

The origins of coalescing binary black holes (BBHs) detected by the advanced LIGO/Virgo are still under debate, and clues may be present in the joint mass-spin distribution of these merger events. Here we construct phenomenological models containing two sub-populations to investigate the BBH population detected in gravitational wave observations. We find that our models can explain the GWTC-3 data rather well, and several constraints to our model are required by the data: first, the maximum mass for the component with a stellar-origin, $m_{\rm max}$, is $39.1^{+2.4}_{-2.7}M_{\odot}$ at 90\% credibility; second, about $15\%$ of the mergers happen in dynamical environments, in which $7-16\%$ of events are hierarchical mergers, and these BHs have an average spin magnitude significantly larger than the first-generation mergers, with ${\rm d}\mu_{\rm a}>0.4 $ at $99\%$ credibility; third, the dynamical component BHs tend to pair with each other with larger total mass and higher mass ratio. An independent analysis focusing on spins is also carried out, and we find that the spin amplitude of component BHs can be divided into two groups according to a division mass $m_{\rm d} = 46.1^{+5.6}_{-5.1}M_{\odot}$. These constraints can be naturally explained by current formation channels, and our results suggest that some of the observed events were likely from AGN disks.


INTRODUCTION
With the help of the Advanced LIGO and Virgo detector network, the number of observed binary black hole (BBH) mergers is rapidly growing. To date, more than 90 BBH merger candidates have been released with the current catalogs of compact binary coalescences (GWTC-3) (Abbott et al. 2019a(Abbott et al. , 2021aThe LIGO Scientific Collaboration et al. 2021a). The origins of these compact objects are still unclear. Various formation channels have been proposed, including, for instance, isolated binary evolution, dynamical capture, and AGN enhancement (see Mapelli 2021 and for recent reviews).
Formation channels may leave their 'fingerprints' in the final population of BHs. One scenario for the origin of BBHs involves the evolution of field binaries. The remnants after the death of very massive stars are expected to fall outside a gap starting at ∼ 40 − 65M and ending at ∼ 125M due to (pulsational) pair-instability supernovae ((P)PISNe) (Belczynski et al. 2016;Woosley 2017;Spera & Mapelli 2017;Stevenson et al. 2019;Mapelli et al. 2020). The BH spins are generally small (Qin et al. 2018;Fuller et al. 2019;Belczynski et al. 2020) under the assumption of efficient angular momentum transport inside their progenitors. On the other hand, in the chain of hierarchical mergers, BHs inside the (P)PISNe mass gap can appear in second or higher-generation mergers, and large spin magnitudes are expected (Miller & Hamilton 2002;Giersz et al. 2015;Gerosa & Berti 2017;Rodriguez et al. 2019;Doctor et al. 2020;Kimball et al. 2020;Arca Sedda et al. 2021;Mapelli et al. 2021;Gerosa & Fishbach 2021). More specifically, for hierarchical mergers happening inside AGN disks, BHs may align their orbits and component spins with the disk, leading to a broad effective spin (χ eff ) distribution with positive mean value (Yang et al. 2019b;Tagawa et al. 2021). Please note that the formation channels discussed above involve conventional concepts, and there are still many theoretical uncertainties that may affect the properties of the remnant BHs (see Sec. 6 for more discussions).
Building phenomenological models to describe the population properties of BBHs is a widely used approach to reveal the key features that could shed light on the formation channels. To capture the main features of the mass distribution for BBHs with a stellar origin, Talbot & Thrane (2018) proposed an astrophysically motivated model consisting of a power-law component and a Gaussian component. A cutoff mass (m max ) for the power law is taken as a parameter to represent the lower edge of the PISN mass gap, and the Gaussian peak is an excess due to PPISNe. While a peak around 35M in the primary mass distribution was revealed by the recent study of The LIGO Scientific Collaboration et al. (2021b) using the latest catalog (see also the in the analysis of some non-parametric approaches (Tiwari & Fairhurst 2021;Li et al. 2021b)), evidence for an upper mass gap was absent. One solution to this discrepancy is that the observed events may originate from multiple channels (Kimball et al. 2021;Galaudage et al. 2021;Wang et al. 2021;Bouffanais et al. 2021;Zevin et al. 2021;Wong et al. 2021;Roulet et al. 2021;Li et al. 2022), thus the mass gap is filled by hierarchical events; finding the exact position of m max remains an important goal in the population studies Wang et al. 2021;Edelman et al. 2021;Baxter et al. 2021). The picture for the spin distribution is even more ambiguous, as it is harder to measure the component spins from the GW signals. In the analysis of GWTC-2 events by Abbott et al. (2021b), the results for the Default spin model (Abbott et al. 2019b) and the Gaussian spin model (Miller et al. 2020) indicated the presence of extremely misaligned spin. However, Roulet et al. (2021) argued that if the model allows for a sub-population with negligible spin, the evidence for negative spins diminishes. Similar results are also obtained by Galaudage et al. (2021), who concluded that a model consisting of a non-spinning sub-population and a rapidly spinning nearly aligned sub-population could better explain the data with respect to LVKC's Default spin model. A recent population study on GWTC-3 by The LIGO Scientific Collaboration et al. (2021b) showed that when the non-spinning sub-population is considered, the data still prefer a negative minimum for the distribution of effective inspiral spins. Very recently, both Callister et al. (2022) and Tong et al. (2022) suggested that there is no evidence for the existence of a non-spinning sub-population, and they found support for the presence of mergers with extreme spin tilt angles. Nevertheless, Vitale et al. (2022) found that the exact shape for the inferred distribution of tilt angle of component BHs strongly depends on the assumed population model and the priors for the model parameters.
The investigations mentioned above may imply that the results obtained from phenomenological models face certain degrees of model dependency due to the finite flexibility of the models. In this work, we carry out the population study from another angle: we assume that the BBHs consist of a "field" sub-population and a "dynamical" subpopulation, and incorporate the pairing functions (Fishbach & Holz 2020a) into the construction of our model. We seek the constraints on the key features and branch ratios of these two components from GWTC-3 events and discuss the requirements on the astrophysical conditions to produce such sub-populations. We also perform a spin-dominant analysis to investigate if there is a correlation between the spins and masses of component BHs to validate the picture of our model. The rest of the paper is arranged as follows: in Sec. 2, we describe our models; in Sec. 3, we introduce the statistical method; the analysis results are presented in Sec. 4; Sec. 5 is our spin-dominant analysis, and we give further astrophysical implication and discussions on our results in Sec. 6 and Sec. 7, respectively.

POPULATION MODELS
One way to model the distribution of the component masses (m 1 , m 2 ) of the merging BBHs is to construct the marginal primary mass distribution and the conditional secondary mass distribution (Abbott et al. 2019b): where Λ m is the hyper-parameters governing the exact shape of the distribution. A form of p(m 2 |m 1 , Λ m ) ∝ m β 2 is generally considered in the literature (Abbott et al. 2019b(Abbott et al. , 2021bThe LIGO Scientific Collaboration et al. 2021b) to reflect the dependency between m 1 and m 2 . On the other hand, Fishbach & Holz (2020a) introduced the 'pairing function' to study the mass distribution. Following their method, the distribution can be written as for m 2 > m 1 .
(2) which can be interpreted as follows: the primary and secondary masses are separately drawn from independent distributions p(m 1 |Λ m1 ) and p(m 2 |Λ m2 ), and the probability of two masses belonging to a merging binary is given by the pairing function w(m 1 , m 2 |Λ p ) (Fishbach & Holz 2020a). We find it convenient to use the pairing function to address diverse preferences for the combination of component masses in multiple channels, so our following modeling is based on Eq.
(2). We consider the distribution model consisting of two sub-populations: one formed through the evolution of field binaries (labeled with 'field'), and the other formed dynamically (labeled with 'dyn'). The model that contains both the mass and spin parameters (Θ = (m 1 , m 2 , a 1 , a 2 , cos θ 1 , cos θ 2 )) can be expressed as where Λ field and Λ dyn are the hyper-parameters for the field and dynamical sub-population respectively, and f dyn is the mixing fraction of the dynamically formed BBHs. Assuming the component BHs in a merging binary are drawn from the same underlying distribution, we extend Eq.
(2) to model each sub-population: where x represents 'field' or 'dyn'. π x is the underlying joint distribution, and we assume that it consists of two independent parts: in which π m,x is the underlying mass distribution and π s,x is the spin distribution for the underlying BHs if they pair to become a component of merging system. For the 'field' channel, π m,field (m i |Λ field ) = P(m i |α, m min , m max , δ m ), and, π s,field (a i , cos θ i |Λ field ) = G(a i |µ a,field , σ a,field )G (cos θ i |µ ct,field , σ ct,field ), where P is a truncated power-law with a smooth tail at low mass as described in Abbott et al. (2019b); G and G are Gaussian distributions truncated at [0, 1] and [−1, 1] respectively. Specifically, under the assumptions of efficient angular momentum transport in stellar evolution (Spruit 2002;Fuller et al. 2019), one would expect µ a,field 1 1 and µ ct,field ∼ 1 (Qin et al. 2018;Fuller et al. 2019;Belczynski et al. 2020) for the binaries formed through the common envelope phase, and the m max should be consistent with the lower edge of PPISN gap (Belczynski et al. 2016). To approximate the preference for symmetric systems in the binary evolution channel, following Fishbach & Holz (2020a), we consider a pairing function of For the dynamical channel, building a phenomenological distribution model is more challenging, as the formation of binaries may happen in different environments, and there might be contributions from hierarchical mergers. To avoid our model being too complicated and redundant, we only consider hierarchical mergers up to the second generation. We first consider a relatively simple case in which the first generation (1G) dynamical BBHs share the same underlying mass distribution with those originated from field binaries, i.e., π m,1G (m i ) = π m,field (m i ). The spin distribution of 1G dynamical BBHs is modeled as π s,1G (m i , a i , cos θ i |Λ field ) = G(a i |µ a,1G , σ a,1G )G (cos θ i |µ ct,1G , σ ct,1G ).
Assuming the component BHs are the remnants of single star evolution, under the efficient angular momentum transport assumption, their spin magnitudes will also be small as the field BBHs (hereafter, we use 'stellar BHs' to represent the field and 1G dynamical BHs); therefore, we let µ a,1G ≡ µ a,field ≡ µ a,ste and σ a,1G ≡ σ a,field ≡ σ a,ste to simplify our model. Guided by the results of some simulations for dynamical processes (O'Leary et al. 2016;Yang et al. 2019a), we consider a pairing function that depends on both the mass ratio and the total mass: In general, dynamical captures in environments such as globular, nuclear, and young star clusters isotropic spin orientation is expected (σ ct,dyn 0) (Mandel & O'Shaughnessy 2010;Rodriguez et al. 2016b;Mapelli et al. 2022), while the evolution of BBHs in AGN disks may lead to certain degrees of alignment, depending on the details of the dynamics in the disk (McKernan et al. 2018(McKernan et al. , 2020Tagawa et al. 2020;Vajpeyi et al. 2022;Mapelli et al. 2022).
The underlying mass distribution of 2G dynamical BBHs can be approximated by the total mass distribution of 1G mergers, which can be calculated by where Λ m,2G = (α, m min , m max , δ m , β dyn , γ dyn ), and we set m * i = 1.05m i to approximate the loss of gravitational mass during the coalescence. The underlying joint distribution for 2G black holes can then be written as Assuming the dynamical environments contain both 1G and 2G BHs, by introducing a parameter r 2G to describe the fraction of 2G underlying BHs, the overall 1G+2G underlying distribution for the dynamical channel is We take the spin tilt hyper-parameters µ ct,field ≡ µ ct,1G ≡ µ ct,2G ≡ 1 in our analysis because a half Gaussian that peaks at 1 and truncates at −1 can approximate both the aligned case (with σ ∼ 0) and the isotropic case (with σ 0). Meanwhile, we also assume that σ ct,1G ≡ σ ct,2G ≡ σ ct,dyn in this work. In addition to the fiducial model (namely the OnePL model) described above, we also consider five additional models for comparison: 1) a Gaussian component representing the pile-up of remnants BHs whose progenitors underwent PPISN is added to the underlying mass distributions ( Talbot

HIERARCHICAL BAYESIAN INFERENCE
We perform hierarchical Bayesian inference to constrain our model parameters. The likelihood for the inference is constructed based on In-homogeneous Poisson process. For a series of measurements of N obs events d, assuming a redshift evolving merger rate of R ∝ (1 + z) 2.7 (as suggested by The LIGO Scientific Collaboration et al. (2021b)), the likelihood for the hyper-parameters Λ can be inferred via (Thrane & Talbot 2019;Abbott et al. 2021b) where N is the expected number of mergers during the observation period and can be derived by integrating the merger rate over the co-moving space-time volume. η(Λ) is the detection efficiency, following the procedures described in Abbott et al. (2021b), we use the injection campaign released in The LIGO Scientific Collaboration et al. (2021b) to estimated this quantity. The n i posterior samples for the i-th event and the default prior π(θ k | ∅) are obtained from the released data accompanying with (Abbott et al. 2019a(Abbott et al. , 2021c; The LIGO Scientific Collaboration et al. 2021a). We use the same criteria that define the detectable events as The LIGO Scientific Collaboration et al. (2021b), i.e., FAR < 1/yr, and 69 BBH events passed the threshold cut 2 . We use the python package Bilby (Ashton et al. 2019a) and the PyMultinest sampler (Buchner 2016) to obtain the Bayesian evidence and posteriors of the hyper-parameters for each model. Note. (i). For all models, we set µ a,field ≡ µa,1G ≡ µa,ste, σ a,field ≡ σa,1G ≡ σa,ste, and σct,1G ≡ σct,2G ≡ σ ct,dyn . We also fix µ ct,field ≡ µct,1G ≡ µct,2G ≡ 1 in the inference; (ii). For the OnePL model, we set α field ≡ α dyn ≡ α; The OnePL & IsoDyn model is the same as OnePL model, except that we fix σ ct,dyn to a very large value to represent the case in which the spin of dynamical BBHs tilt isotropically; (iii). see Appendix. A for the details of PPISN Peak model; (iv). please remind that f dyn and r2G (describing the mixing fractions) are also free parameters in the models; (v). we assume a redshift evolving merger rate of R = R0(1 + z) 2.7 and infer R0 in the analysis.

RESULTS
We summarize in Tab. 3 the (logarithmic) Bayes factors of the models compared to the OnePL model. We interpret ln B > 3.5 (ln B < −3.5) as strong evidence for the corresponding model being more (less) supported by the data than the OnePL model. We also adopt the Akaike information criterion (AIC=2n − 2 ln L max ; Akaike 1981) in the model comparison. The smaller the AIC, the greater the data's preference for the model, and values of ∆AIC > 6 are considered significant for model selection (Liddle 2009). The results show more support for the OnePL model, the TwoPL model, and the PPISN model over the others, and their constraints on the hyper-parameters are shown in Appendix. B.
One intriguing result is that the m max is tightly constrained in the inference. As demonstrated in Fig. 1, the marginal posterior distributions of m max are consistent with each other across the three models. Combing the three posterior sample sets, we obtain an overall constraint of m max = 39.1 +2.4 −2.7 M at 90% credibility. Such tight and consistent constraints may originate from the fact that m max is also constrained by the secondary masses of BBHs, which are also forbidden from the high mass gap. We also investigate the possibility that m max is mainly constrained by the measurement of the heaviest sample (GW190521), whose primary mass is about twice as large as m max . By leaving width of mass range that smoothing function impact on U(0,10) β field power-law index in the pairing function of field BBHs U(0,6) β dyn power-law index for mass ratio in the pairing function of dynamical BBHs U(0,6) γ dyn power-law index for total mass in the pairing function of dynamical BBHs U(-4,12) f dyn the mixing fraction of the dynamically formed BBHs U(0,1) log 10 r2G logarithm of the underlying fraction of the 2nd generation BHs U(-4,0) µa,ste the Gaussian mean value for spin magnitude of stellar BHs U(0,1) σa,ste the Gaussian width for spin magnitude of stellar BHs U(0.05,0.5) µa,2G the Gaussian mean value for spin magnitude of 2nd generation BHs U(0,1) σa,2G the Gaussian width for spin magnitude of 2nd generation BHs U(0.05,0.5) σ ct,field the Gaussian width for spin orientation of the field BBHs U(0.1,4) σ ct,dyn the Gaussian width for spin orientation of the dynamical BBHs U(0.1,4) log 10 r0 logarithm of the local merger rate U(0,3) Note. Here, 'U' means uniform distribution. GW190521 out of the inference, we find that the change on the m max posterior distribution is negligible. Thus, the constraints may come from the overall trend of data rather than a particular event.
To compare with the results in Abbott et al. (2021b), we plot the marginalized primary mass distribution in Fig. 2. Unlike the m 1 distribution shown in Fig. 10 of The LIGO Scientific Collaboration et al. (2021b), the distribution inferred with our model shows a clearly sharp decay in m 1 ∼ 40M . The shallow bump extended to a higher mass is contributed by the 2G dynamical BHs. In some analyses on the GWTC-2 events, a similar sharp decay in the m 1 spectrum is also obtained using different approaches (Wang et al. 2021;Wong et al. 2021;Baxter et al. 2021); now with GWTC-3, we confirm that our model with such a feature can better fit the data with respect to the PP & Default model. To avoid model misspecification (Romero-Shaw et al. 2022b), we perform posterior predicted checks following the procedures described in Abbott et al. (2021b). As shown in Fig. 3, the observed accumulated distribution of primary mass and mass ratio are within the predicted region of the model.
The inference also reveals a significant difference between the spin magnitude of the 2G dynamical BHs and stellarorigin BHs. For the OnePL model, we inferred that µ a,ste = 0.11 +0.08 −0.09 and µ a,2G = 0.87 +0.12 −0.24 ; these two parameters have consistent constraints in other models. We plot the difference between µ a,ste and µ a,2G in Fig. 1. Accounting for all models we have considered, the constraint on the difference is dµ a > 0.4 at 99% credibility. We note that for µ a,ste , the posterior support does not go down to zero at µ a,ste = 0, while for µ a,2G , having an average spin of ∼ 0.7 (which is the case for the remnant of a coalescing system with non-spinning equal-mass components) is within its 90% credible interval. We plot the 90% credible region of the spin amplitude distribution in green for the field + 1G 'dyn' mergers and in purple for the 2G 'dyn' mergers in Fig. 4. We also demonstrate in the figure the reconstructed distributions  indicates that the isotropically tilted case for the dynamical BHs is also strongly disfavored. The posterior predicted check for the effective inspiral spin and the effective precessing spin are demonstrated in Fig. 3.
We address some caveats in interpreting the spin results: first, some astrophysical models predict more complex spin distributions (see Sec. 6. for more discussions), and hence the simple truncated Gaussian has a limitation on recovering the exact shape of the distributions; however, our empirical model can still reflect the rough central tendency of data. Second, we have utilized parameter estimation samples derived with the default prior (Uniform distributions for a i and cos θ i ) for individual events. Thus our inference has limitations in identifying sub-populations with negligible spins, as demonstrated in Galaudage et al. (2021)) (see however Callister et al. (2022) and Tong et al. (2022) for the arguments about this sub-population). We leave further investigation on this issue in future work since it will not change our current conclusion.
The existence of a distinct underlying mass distribution for the dynamical BHs is ambiguous. The logarithmic Bayes factor of the TwoPL model compared to the OnePL model is −1.3, indicating the introduction of the additional parameter, α dyn , does not significantly improve the model in describing the data. The posterior distribution of α dyn implies the underlying mass spectrum of this sub-population may be much flatter than that of the field BBHs, though their 90% credible intervals overlap with each other in the range of 3.17 − 4.16.
The results also reveal different pairing functions between the 'field' BHs and 'dyn' BHs. For the OnePL model, we infer γ dyn = 8.41 +2.15 −2.60 ; for the TwoPL model, the constraint is much weaker, and γ dyn strongly degenerates with α dyn and log 10 r 2G in the inference. We show the corner plot for these parameters in Fig. 5 to illustrate the degeneracy. The parameter β dyn is constrained to be 3.56 +2.41 −1.71 (3.54 +3.39 −1.89 ) for the OnePL (TwoPL) model; for the pairing function of field binaries, the posterior distributions of β field skew toward small values. This finding is consistent with the results obtained by Li et al. (2022) that suggests the events with larger masses (likely dominated by dynamical mergers) prefer more symmetric systems. As a conclusion for the pairing function, one requirement is clear for our model to fit the data: the underlying dynamical BHs tend to pair together with more equal masses and larger total mass.
We further use the results to derive branch ratios for different types of mergers. Combing the posterior distributions from OnePL, TwoPL, and PPISN models, we obtain f dyn = 0.15 +0.12 −0.6 , which means ∼ 15% of the astrophysical mergers would form through the dynamical process. The parameter log 10 r 2G describes the fraction of underlying 2G BHs in the dynamical environments. Its posterior distribution for the TwoPL model deviates from the results of the OnePL model and the PPISN model due to its degeneracy with γ M and α dyn . We further translate the constraint on r 2G into the fraction of mergers that contain at least one 2G BH in the dynamical sup-population, and the result is shown in Fig. 5. Combing the constraints from different models, the fraction is 0.10 +0.06 −0.03 at 90% credible level. After accounting for the selection effect in the observation, our result suggests that such events make up 16 +8 −7 % of the detectable mergers, and 6 +6 −4 % of the systems have both components being 2G BHs. We also compute the probability of the formation channels for each event according to the constrained OnePL model (details are presented in Appendix. C), and find that 11 out of the 69 events are most likely (probability > 50%) containing 2G BHs. In particular, GW190521 030229 and GW191109 010717 have probabilities of ∼ 1 and 0.74 to be 2G+2G BBHs.

THE SPIN-DOMINANT ANALYSIS
We have shown in Sec.4 that the current data can be explained by our model with two sub-populations. Two important constraints on our model were derived: the m max ∼ 40M and the significant difference between the spin magnitudes of the two sub-populations. To investigate the model dependency of these results, here we perform the "spin-dominant analysis". The model used in the analysis has consistent astrophysical motivation with the more complex models in Sec.2, while its parameters are free from the assumptions of the mass model.
We assume that the BHs with masses above a certain division (m d ) have a different spin distribution from that of BHs below m d , then our model describing an individual BH can be written as: the parameters included in γ s are summarized in Tab. 5 along with their descriptions and priors. Note that both BHs in a system follow this model, so the likelihood for a particular posterior sample of an event is π s (a 1 , cos θ 1 , m 1 |Λ s ) × π s (a 2 , cos θ 2 , m 2 |Λ s ). Though different sub-populations may have overlapping BH masses, the model can represent a switch of the dominant sub-population at m d .
The inferred hyper-parameters for this model are shown in Fig. 6. We find that the constraints for (µ a,1 , σ a,1 , σ ct,1 , µ a,2 , σ a,2 , σ ct,2 ) are consistent with those for (µ a,ste , σ a,ste , σ ct,field , µ a,2G , σ a,2G , σ ct,dyn ) in the OnePL model. We derive a division mass of m d = 46.11 +5.59 −5.11 M , which is in agreement with the inferred m max but with a slightly larger median. It is worth noting that neglecting the spin data has a very small impact on the constraints on m max , so we conclude that both mass and spin data show evidence that the dominant sub-population may change at a mass of ∼ 36 − 52M .  Figure 5. Posterior distribution the parameters about the dynamical channel for the OnePL and the TwoPL models; note that r2G is the fraction of 2G BHs in the dynamical environment, and f2G is the fraction of BBHs involving a least one 2G BH.  Our phenomenological models provide a possible way to fit the observed BBH population. Here we discuss the question: what astrophysical conditions are required to produce the resulting sub-populations (described by the constrained hyper-parameters)? This is a self-consistency check for our model; alternatively, assuming our model is correct, this might provide insights into some key aspects regarding the formation of BBHs.

ASTROPHYSICAL IMPLICATION
In principle, the position of the maximum BH mass can be explained by a combination of PI physics at low metallicities (Z) and Fe-dependent Wolf-Rayet (WR) mass loss at high Z without needing to change the uncertain nuclear reaction rate. As shown in Farmer et al. (2019Farmer et al. ( , 2020, a M max ∼ 40M can be attained by varying the 12 C(α, γ) 16 O thermonuclear reaction rate within its 1σ uncertainties for the progenitors with Z ≤ 3 × 10 −3 . On the other hand, at higher Z, iron (Fe) dependent WR winds (Vink & de Koter 2005;Sander & Vink 2020) could be the leading factor that decides M max . When the recent WR theoretical mass-loss recipe of Sander & Vink (2020) was employed in stellar models by Higgins et al. (2021), it was found (Fig. 13 in that paper) that above Z/Z = 0.5 the maximum final mass was of order 40 M , and close to the maximum BH mass uncovered by our phenomenological model.
When a more complex function of the (maximum) BH masses becomes available as a function of Z and redshift z, the issue of constraining the 12 C(α, γ) 16 O thermonuclear reaction rate could be revisited. In addition, the sub-population of very hot WR stars and the WO stars could help constrain this elusive nuclear reaction rate (Tramper et al. 2015). In other words, multi-messenger astrophysics, including both GW events and electromagnetic (EM) information, could be combined to determine both the 12 C(α, γ) 16 O rate as well as stellar wind physics, by disentangling high Z and low Z GW events.
The BH spins have been widely considered to distinguish the BBH formation channels. In the isolated evolution of massive stars via common-envelope phase (Ivanova et al. 2013), the first-born BH originating from an initially massive star has negligible/small spin due to the standard angular momentum transport (Qin et al. 2018;Fuller et al. 2019;Belczynski et al. 2020); the second-born BH from the less massive star can have spin from zero to maximally spinning due to tidal spin-up (Qin et al. 2018;Hu et al. 2022), while this effect is limited within the horizon of current gravitational-wave detectors (Bavera et al. 2020(Bavera et al. , 2022a. These predictions seem to be in line with our inferred result which allows for the 'field' sub-population to have µ a,ste ∼ σ a,ste ∼ 0.1. Nevertheless, within their 90% credible region, µ a,ste and σ a,ste have more posterior supports at slightly larger values. Higher spins in isolated binary evolution can be attained by various processes, like the tidal spin-up of two He stars following the double-core common envelope phases (Olejak & Belczynski 2021), the Eddington-limited accretion onto BHs, the efficient angular momentum transport within massive stars (Qin et al. 2022a;Zevin & Bavera 2022), and the chemically homogeneous evolution Marchant et al. 2016). In addition, both the OnePL/TwoPL/PPISN Peak model and the spin-dominant analysis show that the spin-tilted angle of "field" BBHs are not perfectly aligned. The spins may be tossed due to the supernova kicks during their formation process in the core collapse of massive stars (Tauris 2022).
We have inferred a total local (z = 0) BBH merger rate of 18.2 +7.65 −5.4 Gpc −1 yr −1 , consisting with the result in The LIGO Scientific Collaboration et al. (2021b). Our models enable us to derive branch ratios for different sub-populations. The analysis has shown that the dynamical channel contributes to 9% − 27% of the mergers, in agreement with the expectations of some dynamical environments (Gerosa & Fishbach 2021). The investigation of spin distributions in this work sheds light on the details of dynamical environments. Dynamical assembly in dense stellar environments like globular clusters predict isotropic spins, while the features for the mergers in the AGN disks can change in a variety of ranges depending on the properties of the disk (Vajpeyi et al. 2022) and the dynamics in play (Tagawa et al. 2020). Though at this stage, we might not be able to clearly resolve the tilted angle distribution from the available data (Tong et al. 2022;Vitale et al. 2022), our analysis disfavors both the perfectly alignment case (σ ct,dyn 1) and the isotropically tilted case (σ ct,dyn 1) of the spin, implying that not all of the dynamical mergers happen in old, dense AGNs, or, young, dilute AGNs/globular clusters (Benacquista & Downing 2013;Vajpeyi et al. 2022). It is possible that multiple dynamical environments are contributing to the observable events, and our results suggest that some of the events are likely from AGN disks. The reconstructed spin amplitude distribution for the 2G BBHs has µ a,2G 0, which can be explained by the inheritance of orbital angular momentum of 1G mergers  in the framework of our model. For a coalescing system with non-spinning equal-mass components, the remnant will have spin amplitude ∼ 0.7; On the one hand, non-zero and isotropically distributed spins of 1G BHs (due to, e.g., natal kicks) will increase the dispersion of spin amplitude distribution of 2G mergers; on the other hand, if the component spins and orbital angular momentum were aligned by the accretion torque and co-rotation torques in the AGN disk, the remnants could have larger spins. Recently, Fishbach et al. (2022) addressed the limits on hierarchical black hole mergers from the most negative χ eff systems; their work mainly focus on classical star clusters, further studies can be carried out to quantitatively constraint the mergers in an AGN disk.
Clues for the Hierarchical sub-populations was suggested by Kimball et al. (2021) using the GWTC-2 data. With the GWTC-3 events and by employing the pairing function in the model, we have inferred sub-populations with distributions different from their work. The pairing function for the dynamical sub-population prefers systems with heavier total mass, though large systematical differences are presented in the posterior distributions of γ dyn between the OnePL and the TwoPL model due to its degeneracy with α dyn in the inference. O'Leary et al. (2016) found that dynamical interactions can enhance the merger rate of BBHs by a boost factor ∝ M γ tot , with γ ≥ 4, which is consistent with the constraints on γ dyn in our model. It is interesting to point out Yang et al. (2019a) addressed the AGN disk will also harden the mass distribution of BBHs, and the power-law index of the mass spectrum will change by ∆α dyn ∼ 1.3. This possibility is included in the results for the TwoPL model, considering the statistical uncertainties. The constraint on β dyn is in agreement with the conclusions from some N-body simulations that the pairing of BHs of comparable masses is energetically favorable (Rodriguez et al. 2016a;Banerjee 2017). Further studies comparing the theoretical predictions with the constraints on γ dyn and β dyn could help better understand the dynamic interactions in dense environments.
Among the 11 events (shown in bold in Tab. 6) plausibly contain 2G BHs, the possible hierarchical origins for GW190519 153544, GW190521 030229, GW190602 175927 and GW190706 222641 have also been suggested in Kimball et al. (2021). GW191109 010717, which has a probability of 0.74 to contain two 2G BHs, has non-negligible orbital eccentricity and precession (Romero-Shaw et al. 2022a;Henshaw et al. 2022). Additionally, evidence for the precession in GW190929 012149 was reported by Hoy et al. (2022). Nevertheless, Vigna-Gómez et al. (2021) proposed that the progenitor of GW170729 185629 is a low-metallicity field triple. Further studies on the signals or properties of these individual events may provide additional clues for their origins.

SUMMARY AND DISCUSSION
In this work, we study the joint mass-spin distribution of coalescing binary black holes by incorporating the pairing functions into the hierarchical Bayesian inference. With astrophysical concerns, we propose a model consisting of two sub-populations having different pairing functions and representing binaries form in the filed and dynamical environments, respectively. We find that this model can well explain the current population of BBHs observed by LIGO/Virgo/KAGRA. Under the framework of our model, some key features of different formation channels are revealed and well constrained in the analysis. For the field binaries and 1G dynamical binaries, we obtain a tight constraint on the maximum mass of m max = 39.1 +2.4 −2.7 M ; we also infer that the 2G dynamical mergers have an average spin magnitude that is significantly larger than other kinds of mergers, with dµ a > 0.4 at 99% credibility. As we have discussed in Sec. 6, the inferred values of hyper-parameters in our phenomenological model can be naturally explained by current astrophysical theories. Revealing the potential features and assembling tendency in the mass and spin distributions can help distinguish the acting mechanisms during the formation of field BBHs, as well as the environments of dynamical BBHs.
Our model in this work is designed to capture the main features in the distributions of field binaries and dynamical binaries if they exist; it could be further improved and extended in the future to identify sub-classes under each sub-population. Very massive BHs lying in the conventional PI gap can also be formed through the preserve of large H envelope for progenitor stars with reduced metallicity (Z < 0.1 Z or below) , and could be responsible for some events like GW 190521. Field binaries may also yield systems with relatively high spins and large masses (below the mass gap) through chemically homogeneous evolution Marchant et al. 2016;du Buisson et al. 2020;Bavera et al. 2022b;Qin et al. 2022b). Further studies with additional considerations on these channels are needed to clarify their existence and branch ratios in the GW events.
With more events from future observations, the redshift evolution of mass function and spin distributions might be well modeled, and there are already some studies Biscoveanu et al. 2022) attempting to reveal the trends with current data. Features in the mass spectrum of compact binaries via GW observations, especially the high-mass gap (Fishbach & Holz 2020b) and the low-mass gap (Li et al. 2021a;Ye & Fishbach 2022), are great ingredients in the "spectrum siren" for cosmological probe (Ezquiaga & Holz 2022). Therefore, extending our method to study the redshift evolution of the maximum mass for the BHs may bring new insight into the constraints on the cosmic expansion history (Farr et al. 2019).

A. MODIFICATIONS OF THE POPULATION MODEL
We also try to take some widely concerned astrophysical motivations into our population models. Firstly, we consider the pile-up at the high-mass cutoff of the BH mass spectrum that results from the PPISNe (labeled 'PPISN Peak') (Talbot & Thrane 2018), where P is the normalized power law with smoothing treatment on the low-mass cutoff, so r represents the fraction of the underlying BHs that could have been in the range of ∆ m above m max , but appeared below m max because of the PPISN. Note that ∆ m should be smaller than the width of the high-mass gap caused by the (P)PISN. Therefore, we set the prior of ∆ m to be uniformly distributed in (0, 100)M , ∆ m = 100M , which is wider than the widely expected mass gap (Spera & Mapelli 2017). As presented in Fig. 7, the contribution of the PPISN channel is negligible, and the high-mass cutoff is also tightly constrained, as shown in Fig.1. However, the ∆ m is not well constrained (i.e., is similar to the prior distribution, and ∆ m = 0M (∆ m = 100M ) can not be ruled out). The BHs formed through different formation channels/environments may follow different spectrum indexes. Therefore, we introduce the second modified model (labeled 'TwoPL'), in which we consider another power-law index α dyn for the mass function of first-generation BHs in the dynamical assembly. Finally, we consider the spin orientation of BBHs formed by the dynamical capture may have isotropic distribution; therefore, in the modified model (labeled "IsoDyn"), we assume that π dyn (cos θ 1,2 ) = 1/2 for cos θ 1,2 ∈ (−1, 1), i.e., cos θ 1,2 uniformly distributed in (-1,1). However, from the results summarized in Tab.3, such a scenario is not favored by the GW observation.

B. FULL CORNER PLOT FOR THE POSTERIOR DISTRIBUTIONS
We compare the inferred parameters obtained by the OnePL model and the TwoPL model in Fig. 8. The PPISN Peak model has similar constraints for the overlapping parameters with the OnePL model. We find that most parameters are nearly identical in both models except for γ dyn , log 10 r 2G , which are strongly degenerate with α dyn . Though the constraint on log 10 r 2G is weak, as shown in Fig. 5, the fraction of the mergers that involve at least one 2nd-generation BH is well constrained.

C. PROBABILITIES OF FORMATION CHANNELS FOR EACH EVENT
Using the constraints on the OnePL model, we calculate the probabilities for the formation channels (field, dynamical 1G+1G, 1G+2G 3 , 2G+1G, and 2G+2G) of each event. With the inferred posterior distribution p(Λ| d), the probability of the i-th event came from K channel is where r K is the mixing fraction of the K channel, Z(d i |Λ, K) is the evidence for K channel, and Z(d i |Λ) is the evidence for the combination of all channels, which read Z(d i |Λ) = dθL(d i |θ)π(θ|Λ), where π(θ|Λ, K) is the joint mass-spin distribution for the K channel under the hyper-parameter Λ, and π(θ|Λ) = K r K π(θ|Λ, K). Therefore, one can find that K P i K = 1, i.e., P i K is normalized. We calculate the probabilities of formation channels for the 69 events and summarize them in Tab. 6 (the events that are most likely containing 2G BHs are highlighted in bold face).