What You Don't Know Can Hurt You: Use and Abuse of Astrophysical Models in Gravitational-wave Population Analyses

We use hierarchical Bayesian inference on the branching fractions between different astrophysical formation channels of BBHs. Our methods mostly follow the analysis developed by Zevin et al. (2021a), adapting their codebase, Astrophysical Model Analysis and Evidence Evaluation (AMA ${ \mathcal Z }$ E), for our work. Here, we outline the essentials of this method, as well as the key differences.

We consider five different formation channels: three isolated evolution (field) channels and two dynamical formation channels. The common envelope (CE) (e.g., Paczynski 1976; van den Heuvel 1976; Tutukov & Yungelson 1993; Bethe & Brown 1998; Belczynski et al. 2002, 2016; Dominik et al. 2012; Eldridge & Stanway 2016; Stevenson et al. 2017; Giacobbo & Mapelli 2018) and stable mass transfer (SMT) (van den Heuvel et al. 2017; Neijssel et al. 2019; Gallegos-Garcia et al. 2021; van Son et al. 2022) scenarios are field channels that involve unstable and stable mass transfer, respectively, following the formation of the first black hole. In the chemically homogeneous evolution (CHE) channel, stars in a close, tidally locked binary rotate rapidly, causing temperature gradients that lead to efficient mixing of the stars' interiors. The stars do not undergo significant expansion of the envelope, preventing significant post-main-sequence wind mass loss and premature merging, resulting in higher-mass BBHs (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016). Finally, the two dynamical formation channels lead to merging black holes via strong gravitational encounters that harden the binary in cluster cores (e.g., McMillan et al. 1991; Hut et al. 1992; Sigurdsson & Phinney 1993; Portegies Zwart & McMillan 2000; Miller & Hamilton 2002; Gültekin et al. 2006; O'Leary et al. 2006; Fregeau & Rasio 2007; Downing et al. 2010; Samsing et al. 2014; Ziosi et al. 2014; Rodriguez et al. 2015, 2016; Antonini & Rasio 2016; Askar et al. 2017; Samsing & Ramirez-Ruiz 2017; Kremer et al. 2020), where heavy black holes migrate toward due to dynamical friction (Lightman & Shapiro 1978; Sigurdsson & Hernquist 1993); we consider dynamical formation of BBHs in globular clusters (GCs) and nuclear star clusters (NSCs). See Zevin et al. (2021a) for a more detailed description of the astrophysical models considered in this work.

Each of these channels is modeled to predict the four-dimensional distribution of the BBHs it forms with parameters ${\boldsymbol{\theta }}=[{{ \mathcal M }}_{{\rm{c}}},q,{\chi }_{\mathrm{eff}},z]$, where ${{ \mathcal M }}_{{\rm{c}}}$ is the source-frame chirp mass, q is the mass ratio (defined to be 0 < q ≤ 1 ), χ_eff is the effective dimensionless spin parameter, and z is the redshift; these are constructed into a probability distribution using a four-dimensional kernel density estimate (KDE) bounded by the physical constraints of each parameter. The models also depend on two additional parameters encoding uncertainties in the physical prescription: χ_b, the dimensionless spin of a black hole formed in quasi-isolation directly following core collapse, and α_CE, which parameterizes the efficiency of common envelope ejection (see, e.g., Ivanova et al. 2013). ⁷ We note that the choice of natal black hole spin χ_b does not set all black holes in a given population to merge with this exact spin, because tidal spin-up processes (Qin et al. 2018; Zaldarriaga et al. 2018; Bavera et al. 2021) and hierarchical mergers (Pretorius 2005; González et al. 2007; Buonanno et al. 2008) can increase the spin of black holes that participate in BBH mergers. We assume these two parameters take on a grid of discrete values (χ_b ∈ [0.0, 0.1, 0.2, 0.5], α_CE ∈ [0.2, 0.5, 1.0, 2.0, 5.0]) over which we compute the models, although α_CE only affects the CE channel in our models. A plot of the detection-weighted marginalized model KDEs for α_CE = 1.0, χ_b = 0.2 is shown in Figure 1 as an example. Further details, including formation models, detection weighting, and mathematical framework, can be found in Zevin et al. (2021a).

Zoom In Zoom Out Reset image size

