1 Introduction↩︎

Since the first direct detection of the GW [1], the LIGO-Virgo-KAGRA (LVK) network has observed \(\sim\)​300 candidates of compact binary coalescences by their third observing run [2]–[5], most of them considered to be mergers of BBH. These observations provide a brand-new approach to inspecting the origin of stellar-mass black holes.

The masses of black holes of the first gravitational wave event, GW150914, were about \(\sim30M_\odot\), and the existence of such BHs was a surprise since there were few predictions of such a massive BH by that time. This observation sparked a discussion about the origin of BBH.

There are two main channels proposed for the formation of merging BBHs: the isolated formation channel and the dynamical formation channel. In the isolated formation channel, the progenitors of binary systems are gravitationally bound from the zero-age main-sequence phase, and they evolve without external disturbances and eventually become merging compact objects. Binary interactions are vital for coalescence within the Hubble time. The channel is finely divided into sub-channels, depending on how the binaries interact. In the classical picture of the channel, a binary undergoes an unstable mass transfer and a common envelope phase, in which the binary drastically loses its orbital energy (for example, [6]–[8]). Recently, however, it was reported that the merging BBHs can also be produced through stable mass transfer (for example, [9], [10]).

The stellar metallicity is one of the elements that affects the stellar evolutionary tracks, which in turn affects the BHs that remain as stellar remnants. The possibility that the isolated binaries of PopIII stars are the origins of the BBHs has been a subject of much discussion (for example, [11], [12]). PopIII stars are the first stars in the universe, with zero metal (or extremely low-metal). The evolution of these stars differs from that of PopI stars or PopII stars. PopI stars are relatively metal-rich stars in the Galactic disk, typically with metallicities ranging from approximately 1/10th to three times that of the Sun [13], [14], while PopII stars are the low metal stars with 0<Z\(\lesssim 0.1Z_\odot\), where \(Z_\odot\) denotes the solar metallicity. The threshold separating PopII and PopIII was reported as \(Z\sim 10^{-5}Z_\odot\) [15]. As the progenitor of BBH, BHs from PopIII stars tend to have higher masses than those from PopI/PopII stars on average, typically \(\sim30M_\odot\). There are reasons for this tendency: first, the initial mass function of PopIII stars is expected to be more top-heavy than that of PopI and PopII stars, due to the inefficient cooling of the primordial gas [16]. Also, binaries of PopIII tend not to lose stellar mass through binary interactions because they tend to evolve via blue giant stars with a small radius [11], [17]–[19]. Moreover, the mass loss by stellar winds is greatly suppressed for no metal stars [20].

On the other hand, several models consider dynamical formation in dense stellar environments. Binary black holes can form efficiently in globular clusters [21]–[23], in young stellar clusters [24]–[26], and in nuclear star clusters [27], [28]. In these environments, repeated dynamical interactions can assemble and modify BBHs prior to merger.

Other formation channels, such as the CHE [29]–[31], dynamical triples [32], [33], the formation in active galactic nuclei [34]–[37], and primordial black holes (PBHs) [38]–[40] have also been proposed. The merger rate predictions of BBHs for various formation channels are summarized in [41].

The BBH formation channels must be consistent with observational results. The merger rate of BBHs, and the distribution of parameters of BBHs inferred from the observation, are key elements to constrain the BBH formation channels. LVK found the overdensities in the mass distribution near \(\sim 10M_\odot\) and \(\sim 35M_\odot\), the population-level anti-correlation between mass ratio and effective spin parameter which was pointed out first in [42], and the redshift evolution of the merger rate density. Those properties need to be explained by formation channels.

Identifying the formation channels of the GW sources is one of the central goals of GW astronomy. A number of studies have attempted to constrain these channels using the growing GW catalogue. On the theoretically-driven approach, population-synthesis models and/or dynamical formation scenarios have been compared against the observed distributions of source parameters. For instance, Ref. [43] combined five state-of-the-art BBH formation pathways and concluded that multiple formation pathways and proper physical prescriptions are needed to interpret the catalogue. In Ref. [44], the PBH channel is explored with a focus on the correlation in the mass-ratio parameter \[q = \frac{m_2}{m_1} (\leq1),\] and the effective inspiral spin parameter of the BBH \[\chi_\text{eff}= \dfrac{(\boldsymbol{\chi}_1+q\boldsymbol{\chi}_2)\cdot\hat{\boldsymbol{L}}}{1+q},\] where \(m_1\) and \(m_2\) are the component masses of the BBH with \(m_1\geq m_2\), \(\boldsymbol{\chi}_1\) and \(\boldsymbol{\chi}_2\) are the dimensionless spin parameters of each BH, and \(\hat{\boldsymbol{L}}\) is the unit vector parallel to the orbital angular momentum. In parallel, data-driven approaches that rely on minimal astrophysical assumptions on the source-parameter distribution shapes, such as [45]. The examples cited above represent only a fraction of the rapidly expanding literature on this topic. For a broader review, see, e.g., [46], [47].

In the GW astronomy context, PopIII stars have attracted interest since the first direct GW observation. The mass range of the first event, GW150914, is both \(\sim 30 M_\odot\) [1]. However, the absence of direct observation of these stars makes the theory uncertain. This uncertainty results in uncertainty in the star formation rate of these stars. Estimated local merger density of BBHs from PopIII stars spans from \(10^{-2}\;\mathrm{Gpc^{-3}yr^{-1}}\) [48] to \(10^{2}\;\mathrm{Gpc^{-3}yr^{-1}}\) [19]. In addition, since the PopIII stars are the first stars in the universe, the merger rate of PopIII BBHs starts to reach a maximum at a high redshift, for which only next-generation gravitational wave detectors can observe. Therefore, there are proposals [11], [19], [49], [50] to verify the channel by using the next-generation detectors such as Einstein Telescope [51].

Although there is large uncertainty in the prediction of the merger rate of the PopIII BBHs in various computations [11], [50], binary population synthesis studies have shown that the BH mass distribution from PopIII binaries typically exhibits a peak around \(\sim 30M_\odot\). In addition, more massive BBH mergers, such as GW190521, can also be interpreted as originating from PopIII BBHs [52], [53]. PopIII stars are difficult to observe directly with optical and infrared telescopes, as they formed and likely ended their lives at very high redshift. Thus, it is important to constrain their properties and the formation history with GW observations.

In this work, we consider a binary population model consisting of a mixture of five BBH population models: isolated binaries from either PopI or PopII, isolated binaries from PopIII, and BBHs dynamical evolution either in GC or NSC, and CHE. We perform a hierarchical Bayesian analysis by using the BBH events listed in the GWTC-3 catalog by LVK. We then constrain the fraction of the contribution of each formation channel to the cosmic merger rate and the observed merger rate of BBHs. We also investigate the distribution of parameters that describe the BBH waveform.

The paper is organized as follows. Sec.2 describes the formation models we adopt in this work and the summary of the hierarchical Bayesian inference. In Sec.3, we describe the result of our inference, namely the posterior distribution for branching fraction and the Bayes factors, in Sec.4, we discuss further the predicted parameter distribution and the constraint on the merger rates. In Sec.5, we summarize our findings.

2 Methods↩︎

In this section, we summarize the formation channel models we adopt and the hierarchical Bayesian analysis framework.

2.1 BBH channels↩︎

This study utilizes sample sets of BBH systems obtained from population syntheses from several published formation models:

• Isolated evolution of PopI and PopII stars

• Isolated evolution of PopIII stars

• Chemically homogeneous evolution

• Dynamical evolution in the globular cluster

• Dynamical evolution in Nuclear stellar cluster

For the isolated evolution of PopI and PopII stars (hereafter Pop. I+II channel), we use the result from [54]. Among their models, we use the model called M30.B, which is the standard model in their work. This model employed state-of-the-art physical models such as star formation rates [55] and an accretion model. The model assumes the efficient angular momentum transport that includes the Tayler-Spruit magnetic dynamo [56] so that the natal spin magnitude for BH is in the range \(0.05\lesssim a\lesssim 0.15\). The BBH local merger rate of the model is 43.7 \(\mathrm{Gpc^{-3} yr^{-1}}\).

For the isolated evolution of PopIII stars (hereafter PopIII channel), we use a model presented in [19]. The model we adopt is the ‘M100’ model in [19], which was shown to fit well with the GWTC-2 results in [57]. The model assumes a flat initial mass function from \(10M_\odot\) to \(100M_\odot\) for Pop III binaries. The local merger rate in this model is 6.36 \(\mathrm{Gpc^{-3} yr^{-1}}\).