Overall, we perform hierarchical inference on seven hyperparameters ⁸ Λ = [ β , χ_b, α_CE], where β = [β_CE, β_CHE, β_SMT, β_GC, β_NSC] are the five astrophysical formation channel branching fractions. The steps involved in performing hierarchical inference on GW populations given a set of posterior samples (${\boldsymbol{\theta }}=[{{ \mathcal M }}_{{\rm{c}}},q,{\chi }_{\mathrm{eff}},z]$ in our case) for N_obs sources in the presence of selection effects have been thoroughly discussed in the literature (Farr et al. 2015; Mandel et al. 2019; Thrane & Talbot 2019; Vitale et al. 2022b; Abbott et al. 2023), and we therefore do not review them here. As in Zevin et al. (2021a), we use a broad, agnostic prior for the hyperparameters: a flat symmetric Dirichlet prior for β and a uniform prior for χ_b and α_CE over the allowed discrete values. The outputs of the inference are samples of the hyperposterior p( β ∣χ_b, α_CE) over the grid of hypermodels (i.e., allowed values of χ_b and α_CE), from which we can compute the marginalized hyperposterior p( β ) as well as the Bayes factors ${\mathrm{BF}}_{b}^{a}$ of hypermodel a compared to hypermodel b. Bayes factors >1 indicate that model a is more supported than model b, with ${\mathrm{BF}}_{b}^{a}\gt 10\,(100)$ indicating strong (decisive) support for model a over model b (Kass & Raftery 1995).

We also recover the detectable branching fractions ${\beta }^{\det }$, which represent the fraction of detectable BBHs originating from each channel. These are computed by rescaling each underlying branching fraction by its detection efficiency, defined as ${\xi }_{j}^{{\chi }_{{\rm{b}}},{\alpha }_{\mathrm{CE}}}=\int {P}_{\det }({\boldsymbol{\theta }})p({\boldsymbol{\theta }}| {\mu }_{j}^{{\chi }_{{\rm{b}}},{\alpha }_{\mathrm{CE}}})\,d{\boldsymbol{\theta }}$, where ${P}_{\det }({\boldsymbol{\theta }})$ is the probability of detecting a BBH with parameters θ and $p({\boldsymbol{\theta }}| {\mu }_{j}^{\chi ,\alpha })$ is the probability of formation channel j producing a BBH with parameters θ , dependent on the model μ_j and choice of physical prescription χ_b and α_CE (Zevin et al. 2021a). We apply this method to both real (Section 2.2) and simulated (Sections 3 and 4) BBH data.

We extend the work of Zevin et al. (2021a) by applying the inference to confident BBH detections up to the second half of the third observing run (O3b). We use the publicly released priors and posterior samples from the GWTC-2.1 (The LIGO Scientific Collaboration et al. 2021a) and GWTC-3 (The LIGO Scientific Collaboration et al. 2021b) analyses; detection probabilities are calculated for the LIGO-Hanford, LIGO-Livingston, and VIRGO network operating at midhighlatelow sensitivities (Abbott et al. 2018). There are two key differences between our analysis and that of Zevin et al. (2021a). First, we apply a more stringent detection threshold of FAR ≤ 1 yr⁻¹ (Abbott et al. 2023); therefore, events GW190424_180648, GW190514_065416, and GW190909_114149, which were used in the previous analysis, are now excluded. Second, we evaluate the prior at the posterior points analytically rather than by constructing a KDE from prior samples; for further discussion on this point, see Appendix A. Finally, as in Zevin et al. (2021a), we opt to exclude GW190521 and GW190814 from the analysis, as their posteriors extend significantly to regions of our model KDEs with little to no support. Furthermore, in the case of GW190814, it is uncertain whether this system is comprised of two black holes or a neutron star and a black hole. We find that the inclusion of GW190521 does not significantly affect our results; see Franciolini et al. (2022) for an analysis that includes GW190521. Overall, we do hierarchical inference on 68 BBHs, compared to the 45 BBHs considered in Zevin et al. (2021a). Unless otherwise specified, we report results as median and 90% symmetric credible intervals.

Figure 2 shows the posterior distributions on the underlying branching fraction β , including detections from O3b; the same plot but for the detectable branching fractions can be found in Appendix B. We find ${\beta }_{\mathrm{CE}}={90}_{-11}^{+5.0} \%$, ${\beta }_{\mathrm{CHE}}={0.6}_{-0.5}^{+1.0} \%$, ${\beta }_{\mathrm{GC}}={3.7}_{-3.4}^{+7.6} \%$, ${\beta }_{\mathrm{NSC}}={1.1}_{-0.8}^{+1.7} \%$, and ${\beta }_{\mathrm{SMT}}={4.0}_{-3.2}^{+6.5} \%$, indicating strong support for the CE channel dominating the underlying astrophysical population in our set of models. However, there is comparable contribution from all five channels to the detectable BBH population: ${\beta }_{\mathrm{CE}}^{\det }={21}_{-10}^{+16} \%$, ${\beta }_{\mathrm{CHE}}^{\det }={6.9}_{-5.2}^{+8.5} \%$, ${\beta }_{\mathrm{GC}}^{\det }={24}_{-21}^{+26} \%$, ${\beta }_{\mathrm{NSC}}^{\det }={20}_{-13}^{+17} \%$, and ${\beta }_{\mathrm{SMT}}^{\det }={24}_{-19}^{+24} \%$. We attribute this difference to the fact that, compared to other channels, the CE channel produces less massive black holes distributed at higher redshifts, which therefore are harder to detect; see Figure 1. With 90% (99%) credibility, no single formation channel contributes to more than 49% (61%) of the detectable BBH population. Additionally, over 98% of posterior samples have significant (${\beta }^{\det }\gt 10$%) contributions from three or more different formation channels. Overall, we find that the CE channel contributes the most to the underlying BBH population; however, as we discuss in Section 4, the posterior for β_CE is also the most uncertain and prior dominated. We additionally find that a mix of formation channels contributes to the detectable BBH population. These results, as well as the overall shapes of the branching fraction posteriors, are consistent with the analysis with GWTC-2 data (Zevin et al. 2021a), which found ${\beta }_{\mathrm{CE}}={71}_{-60}^{+19} \%$ and that no single channel contributed to more than 70% of the detectable BBH population with 99% confidence. No strong correlations are apparent upon examining a corner plot of the branching fractions.

Zoom In Zoom Out Reset image size

Turning to the selection of physical prescription hyperparameters χ_b and α_CE, we find no posterior support for models with χ_b > 0.1; there is no significant preference between the χ_b = 0.0 and χ_b = 0.1 models, with ${\mathrm{BF}}_{{\chi }_{{\rm{b}}}=0.0}^{{\chi }_{{\rm{b}}}=0.1}=1.18$. We favor high common envelope efficiencies, with ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=1.0}^{{\alpha }_{\mathrm{CE}}=5.0}=249$, and strongly disfavor low common envelope efficiencies, with ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=1.0}^{{\alpha }_{\mathrm{CE}}=0.2}=5\times {10}^{-5}$. These results are also consistent with Zevin et al. (2021a). The main difference is that we obtain stronger constraints, both in singling out preferred hypermodels and in the uncertainties (widths) of the branching fraction hyperposteriors. We attribute this effect to the increase in sample size from GWTC-2 to GWTC-3, as well as the difference in method in evaluating the prior at the posterior points during the inference; see Appendix A for further discussion.

Next, we perform the same inference using simulated BBH observations in a universe where the BBH population exactly follows our models. The motivation for this analysis is to to quantify how uncertainties in the hyperposterior scale with the number of observed events.

We first create a mock "universe," i.e., a set of true values for our hyperparameters Λ = [ β , χ_b, α_CE]. From each formation channel j, we draw n_j = β_j n BBHs from the population model $p(\theta | {\mu }_{j}^{{\chi }_{{\rm{b}}},{\alpha }_{\mathrm{CE}}})$, for a total of n = 5 × 10⁴ BBHs that form our underlying population. Next, we draw from this underlying population, assigning extrinsic parameters (sky location and inclination) from an isotropic distribution. For each BBH system, we calculate its optimal signal-to-noise ratio ρ_opt assuming a network consisting of LIGO-Hanford, LIGO-Livingston, and Virgo operating at O4 low (LIGO) and high (Virgo) sensitivities (Abbott et al. 2018; O'Reilly et al. 2022), and keep only the mock signals with ρ_opt ≥ 11; we repeat this process until we have a mock catalog of N_obs detections. Then, we perform parameter estimation on the N_obs BBHs. For both the S/N calculation and parameter estimation, we use the Bayesian inference software bilby (Ashton et al. 2019) and the IMRPhenomXP waveform approximant (Pratten et al. 2021). Finally, we use these posterior samples for the hierarchical inference analysis outlined in Section 2.1. For consistency, we use the same O4 sensitivities for the detection weighting during the inference as the above computation of mock signal S/Ns.

The various mock universes that we use in this paper, along with the sections in which they are discussed, are summarized in Table 1. For this analysis, we choose a fiducial unequal mixture of formation channels with underlying branching fractions β_CE = 40%, β_CHE = 5%, β_GC = 25%, β_NSC = 5%, and β_SMT = 25%, quasi-isolated natal spins of χ_b = 0.0, and a CE efficiency of α_CE = 1.0.

To investigate the scaling of the hyperposterior with N_obs, for each universe we repeat the inference with N_obs = 50, N_obs = 150, and N_obs = 250. N_obs = 250 represents an estimate on the total number of BBH detections anticipated by the end of O4 (Weizmann Kiendrebeogo et al. 2023).

In Figure 3, we plot for our chosen fiducial universe the overall posterior distribution on the underlying branching fractions for different values of N_obs (left column), as well as contributions (to scale) from the most favored values of χ_b and α_CE (columns 2 to 5). First, we do recover the true values of β and χ_b. For N_obs = 250, we find ${\beta }_{\mathrm{CE}}={28}_{-21}^{+26} \%$, ${\beta }_{\mathrm{CHE}}={3.0}_{-1.5}^{+2.2} \%$, ${\beta }_{\mathrm{GC}}={40}_{-16}^{+18} \%$, ${\beta }_{\mathrm{NSC}}={3.6}_{-2.1}^{+3.1} \%$, and ${\beta }_{\mathrm{SMT}}={24}_{-13}^{+16} \% ;$ the true value of the branching fraction falls within the 90% symmetric credible interval of the posterior for all five channels. Furthermore, χ_b = 0, the chosen true value for χ_b, is favored over the next best model, χ_b = 0.1, by a Bayes factor of 3.5 × 10⁵.

Zoom In Zoom Out Reset image size

Unlike χ_b, the model with α_CE equal to the true value is not favored, with a marginal preference for higher values α_CE = 5.0 and α_CE = 2.0 over the true value α_CE = 1.0 with Bayes factors ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=1.0}^{{\alpha }_{\mathrm{CE}}=5.0}=8.9$ and ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=1.0}^{{\alpha }_{\mathrm{CE}}=2.0}=2.5$. Because α_CE only affects the CE channel, it is not too surprising that the inference did not decisively favor any one value, in a universe where most detected BBHs do not come from the CE channel. We note two possible contributing factors to favoring higher values of α_CE over the true value. If BBHs formed in the CE channel disperse their envelopes more efficiently (i.e., have higher values of α_CE ), the resulting BBHs are