For the remaining formation channels, the results of population syntheses are taken from the previous work [43]. More particularly, [31] for CHE, [58] for dynamical evolution in GC, and [28] for dynamical evolution in NSC. All the model assumes the natal spin of the BH is zero. We will denote them as CHE, GC, NSC channel resepectively.

2.2 Population Inference↩︎

In this work, we assume that BBHs in the universe are formed in several different formation channels. We then consider a mixed model in which the total merger rate is the sum of the merger rates of each channel weighted by the branching fractions \(\{f_j\}\) with \(\sum_{j} f_j = 1\), where \(j\) denotes each formation channel. Here, \(f_j\) represents the underlying contribution of the \(j\)-th channel to the astrophysical BBH population, corrected for observational selection effects. We apply a hierarchical Bayesian method [59], [60] to estimate the branching fractions by using the BBH events in the GWTC-3 catalog [4]. We consider BBHs in GWTC-3 with the False Alarm Rate less than 1 \(\mathrm{yr^{-1}}\). This gives us 69 events [61]. We use the posterior/prior samples that are available at the Gravitational Wave Open Science Center (GWOSC) [62].

In this work, we don’t estimate the total merger rate itself. Therefore, we marginalize the total coalescence rate of BBHs by assuming a log-flat prior [59]. In that case, we can deduce the posterior probability density of branching fraction \(\require{physics} \Lambda = \qty{f_j}\), \(\require{physics} p(\Lambda\mid\qty{x}),\) given observation set \(\require{physics} \qty{x}\): \[\require{physics} \label{eq:final} p(\Lambda\mid\qty{x})\propto\pi(\Lambda)\prod_{i=1}^{N_\mathrm{obs}}\dfrac{1}{\alpha(\Lambda)}\int\dd{\theta_i}\dfrac{p(\theta_i\mid x_i)}{\pi(\theta_i)}p(\theta_i\mid\Lambda),\tag{1}\] The posterior probability density describes the probability density given the observations \(\require{physics} \qty{x}\) to infer the unknown population parameter \(\Lambda\) through the binary parameters \(\theta\). A detailed derivation of this equation is given in the Appendix 6. The \(\alpha(\Lambda)\) in this equation is a quantity that can be called "detection efficiency" and will be explained later. \(p(\theta_i\mid x_i)\) is the posterior probability density for estimating the source parameters based on the observed data \(x_i\), and \(\pi(\theta_i)\) is the prior probability density of the parameter \(\theta_i\). \(p(\theta_i\mid\Lambda)\) is the probability density that describes the probability of an event occurring when the hyperparameter is \(\Lambda\) and the true event parameter is \(\theta\). Since the integral on the multidimensional parameter space is computationally expensive, we approximate the integral by the discrete sum of the posterior samples of each event available in the LVK GWTC-3 at GWOSC. We have \[\require{physics} \label{eq:approxed} p(\Lambda\mid\qty{x})\propto\pi(\Lambda)\prod_{i=1}^{N_\mathrm{obs}}\dfrac{1}{\alpha(\Lambda)}\dfrac{1}{S_i}\sum_{k=1}^{S_i}\dfrac{p(\theta_i^k\mid\Lambda)}{\pi(\theta_i^k)},\tag{2}\] where \(S_i\) denotes the total number of posterior samples used in calculating the sum, and \(\require{physics} \qty{\theta_i^k}\) are the LVK posterior samples.

The last term \(p(\theta_i\mid\Lambda)\) in (1 ) is given as \[p(\theta\mid\Lambda) = \sum_j f_j p(\theta\mid \mu^j),\] where \(p(\theta\mid \mu^j)\) is a probability density distribution that describes the probability that a source with the parameter \(\theta\) is generated in \(j\)-th formation channel, represented by \(\mu^j\).

Under this formula, we can express \(\require{physics} p(\Lambda\mid\qty{x})\) as \[\require{physics} p(\Lambda\mid\qty{x})\propto\prod_{i=1}^{N_\mathrm{obs}}\dfrac{1}{\alpha(\Lambda)}\sum_j\dfrac{f_j}{S_i}\sum_{k=1}^{S_i}\dfrac{p(\theta_i^k\mid\mu^j)}{\pi(\theta_i^k)}.\]