• (a)  
  Less massive, leading to lower detection efficiencies and therefore larger measurement uncertainties for this channel. This is because models with larger α_CE have less low-mass binaries merging within the CE itself, leading to more low-mass BBHs being able to form and merge (Bavera et al. 2021).
• (b)  
  Lower spinning, due to the post-CE separations being wider for higher α_CE and therefore less susceptible to tidal spin-up (Qin et al. 2018; Zaldarriaga et al. 2018; Bavera et al. 2021); this leads to more detections with χ_eff closer to χ_b = 0.0, the chosen fiducial value for this universe. In χ_eff space, this is a feature degenerate with dynamical channels such as the GC channel, which also produces BBHs with χ_eff ≈ 0 due to the BBHs produced having isotropic spin orientations.

This result suggests that our result in Section 2.2, that higher common envelope efficiencies are favored, should be viewed with some caution.

More notably, we see convergence of the hyperposteriors toward their true values as we increase the number of detections. First, we highlight the narrowing of the branching fraction posteriors from the first row to the third row of Figure 3. We quantify uncertainties in the branching fraction posteriors by the widths of the 90% symmetric credible intervals, and will hereafter use the two phrases interchangeably. These uncertainties decrease by up to 69% as we go from 50 to 250 mock events (see Table 2). The Bayes factor in favor of the correct χ_b increases by nearly five orders of magnitude, strongly selecting the true value χ_b = 0.0. This is illustrated by the empty plot corresponding to the contribution from χ_b = 0.1 for N_obs = 250 (third row, middle column), indicating negligible support for competing values of χ_b. On the other hand, the inference has increasing support for α_CE = 5.0, but even at 250 mock detections, we only weakly prefer it over the true model α_CE = 1.0.

Table 2. Uncertainties in the Measurement of the Branching Fractions with N_obs Shown in Figures 3 and 4

Nobs	ΔβCE	ΔβCHE	ΔβSMT	ΔβGC	ΔβNSC	${\rm{\Delta }}{\beta }_{\mathrm{CE}}^{\det }$	${\rm{\Delta }}{\beta }_{\mathrm{CHE}}^{\det }$	${\rm{\Delta }}{\beta }_{\mathrm{SMT}}^{\det }$	${\rm{\Delta }}{\beta }_{\mathrm{GC}}^{\det }$	${\rm{\Delta }}{\beta }_{\mathrm{NSC}}^{\det }$	${\mathrm{BF}}_{{\chi }_{{\rm{b}}}=0.1}^{{\chi }_{{\rm{b}}}=0.0}$	${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=5.0}^{{\alpha }_{\mathrm{CE}}=1.0}$
50	49%	11%	52%	17%	56%	17%	19%	55%	34%	49%	8.0	0.81
150	46%	6.2%	48%	9.7%	43%	6.4%	11%	47%	24%	42%	27	0.36
250	46%	3.7%	33%	5.2%	28%	4.6%	7.2%	30%	14%	25%	3.5 × 105	0.11
% change	5.9%	66%	36%	69%	49%	73%	62%	46%	60%	49%

Notes. The rightmost two columns give the Bayes factors of the true value of χ_b and α_CE over the most favored competing value for 50, 150, and 250 mock detections (first, second, and third rows). The rest of the table shows the underlying and detectable branching fraction uncertainties, as well as their percent changes between the first and third rows. We note that these are credible 90% interval widths in percentage points, not percent uncertainties.

Download table as:  ASCIITypeset image

It is worth noting that the scaling of the uncertainty with N_obs varies significantly between formation channels. Some channels have more distinctive features in parameter space (see Figure 1) that make them easier to distinguish during the inference. Additionally, differences in the detection efficiencies of different formation channels likely play a role as well. Figure 5 shows the percent decrease in the uncertainty of branching fraction posteriors from N_obs = 50 to N_obs = 250 against the detection efficiency of its channel ${\xi }^{{\chi }_{{\rm{b}}},{\alpha }_{\mathrm{CE}}}$ for our fiducial values χ_b = 0.0 and α_CE = 1.0. Channels with lower detection efficiencies, most notably the CE channel, appear to scale much more poorly with N_obs than channels with higher detection efficiencies. Indeed, the CE channel tends to produce lower-mass black holes at higher redshifts (due to typically shorter delay), leading to fewer of them being detectable (see Figure 1). We emphasize the effect of the low detection efficiency of the CE channel: despite the fact that 40% of the mock underlying population originates from this channel, only 2, 5, and 17 CE BBHs end up in the mock observations, out of the 50, 150, and 250 total detections, respectively. Not only does the β_CE posterior have the largest uncertainty of all formation channels, it also barely decreases in uncertainty as we increase N_obs. If we instead examine how the detectable branching fractions scale with N_obs in Figure 4, we can see that all detectable branching fractions narrow at similar rates as we increase N_obs from 50 (first row) to 250 (third row). Indeed, we can see in Figure 5 that the detectable branching fractions (dotted line) do not exhibit the dependence of the convergence rate on detection efficiency that the underlying branching fractions (solid line) have.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Already, we can see some of the difficulties that arise from hierarchical Bayesian inference in the face of large measurement uncertainties, selection effects, and degenerate features in parameter space. This is a common theme, and we will explore these problems in further detail in Section 4.

Finally, we repeat this analysis for a different set of hyperparameters consisting of equal branching fractions between all formation channels (β_j = 0.2 for all j), χ_b = 0.2, and α_CE = 1.0 (see the second row of Table 1). We find similar results in the convergence with N_obs. Uncertainties in the underlying branching fraction decrease by up to 47% from N_obs = 50 to N_obs = 250, and consistently with our previous results, the uncertainty in β_CE does not decrease. We find also that the shrinking of the detectable branching fraction (${\beta }^{\det }$) posterior uncertainty is more consistent across different channels than the underlying branching fraction. Increasing N_obs causes a strong preference for the true value of χ_b = 0.2; while we have ${\mathrm{BF}}_{{\chi }_{{\rm{b}}}=0.1}^{{\chi }_{{\rm{b}}}=0.2}=1.0$ for N_obs = 50, there is no posterior support for χ_b ≠ 0.2 at N_obs = 150 and N_obs = 250. Similarly to the previous example, there is no strong preference for the true value of α_CE, although it is slightly favored with ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=0.5}^{{\alpha }_{\mathrm{CE}}=1.0}=1.7$. Figures showing the marginalized posteriors on β and ${\beta }^{\det }$ for this set of hyperparameters can be found in Appendix B.

For each channel j, we perform hierarchical inference on N_obs = 100 mock detections in a universe where the entire underlying BBH population originates from channel j (i.e., β_j = 1). We choose for our true values of the physical prescription χ_b = 0.0 and α_CE = 1.0, although the latter choice only affects the CE-dominated universe. Figure 6 shows the branching fraction posteriors and support for different values of χ_b for each single-channel dominated universe.

Zoom In Zoom Out Reset image size

In general, none of the branching fraction posteriors have significant support for β = 1, even for the channel that is actually producing the entirety of the BBHs. This happens because we use a flat, symmetric Dirichlet distribution for our prior, which results in a prior preference for a mixture of channels rather than a single dominating channel; this is illustrated by the gray curves in Figure 6.

The fifth percentiles of the dominating-channel branching fraction posteriors for the CE, CHE, GC, NSC, and SMT-dominated universes are β_CE = 94%, β_CHE = 62%, β_GC = 41%, β_NSC = 31%, and β_SMT = 51%, respectively. The degree to which we underestimate the contribution from the dominating channel varies significantly depending on the channel. We again highlight the effect of detection efficiency. As we saw in Figure 5, the CE, GC, and SMT channels have the lowest detection efficiencies, especially the CE channel. Across the different universes under consideration, β_CE, β_GC, and β_SMT (first, middle, and last columns, respectively) have larger uncertainties that compete with and take away from the branching fraction posterior of the dominating channel; because of the lower detection efficiencies of these channels, it is difficult to discern whether this channel is nonexistent or if its contribution to the full set of detected observations is minor. The case of β_CE is particularly severe: as it does not strongly deviate from the prior, the 95th percentile for β_CE is greater than 32% in all universes for which there is no CE contribution to the underlying population. The CE-dominated universe (first row) does not suffer from this effect; as a result, only in that universe do we recover with narrow precision that the dominating channel is indeed dominating. In the opposite case, the NSC channel has the highest detection efficiency. In the NSC-dominated universe (fourth row), the median value of the β_NSC posterior is less than 0.5, with large contributions to the BBH population from the channels with lower detection efficiencies (β_CE, β_GC, and β_SMT).