In this work, we use four parameters for \(\theta\): the source frame chirp mass \(\mathcal{M},\) the mass ratio \(q,\) the effective inspiral spin parameter \(\chi_\mathrm{eff}\), and the merger redshift \(z\).

The detection efficiency \(\alpha(\Lambda)\) is computed as \[\require{physics} \begin{align} \alpha(\Lambda) &= \int\dd{\theta} p(\theta\mid\Lambda)p_\mathrm{det}(\theta) \notag\\&= \sum_{j}f_j\int\dd{\theta} p(\theta\mid\mu^j)p_\mathrm{det}(\theta) \eqqcolon \sum_j f_j\alpha_j, \end{align}\] where we set \[\require{physics} \alpha_j \mathrel{\vcenter{:}}= \int\dd{\theta} p(\theta\mid\mu^j)p_{\mathrm{det}}(\theta).\] In this equation, \(p_\mathrm{det}(\theta)\) is the probability that a binary merger with a true event parameter of \(\theta\) can be detected in noisy observations: \[\require{physics} p_\mathrm{det}(\theta) = \int\dd{x} p(x\mid\theta),\] where the integration sums up all detectable GW signals.

In this work, we use the injection set from GWTC-3 [62] to evaluate the detection efficiency by using the importance sampling method. We also tested the evaluation method applied in [43] and found that the evaluated values are consistent with our method.

Recent population studies have some criteria for \(\alpha(\Lambda)\) to reject unphysical hyperparameters (e.g., [63]), but we do not have them because imposing those restrictions would severely limit the possible Hyperparameters. For \(p(\theta\mid\mu^j)\), we construct the kernel density estimator (KDE) from the calculations of the \(j\)-th astrophysical channels. Some of the formation channels have a peak at \(q=1\), which means that several formation channels tend to have equal-mass binary systems. However, since \(q\leq 1\) is a physical constraint of the parameter, the simple KDE method fails to represent the distribution around \(q=1\). Therefore, we applied the reflection method [64] to the data point \(q=1\). Finally, we calculate the prior distribution \(\pi(\theta)\), which was used by LVK for each BBH candidate. For the event prior to GWTC-3, \(\pi(\theta)\) is uniform in redshifted component mass, uniform in spin magnitude, and isotropic in spin orientation. There are two variants of posterior samples for each event, the difference being the distance prior. As in [61], we use the ‘nocosmo’ file to rely on the abundant number of posterior samples [2]–[4]. We utilize the Dynesty package for plotting the posterior distribution.

3 Results↩︎

The result from a synthesis of five formation channels is shown in Fig. 1. This shows the inferred underlying fraction, representing the intrinsic contribution of each formation channel to the astrophysical BBH population in the Universe. In this mixing scheme, the dominant formation channel is Pop. I+II , with an underlying fraction of \(0.77^{+0.11}_{-0.15}\) at the 90% confidence interval. On the other hand, if one takes the PopIII channel into account, the branching fraction of this formation channel accounts for \(0.14^{+0.16}_{-0.11}\) of the total. CHE channel contributes quite a low percentage in this synthetic scheme, with \(<0.6\%\) contribution estimated at the 90% confidence interval. The contribution from dynamical evolutions, evolution in GC or NSC, accounts for \(0.08^{+0.07}_{-0.04}\), and the contribution from the GC is greater than that of the NSC.

From the sensitivity estimation of the model and the underlying fraction parameter, one can define the detectable fraction of a channel [43]: \[\label{eq:32detectable95fraction} f^\mathrm{det}_j = \frac{f_j\alpha(\mu_j)}{\sum_kf_k\alpha(\mu_k)},\tag{3}\] which indicates the contribution to the observed GW catalog from each population channel. Fig. 2 demonstrates the 5-dimensional distribution in terms of detectable fraction. We see that the contribution from Pop. I+II channel is disfavored in the detectable fractions. This is because most of the BBHs that can be made in this channel are less massive, so the selection effect is relatively small. The selection effect in this calculation for this channel is \(\alpha(\mu_\mathrm{I+II}) = 3.4\times10^{-4}\), which is at least 10 times smaller than other formation channels. Instead, the contribution from dynamical formation channels is pushed in terms of detectable fractions. The contribution from the CHE channel to the detectable fraction is zero-consistent, but it rises to \(0.03^{+0.06}_{-0.02}\) at 90% credibility.