There are also several interesting features of the spin model selection. First, in the CHE-dominated universe (second row), the inference is not able to select the true natal black hole spin of χ_b = 0.0, and instead gives approximately equal support to χ_b = 0.0, 0.1 and 0.2, as illustrated by the similar heights of the different colored curves. Let us recall that, while the aligned-spin CE and SMT field channels mostly produce BBHs with χ_eff ≈ χ_b (although CE has a tail for higher spins from tidal spin-up) and the isotropic-spin GC and NSC dynamical channels produce BBHs scattered around χ_eff ≈ 0, the CHE channel uniquely produces BBHs with χ_eff ≳ 0.2 irrespective of χ_b, due to strong tidal spin-up effects (see Figure 1). As a result, the inference is not able to discern between values of χ_b between 0 and the true value 0.2. In the universes dominated by dynamical formation channels (GC, third row, and NSC, fourth row), the selection of the true value of χ_b = 0.0 is less strong, with non-negligible support for χ_b = 0.1 as shown by the green curves. In these universes, most detections have χ_eff ≈ 0, and χ_b only affects the width of this distribution, which is a weaker feature to detect. Additionally, the subpopulation of hierarchical mergers in these channels help to drive the mild support for χ_b > 0; at lower values of χ_b, hierarchical mergers occur more readily, due to weaker gravitational recoil kicks. Hierarchical mergers can have ∣χ_eff∣ significantly greater than zero (Pretorius 2005; González et al. 2007; Buonanno et al. 2008; Kimball et al. 2020), as illustrated by the wings of the marginalized χ_eff distribution for the GC and NSC channels in Figure 1. Because hierarchical mergers form a larger subpopulation in NSC, due to their deeper potential wells and ability to retain post-merger black holes (Miller & Lauburg 2009; O'Leary et al. 2009; Hong & Lee 2015; Antonini & Rasio 2016), this effect is greater for the NSC-dominated universe. Additionally, in the NSC-dominated universe, the β_NSC posterior for χ_b = 0.1 (green curve) peaks at a higher value than χ_b = 0.0 (blue curve). This is because, when χ_b = 0.1, field channels produce fewer of the χ_eff ≈ 0 BBHs that make up the bulk of the population.

Finally, although it is not shown in Figure 6, we comment on the inference on α_CE. First, we do favor the true value of α_CE = 1.0 for the CE-dominated channel, preferring it over the next most favored model with ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}=2.0}^{{\alpha }_{\mathrm{CE}}=1.0}=72$, as expected from a universe where all detections are from the CE channel. We also note consistency with the results of Section 3.2 in the recovery of α_CE for universes where there is no contribution from the CE channel; there is a small bias toward higher values of α_CE. We find ${\mathrm{BF}}_{{\alpha }_{\mathrm{CE}}\lt 1.0}^{{\alpha }_{\mathrm{CE}}\gt 1.0}=1.52,1.54,1.68,$ and 1.58 for the CHE, GC, NSC, and SMT-dominated universes, respectively. ⁹

To examine the effect of χ_b on the hyperposterior, we perform hierarchical inference on different universes with the same underlying branching fractions but different true values of χ_b, using 250 mock detections. To isolate the effect of χ_b, we choose an equal mixture of formation channels in the underlying population (β_j = 0.2 for all channels j). We choose α_CE = 1.0.

Figure 7 shows the marginalized branching fraction posteriors for universes with χ_b = 0.0 and χ_b = 0.2. Consistently with our results in Section 3.2, we recover the true values of the branching fractions as well as the true value of χ_b for both universes. Here, we can see that the selection of χ_b is nonlinear: it is harder to distinguish between lower black hole spins (i.e., χ_b = 0.0 versus χ_b = 0.1 ) and higher ones. While there is no support for other values of χ_b in the posterior of the χ_b = 0.2 universe, there is still non-negligible support for χ_b = 0.1 in the χ_b = 0.0 universe, with ${\mathrm{BF}}_{{\chi }_{{\rm{b}}}=0.1}^{{\chi }_{{\rm{b}}}=0.0}=30$. This is as expected; it is difficult to discern between slowly spinning and nonspinning populations, due to the inherent measurement uncertainty of GW observations (Baird et al. 2013; Vitale et al. 2014; Pürrer et al. 2016; Vitale et al. 2017a).

Zoom In Zoom Out Reset image size

Figure 8 shows a corner plot of the branching fraction posteriors for both universes. Here, we can see the same effect: despite having the same number of mock detections, the data are more informative (yielding a posterior more different from the prior) in the universe with nonzero χ_b.

Zoom In Zoom Out Reset image size

Finally, we perform hierarchical inference with one formation channel excluded, such that, while all five channels are contributing sources, the inference is only performed with four channels. This analysis is motivated by the fact that any population analysis done on real BBH data likely does its inference with an incomplete set of formation channels; we most likely do not know the totality of all possible BBH formation channels, nor can we model them all accurately and self-consistently. We again consider an equal-mixture branching fraction universe (β = 0.2 between all formation channels) and true values for our physical prescription of α_CE = 1.0 and χ_b = 0.2. Figure 9 shows the marginalized branching fraction posteriors and support for different values of χ_b for the full inference as well as for inferences with one channel excluded. We defer discussion of biases in α_CE selection to Appendix B and instead focus on χ_b and the formation channel branching fractions in this section.

Zoom In Zoom Out Reset image size

By examining which channels receive more or less posterior support as a result of the inference's incomplete knowledge of formation channels, we can infer the correlations between the branching fractions of different formation channels. For example, when the CE channel is excluded (second row), the increase in branching fraction is spread approximately equally over the other channels, suggesting that β_CE is negatively correlated with the other branching fractions, as is the prior. On the other hand, when the SMT channel is excluded from the inference (last row), the β_CE posterior shifts toward higher values, while support for the other channels decreases.

Such correlations are consistent with the branching fraction corner plot from the full inference (Figure 8); the pink contours are relevant to the set of mock detections discussed in this section. The uncertainties in β_CE are the greatest of all channels, and the posterior the least constrained; as such, the marginalized β_CE posterior closely follows the prior. As a result, β_CE is negatively correlated with all channels, which leads to systematic overestimation when doing inference with a channel excluded. In general, we can see that some branching fraction posteriors (i.e., β_NSC and β_CHE) are constrained much better than others and are affected the least by the choice of prior. As seen in Section 4.1, the prior can cause a bias toward higher values of β for channels with lower detection efficiencies and greater uncertainties.

Next, we highlight the effect of incomplete formation channel knowledge on the selection of the natal black hole spin χ_b. We see two cases in which the wrong value of χ_b is selected. When the CHE channel is excluded (third row), we infer a high natal black hole spin of χ_b = 0.5, as indicated by the yellow curve, with a Bayes factor over the true value χ_b = 0.2 of ${\mathrm{BF}}_{{\chi }_{{\rm{b}}}=0.2}^{{\chi }_{{\rm{b}}}=0.5}=373$. This is due to the high-χ_eff BBHs that the CHE channel produces. When the inference does not account for the CHE channel, it tries to explain these highly spinning CHE black holes with other field-channel BBHs spinning at a higher χ_b. Then, in order to still account for the lower χ_eff non-CHE BBHs, the branching fractions for the GC and SMT channels are increased and decreased, respectively. We remark that no such adjustment is seen for the NSC and CE channels, despite them having features in χ_eff space similar to those of the GC and SMT channels, respectively.

An opposite effect is seen when excluding the NSC channel (fifth row), which, as noted above, produces BBHs scattered around χ_eff = 0 with tails that extend to more positive and negative values, due to the presence of hierarchical mergers. The inference has constructed two competing explanations in order to explain these lower-spinning BBHs: low-χ_b CE BBHs (represented by the green and blue curves with posterior support for high β_CE and low β_GC) and higher-χ_b GC BBHs (represented by the pink curve with support for low β_CE and high β_GC). This highlights the degeneracy between low-spinning field channels and high-spinning dynamical channels in producing similar features in the population χ_eff distribution. This case is especially remarkable because the green and blue curves (χ_b = 0.0 and χ_b = 0.1 ) for the β_CE posterior bear a striking resemblance to the β_CE posterior inferred from GWTC-3 data (Figure 2), even though these features are purely an artifact of the inference neglecting a single BBH formation channel. We again remark on the difference between the β_CE and β_SMT posteriors: despite the fact that both CE and SMT are field channels with the similar features in χ_eff space, support for β_SMT does not increase (and rather decreases) in the low-χ_b model as a result of the exclusion of the NSC channel. One reason why this may occur is that SMT BBHs cannot go through tidal spin-up, unlike CE BBHs, and hence they account for a narrower range of χ_eff concentrated around χ_b. Therefore, β_SMT can be strongly affected by the exclusion of a channel with strong features in χ_eff space, especially when the wrong model of χ_b is inferred.