To evaluate the presence of the PopIII formation channel, we compared two nested models: A model that incorporates five channels, and a model in which the PopIII formation channel does not contribute to the catalog at all. We evaluate the likelihood ratio among these models and find that the existence of PopIII star formation channel improves the likelihood by \(11.87,\) or \(1.07\) in \(\log_{10}\) scale. This big difference suggests that the PopIII formation channel captures important aspects of GWTC-3 that cannot be adequately represented by other channels.

In this alternative scheme, the domination of the isolated evolution of Pop. I+II is much stronger, with a branching ratio of \(0.89^{+0.05}_{-0.08}\) at the 90% confidence interval, while the contribution from dynamical evolutions is also greater, \(0.07^{+0.07}_{-0.04}\) and \(0.04^{+0.03}_{-0.02}\) for GC and NSC, respectively.

4 Discussion↩︎

We analyze the latest BBH merger catalog to infer the fractions of astrophysical BBH formation channels, including the isolated evolution of PopIII stars. The analysis incorporates five binary formation channels: Pop. I+II , PopIII, CHE, GC, and NSC.

The mixed model results demonstrate that isolated formation channels dominate BBH production, accounting for approximately 90% of astrophysical merger events. This finding is consistent with a previous study [43], which also demonstrated that isolated dominance can break down when natal BH spins are large. In our analysis, we assume natal spins for isolated BBHs are less than 0.2, making the observed Pop. I+II dominance consistent with theoretical expectations.

In our mixed model analysis, the CHE channel contributes only a small fraction to the detected population, despite its characteristic mass scale around \(\sim 30\,M_\odot\) like the PopIII channnel. This is mainly because the CHE channnel predicts a broad positive \(\chi_{\rm eff}\) distribution with a peak around \(\chi_{\rm eff}\sim 0.5\), whereas the current GWTC-3 BBH sample does not favor such high effective spins. Therefore, the relative detectable fractions are governed by the combined effects of the mass distribution, spin distribution, and detector selection effects. In addition, we find that the PopIII channel contributes more significantly than the dynamical formation channels (GC and NSC). However, our results show a stronger preference for the Pop. I+II channel compared to [43]. This discrepancy may arise from differences in the underlying population synthesis models employed in each study. This highlights a fundamental challenge in population inference studies: results can vary substantially depending on the choice of population synthesis framework. The sensitivity to model assumptions underscores the importance of continued refinement of theoretical predictions and the need for systematic comparison across different modeling approaches.

Our likelihood ratio analysis demonstrates that incorporating the PopIII channel increases the model likelihood by a factor of \(\sim\)​10. According to Jeffreys’ scale for Bayes factors, this constitutes strong evidence, though it falls short of being decisive. However, this likelihood improvement must be interpreted cautiously due to the nested nature of our model comparison framework. The observed enhancement can partially reflect the inherent advantage conferred by additional parameter flexibility rather than necessarily indicating superior theoretical validity. The inclusion of an additional formation channel naturally increases the model’s capacity to accommodate the observed data through increased degrees of freedom. Consequently, the likelihood difference alone cannot definitively establish whether the PopIII channel provides a more accurate or physically motivated explanation of the GWTC-3 observations.

Using the inferred relative fraction parameters, one can obtain astrophysical predictive distributions within this framework for physical quantities.

The mass distribution of formation channels that are used in this analysis is depicted in Fig. 3. Pop. I+II channel mainly contributes to the BHs with mass \(\leq20\;M_\odot\). The contribution from CHE models is concentrated at \(m_1\approx30\;M_\odot\). PopIII, GC, and NSC channels contribute to the mass distribution for a wide mass range. Fig. 4 shows the recovered distribution with 90% credibility from the synthesis of the formation channels. The \(m_1\) distribution is highly peaked at \(m_1\approx10\;M_\odot\), from Pop. I+II channel. Distribution of \(m_1\gtrsim20\;M_\odot\) is more uncertain compared to the less massive region, reflecting the uncertainty of the compositions of PopIII, GC, and NSC channels.