As mentioned in the previous paragraph, we also find it noteworthy that the inference favors two separate competing explanations when NSC is excluded. The GC and NSC marginalized χ_eff KDEs are similar due to the isotropy of spin orientations in these two channels; upon examination of the marginalized KDEs of the other three BBH parameters (${{ \mathcal M }}_{{\rm{c}}}$, q, z), the GC channel appears still to have the closest resemblance to the NSC channel. Despite this, the inference does not simply overcompensate for the exclusion of the NSC channel by correspondingly increasing β_GC, suggesting the influence of higher-dimensional features in parameter space. On the other hand, no such effect is seen when excluding the GC channel, whose posterior has a simple overcompensation in β_NSC and β_CE.

Indeed, although we have been able to broadly interpret these posteriors by focusing on different channels' features in χ_eff space, there must be other subtle effects at play. We have shown in multiple ways the differences between the β_CE and β_SMT posteriors and the β_GC and β_NSC posteriors, despite having them similar features in χ_eff space. Although the exclusion of each channel has its own unique and interesting consequences, there appears to be a bias for the CE and GC channels, systematically underestimating the CHE and NSC channels. We point again to detection efficiency: of the field and dynamical channels, respectively, the CE and GC channels have by far the lowest detection efficiencies. With the uniform branching fraction spread in our current set of hyperparameters, only 2% of detections are expected to be from the CE channel (versus 30% and 16% from the CHE and SMT channels) and 16% from the GC channel (versus 36% from the NSC channel).

Finally, we have repeated this analysis for the set of hyperparameters used in the convergence analysis in Section 3, with an unequal mixture of channels and χ_b = 0.0. In this universe, the bias toward overestimating β_GC is more severe due to the lack of spin information from our choice of χ_b. Consistently with our previous discussion, there again seems to be a bias toward the CE and GC channels when one channel is excluded. Due to the CHE and NSC channels' high detection efficiencies, the β_CHE and β_NSC posteriors support their low true values (5% for both channels) with small uncertainties. The corresponding plot of the branching fraction posteriors for this universe can be found in Appendix B.

In this section, we conducted several investigations of different biases that arise in hierarchical Bayesian inference based on astrophysical formation models of BBHs. We summarize the key takeaways as follows:

• 1.  
  Single-channel dominated universes (Section 4.1)

  • (a)  
    Even when our models of the underlying formation channels are perfectly accurate, at N_obs = 100 the data are still relatively uninformative due to information loss from parameter estimation and low detection rates. Thus, results of inference are still influenced by the choice of a flat prior and exhibit bias toward a mixture of formation channels. This puts a caveat on our result from Section 2.2 that no single channel dominates the underlying BBH population, from the inference on GWTC-3 data.
  • (b)  
    It is difficult to precisely infer the true value of α_CE, as it only affects the CE channel, which produces very few detectable BBHs, due to its low detection efficiency relative to the other channels considered. Only with mock catalogs where the CE channel dominates the detections can we recover the true value of α_CE, as expected.
• 2.  
  Biases in spin inference (Section 4.2)

  • (a)  
    It is easier to infer higher values of χ_b than lower values. When χ_b is low (0.0 or 0.1), one has less spin information, and uncertainties in both the branching fraction posteriors and the selection of the true value of χ_b are greater.
• 3.  
  Inference with incomplete populations (Section 4.3)

  • (a)  
    Both the branching fraction posteriors and the inferred value of χ_b can be heavily affected if the inference is performed without knowledge of all contributing formation channels. Some branching fractions are overestimated or underestimated by a factor of ∼3 or more from the exclusion of a formation channel, and ignoring channels that produce particularly high or low χ_eff BBHs can cause the inference to strongly select an incorrect value of χ_b.
  • (b)  
    Degeneracies exist between different sets of hyperparameters that can make it difficult for the inference to discriminate between them. For example, it can be difficult to distinguish between field BBHs with low χ_b and dynamical BBHs with higher χ_b, due to χ_eff being the only spin information used in the inference (which in turn is due to the fact that χ_eff is arguably the only spin parameter that can be measured for all BBHs with advanced detectors). There are likely more subtle degeneracies and correlations in higher-dimensional parameter space that are difficult to explain from the marginalized distributions but nonetheless play a role in the inference.
  • (c)  
    Inference on the underlying branching fractions can be biased due to the varying detection efficiencies of different channels. In particular, the CE channel has a detection efficiency nearly an order of magnitude below the other channels, causing large measurement uncertainties in β_CE and a relatively uninformed (i.e., close to the prior) posterior for β_CE. As a result, the exclusion of a channel from the inference usually results in the β_CE posterior support extending to higher values; similar effects can be seen in other low detection efficiency channels, such as the GC and SMT channels.