Fig. 5 illustrates the \(\chi_\text{eff}\) distribution of formation channels. Though the treatment of the spin in population synthesis of isolated evolution has uncertainty, they tend to have positive \(\chi_\text{eff}\). This positive alignment is expected because the stellar spins of the progenitor stars are generally aligned with the orbital angular momentum in isolated binaries. The \(\chi_\text{eff}\) distribution for the CHE channel has a peak around \(\chi_\text{eff}\approx0.5\), and it is widely distributed with positive \(\chi_\text{eff}\) spanning from \(\chi_\text{eff}\approx0\) to \(\chi_\text{eff}\approx1\). This broad but positive distribution reflects the efficient angular momentum transfer in CHE channel.

The \(\chi_\text{eff}\) distribution of dynamical evolutions can be decomposed as a peak at \(\chi_\text{eff}=0\) and a flat distribution around \(\abs{\chi_\text{eff}}\approx0.5\). The former can be explained by the spin orientation being randomized through many-body interactions in dense stellar environments, leading to no preferred spin alignment. The latter is the contribution from the hierarchical mergers, as the massive BH from a BH merger has a spin magnitude around \(\abs{\chi}\approx 0.6\) due to the angular momentum conservation during the merger process. The NSC channel tends to have more hierarchical mergers compared to the GC channel, primarily because the deeper gravitational potential and higher escape velocity in NSC allow them to retain more merger products.

Fig. 6 shows the \(\chi_\text{eff}\) distribution obtained from the population inference. The uncertainty of spin distribution at \(\chi_\text{eff}>0.2\) can be mainly attributed to the uncertainty of the fraction of the PopIII channel, while no event to date has such a high \(\chi_\text{eff}\) parameter with high credibility. Thus, further subdivision of the isolated evolution channel would allow for more precise estimates. This approach could provide more detailed physical insight.

Compared to the latest population inference from LVK collaboration [65], the predicted distribution for the low-mass end of \(m_1\) exhibits the structure around \(\sim 20M_\odot\) noted in GWTC-4 in a more pronounced form, whilst a dip around \(\sim15 M_\odot\) has newly emerged. Within our framework, these characteristics are largely attributed to the Pop. I+II model. Stars with low metallicity (\(\sim0.1Z_{\odot}\)) tend to produce black holes with masses around \(\sim20M_{\odot}\) [66], which gives rise to a local enhancement in the black-hole mass distribution. For the high-mass end, the predictive distribution from our framework drops earlier than the LVK’s one, but a relatively small number of such events makes the uncertainty of this region bigger. Further observations in this region are awaited. For the \(\chi_\text{eff}\) distribution, the discrepancy between this framework and GWTC-4 is much clearer. The distribution for \(\chi_\text{eff}\) in GWTC-4 analysis is fitted with a (skewed-)Gaussian distribution, while the predictive distribution from this framework shows multiple features because of the contributions of each subpopulation. If one considers the integrated subpopulations to be the true picture, the shape of the effective spin distribution appears highly jagged whch may have been smoothed out by the fit by a Gaussian-like distribution. However, further accumulation of events is awaited to clarify this.

We note that our results are subject to several sources of uncertainty. For PopIII stars, the assumed initial mass function (IMF) can affect the BH mass distribution, particularly the slope of the high-mass tail, although the characteristic mass scale around \(\sim 20\)–\(30M_\odot\) is expected to be relatively insensitive to the IMF [11], [57]. For Pop. I+II binaries, uncertainties in the treatment of the common-envelope phase can significantly impact the formation efficiency and orbital properties of merging systems [43], [67]–[69]. In addition, for dynamically formed binaries such as those in globular clusters, the assumed initial spins of BHs can affect the predicted spin properties of the population [43]. These uncertainties may quantitatively affect our results, and a more systematic exploration is left for future work. In this study, we examine five formation pathways, but note that these selections are not determinative. Various other intriguing formation pathways exist in reality. One notable example of a formation channel not included in our analysis is the active galactic nucleus (AGN) channel. This channel considers the possibility that binary black holes (BBHs) form and merge within the dense gas disks surrounding supermassive black holes. We did not include the AGN channel in our mixture model since publicly available population synthesis data for this channel are not currently accessible. Nevertheless, several studies have suggested that the AGN channel could make a non-negligible contribution to the observed BBH population, particularly for systems with high masses or moderate eccentricities [70], [71]. Due to the migration and interaction processes in AGN disks, BBHs formed in such environments can exhibit both aligned and misaligned spin configurations. This would lead to an intermediate distribution of the effective spin parameter, \(\chi_{\mathrm{eff}}\), which differs from that predicted by isolated binary evolution or dense stellar dynamics. Recent theoretical work by Tagawa et al. [37] demonstrates that repeated mergers of stellar-mass black holes can occur efficiently in AGN disks, leading to the formation of massive black holes through hierarchical mergers. Such processes could naturally produce black holes in the pair-instability mass gap and may also result in observable signatures in spin and mass distributions. Furthermore, recent hierarchical Bayesian analyses indicate that the AGN channel might explain some of the features seen in the GWTC-3 catalog, including possible hints of hierarchical mergers or spin misalignment [72]. Although our current framework cannot directly constrain the AGN contribution, the inclusion of this channel in future studies—once detailed population models become available—will be important for a more complete understanding of BBH formation channels. Another dynamical formation channel not included in our current mixture model is the young dense cluster (YDC) channel. YDCs are stellar systems with high densities and a rich population of massive stars during their early evolutionary phases. Frequent dynamical interactions in such environments—particularly binary-single and binary-binary encounters—facilitate the efficient formation and hardening of binary black holes (BBHs), some of which may merge within the cluster’s lifetime [73]. A distinctive feature of YDC-origin BBHs lies in their spin distribution. Black holes that form directly from massive stellar collapse in YDCs are generally expected to have low natal spin magnitudes, assuming efficient angular momentum transport in progenitor stars. However, YDCs are also favorable environments for hierarchical mergers: a black hole formed in an earlier merger may participate in subsequent mergers due to the high interaction rates. Such hierarchical mergers typically produce remnant black holes with moderate to high spin (\(a \sim 0.6\)–\(0.8\)), depending on the mass ratio and spin alignment of the progenitor BBHs [74]. The mass distribution of merging BBHs in young dense cluster models appears to be an approximately log-flat distribution, i.e., \(p(M) \propto M^{-1}\) from 5 \(M_{\odot}\) to 50 \(M_{\odot}\) [73]. Therefore, YDCs can host a population of massive black holes with higher spin than what is predicted by isolated binary evolution, while still allowing a broad distribution in the effective inspiral spin parameter \(\chi_{\mathrm{eff}}\) due to random spin orientations acquired through dynamical processes.

Although we do not include the YDC channel in our current analysis due to the lack of publicly available population synthesis datasets specifically targeting this environment, the potential contribution of YDC-origin BBHs—especially in the form of massive, high-spin systems—makes this channel a strong candidate for inclusion in future work.

Finally, since Pop III stars are the first stars in the universe, they formed much earlier than Pop I/II stars. Previous studies show that the peak of the Pop III BBH formation is located at \(z \sim 10\) [19], [75]. Therefore, next-generation GW detectors will reveal the nature of the BBH population (such as the redshift-spin correlation) more clearly. The prediction of binary parameter distribution under the mixed model will be informative for the upcoming next-generation GW detectors era.

5 Summary↩︎

There are many theoretical pathways proposed to explain the formation of BBHs. In this paper, we take five formation channels: Pop. I+II , PopIII, CHE, GC, and NSC channels. From our mixing scheme revealed that Pop. I+II channel dominates the origin of BBHs, but it also suggested that PopIII channel contributes approximately 10% to BBH mergers across the entire universe. From the perspective of observed catalogues, contributions from diverse formation channels, including channels beyond those mentioned above, were suggested. These findings underscore the importance of simultaneously considering the contributions of various formation processes. It should be noted, however, that the inferred branching fractions \(\{f_j\}\) and the resulting posterior predictive distributions of the source parameters, such as the mass distribution, are inherently conditioned on the choice of formation channel models considered in this work. Each model encodes a set of astrophysical assumptions and idealizations, and the results may shift quantitatively as these prescriptions are refined or as additional formation channels are incorporated. Future improvements in both theoretical models and observational constraints will be essential to refine our understanding of BBH formation pathways. Further observations and the advent of next-generation detectors will lead to the accumulation of events with increasingly precise analysis results in the future. These events will enable the imposition of stricter constraints on the contribution rates from each channel. Furthermore, incorporating additional formation channels and reducing model uncertainties through multi-messenger observations and refined stellar evolution physics will enable us to construct a more complete picture of the astrophysical BBH formation. Ultimately, disentangling the origins of BBHs will not only illuminate the diverse evolutionary paths of massive stars and compact objects but also provide crucial insights into the star formation history and chemical evolution of the universe across cosmic time.

Acknowledgement↩︎

M. I. was supported by Forefront Physics and Mathematics Program to Drive Transformation (FoPM), a World-leading Innovative Graduate Study (WINGS) Program, the University of Tokyo. This research has made use of data obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO Laboratory and Advanced LIGO are funded by the United States National Science Foundation (NSF) as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; National Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea; Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan. This work was supported by JSPS Grant-inAid for Scientific Research on Innovative Areas 2905: JP17H06358, JP17H06361 and JP17H06364, JSPS Core-to-Core Program A. Advanced Research Networks, JSPS Grant-in-Aid for Transformative Research Areas (A) 20A203: JP20H05854, the joint research program of the Institute for Cosmic Ray Research, University of Tokyo.

T. K. acknowledges support from JSPS KAKENHI Grant Numbers JP21K13915 and JP22K03630 and the financial support from the Science Moves Award.

6 Hierarchical Bayesian Analysis↩︎

Here, we describe a summary of the hierarchical Bayesian analysis method. In a Bayesian manner, given a set of data \(\require{physics} \qty{x}\) from observation, one can calculate the posterior distribution for population parameter \(\Lambda\), \(\require{physics} p(\Lambda\mid\qty{x})\), as \[\require{physics} p(\Lambda\mid\qty{x}) \propto \pi(\Lambda)p(\qty{x}\mid\Lambda),\] where \(\pi(\Lambda)\) is the prior distribution for the population parameter \(\Lambda\). If we assume that all events are independent of each other, this probability \(\require{physics} p(\Lambda\mid\qty{x})\) can be expressed as \[\require{physics} p(\Lambda\mid\qty{x}) = \pi(\Lambda)\prod_{i=1}^{N_\mathrm{obs}}p(x_i\mid\Lambda),\] where \(x_i\) denotes the observation data of \(i\)-th event. For each GW event, we estimate the source parameter, such as masses, and we denote this by \(\theta.\) The probability of detecting the GW event \(x\) from population parameter \(\Lambda\) can be described by marginalizing over \(\theta\): \[\require{physics} \label{eq:naive95p} p(x\mid\Lambda)\propto \int \dd{\theta}p(x\mid\theta)p(\theta\mid\Lambda).\tag{4}\] Since we have the detection, \(x\) must be ‘detectable.’ In other words, \(x\) must pass the detection criteria. Therefore, the normalization for equation (4 ) must be \[\require{physics} \int \dd{x}p(x\mid\Lambda) = 1,\] where the integral domain is the entire \(x\) satisfying detection criteria. This requires a parameter so-called ‘detection efficiency’ \(\alpha(\Lambda)\) as a normalization factor: \[\require{physics} \alpha(\Lambda) = \int\dd{\theta}p_\mathrm{det}(\theta)p(\theta\mid\Lambda),\] where \(p_\mathrm{det}(\theta)\) is the probability that a gravitational wave from a binary system with binary parameter \(\theta\) is observed by the detectors. With the normalization factor, the equation to assess \(p(x\mid\Lambda)\) becomes \[\require{physics} p(x\mid\Lambda)= \dfrac{1}{\alpha(\Lambda)}\int \dd{\theta}p(x\mid\theta)p(\theta\mid\Lambda).\]

For every event in the GWTC-3, LVK collaboration estimates the source parameter on it, applying Bayesian analysis [4]. Hence, \(p(x\mid\theta)\) can be converted to \(p(\theta\mid x)p(x)/\pi(\theta)\) by applying Bayes’ theorem: \[\require{physics} p(x\mid\Lambda)= \dfrac{1}{\alpha(\Lambda)}\int \dd{\theta}\dfrac{p(x)p(\theta\mid x)}{\pi(\theta)}p(\theta\mid\Lambda).\]

Finally, we have \[\require{physics} p(\Lambda\mid\qty{x})= \pi(\Lambda)\prod_{i=1}^{N_\mathrm{det}}\qty[\dfrac{1}{\alpha(\Lambda)}\int \dd{\theta}\dfrac{p(\theta\mid x_i)}{\pi(\theta)}p(\theta\mid\Lambda)].\]

References↩︎