Two of a Kind: Comparing Big and Small Black Holes in Binaries with Gravitational Waves

When modeling the population of merging binary black holes, analyses have generally focused on characterizing the distribution of primary (i.e., more massive) black holes in the binary, while using simplistic prescriptions for the distribution of secondary masses. However, the secondary mass distribution and its relationship to the primary mass distribution provide a fundamental observational constraint on the formation history of coalescing binary black holes. If both black holes experience similar stellar evolutionary processes prior to collapse, as might be expected in dynamical formation channels, the primary and secondary mass distributions would show similar features. If they follow distinct evolutionary pathways (for example, due to binary interactions that break symmetry between the initially more massive and less massive stars), their mass distributions may differ. We present the first analysis of the binary black hole population that explicitly fits for the secondary mass distribution. We find that the data is consistent with a ∼30 M ⊙ peak existing only in the distribution of secondary rather than primary masses. This would have major implications for our understanding of the formation of these binaries. Alternatively, the data is consistent with the peak existing in both component mass distributions, a possibility not included in most previous studies. In either case, the peak is observed at 31.4−2.6+2.3M⊙ , which is shifted lower than the value obtained in previous analyses of the marginal primary mass distribution, placing this feature in further tension with expectations from a pulsational pair-instability supernova pileup.

1. INTRODUCTION Dozens of gravitational-wave (GW) events have been observed by the LIGO (Aasi et al. 2015), Virgo (Acernese et al. 2014) and KAGRA (Akutsu et al. 2021) detector network (Abbott et al. 2021a), and many more detections are anticipated by the end of the fourth observing run.However, the formation mechanism of binary black hole (BBH) mergers, which source the majority of detected GW events, is still largely uncertain.While it is not possible to know the formation history of any one BBH with certainty, the population of all merging BBHs encodes information about which astroafarah@uchicago.edu physical processes give rise to the bulk of these systems (e.g.Stevenson et al. 2015;Zevin et al. 2017).
The stellar-mass BBHs detectable by LIGO, Virgo and KAGRA are likely created from the collapse of massive stars.These massive stars may be born as a binary system in the galactic field that then evolves into a BBH system.Alternatively, they may be born in a dense stellar environment, such as a star cluster, in which the BH stellar remnants dynamically assemble into tightly-bound binaries.The population of GW sources contains signatures of the initial conditions of their progenitor stellar systems, as well as many of the evolutionary processes that occur between star formation and BBH merger.The initial conditions that impact the GW source population include the initial mass function of either binary or single stars (binary IMF and IMF, respectively) in different environments, the birth metallicities of the progenitor stars, and other aspects of their formation environments.Depending on the BBH formation scenario, the evolutionary processes that are imprinted on the BBH population include stellar mass loss, transfer of matter between two component stars, the supernova process for massive stars, and dynamical interactions in star clusters or the disks of active galactic nuclei (see reviews by Mapelli 2020;Mandel & Farmer 2022, and references therein).
In general, these uncertain formation processes affect the masses, spins, redshifts, eccentricities, and merger rates of BBH systems.Here we focus on the mass distribution, since BH masses are well-measured from the GW signal.
The BBH mass distribution is typically parameterized by the primary mass m 1 , the larger of the two component masses in the binary, and either the secondary mass m 2 , or the mass ratio q = m 2 /m 1 .Within these parametrizations, it is typically assumed either that the primary and secondary masses follow the same underlying distribution (Fishbach & Holz 2020;Doctor et al. 2020;Farah et al. 2022;Edelman et al. 2022b;Abbott et al. 2023a;Sadiq et al. 2023), or that the primary mass distribution has distinct features from the secondary mass distribution.The latter assumption implies that primary and secondary mass are physically meaningful labels, and is typically achieved by modeling the secondary mass distribution as a single power law between some minimum mass and m 1 (Fishbach & Holz 2017a;Kovetz et al. 2017;Talbot & Thrane 2018;Baxter et al. 2021;Abbott et al. 2021b;Tiwari 2021;Abbott et al. 2023a;Edelman et al. 2022a;Callister & Farr 2023;Godfrey et al. 2023).The former assumption-that m 1 and m 2 follow the same underlying distribution-implies that the two-dimensional mass distribution is symmetric in primary and secondary mass, meaning that m 1 and m 2 can be interchanged without changing the mass distribution.Neither of these assumptions have been explicitly verified with the available data.In this work, we aim to determine if both components in merging BBHs follow the same underlying distribution or if there is a physical distinction between primary and secondary masses.
Understanding whether the primary and secondary masses in BBH mergers follow the same distribution will provide insight into their formation histories.For BBH mergers that are dynamically assembled in dense environments, we expect that both component BHs are drawn from the same population of stellar remnants (i.e.m 1 and m 2 are not physically meaningful labels), and that the two-dimensional BBH mass distribution will therefore be symmetric in m 1 and m 2 .For BBH mergers that originate from binary stars that formed and evolved in relative isolation ("field binaries"), there may be a physical distinction between primary and secondary BH masses.To start off, the progenitor stars in the binary are drawn from the binary IMF, which may not be symmetric between the two components (e.g.Grudić et al. 2023).Furthermore, the two stars exchange mass during binary stellar evolution.In each phase of mass transfer, one component acts as the donor and the other as the accretor, depending on their initial masses.Mass transfer affects the donor and the accretor stars in different ways, which can impact the masses of the resulting BHs following stellar collapse (Laplace et al. 2021).On the other hand, both BH progenitors are expected to have undergone binary stripping and therefore experience similar supernova physics, potentially washing out any significant differences between the mass distributions of firstand second-born BHs (van Son et al. 2022a;Schneider et al. 2023).However, the degree to which these distributions are similar or different depends on uncertain mass loss and accretion physics (van Son et al. 2022b).
Supernova kicks may also be different for first versus second-born components in merging binaries, because kicks determine whether or not the binary can merge within the age of the Universe (Kalogera 1996;Gallegos-Garcia et al. 2022).Because natal kicks are related to the remnant mass via the supernova prescription (Fryer et al. 2012;Mandel et al. 2021), different preferences for the natal kick magnitudes between first-born and second-born BHs may cause the primary and secondary mass distributions to differ in merging binaries formed through isolated binary evolution (e.g.Oh et al. 2023).
In short, binary stellar evolution consists of several processes that can break the symmetry between the population of initially more massive stars, which generally correspond to the first born and more massive (primary) BHs, and the population of initially less massive stars, which generally correspond to the second born and less massive (secondary) BHs.However, if mass inversion occurs in some systems, some initially less massive stars will end up as the more massive BH by the time of merger, and the distribution of primary BH masses will have contributions from both the second-born and firstborn BHs (Olejak & Belczynski 2021;Hu et al. 2022;Broekgaarden et al. 2022;Zevin & Bavera 2022).If mass inversion happens in exactly half of merging BBH systems, the primary and secondary component mass distributions may be indistinguishable even if the first and second born distributions differ.
From a data analysis perspective, knowing that the primary and secondary mass distributions are the same allows us to measure a single set of model parameters: those of the shared distribution.This may allow features of the distribution to be measured with higher precision since both components in the binary will contribute to the measurement of each feature, rather than just the primary mass.Furthermore, disentangling the role of primary and secondary masses aids the physical interpretations of such features.
For example, Abbott et al. (2023a) found that the mass distribution exhibits a peak at primary masses m 1 ∼ 35 M ⊙ .The astrophysical origin of this overdensity in the mass distribution is uncertain, although it may be related to (pulsational) pair-instability supernovae (e.g.Heger & Woosley 2002;Fishbach & Holz 2017b;Talbot & Thrane 2018;Farmer et al. 2019).A necessary ingredient towards understanding the origin of the m 1 ∼ 35 M ⊙ peak is to first understand whether the peak appears exclusively among primary BBH masses, or whether secondary masses also display a peak, indicating that secondary BHs also experience the astrophysical process that leads to a mass pileup.
This paper is organized as follows.In Section 2 we describe the different population models considered in this work.In Section 3 we present the results of fitting each of the models to the Third Gravitational-Wave Transient Catalog (GWTC-3 Abbott et al. 2021a), finding that while the primary and secondary masses of merging BBHs are consistent with following the same underlying distribution, it is also possible that the secondary mass distribution has a more prominent peak than does the primary mass distribution.In Section 4 we discuss possibilities for future observations.In Section 5 we summarize our conclusions and present possible astrophysical interpretations of our results.The appendices discuss the fundamental differences between commonlyused parametrizations for the two-dimensional mass distribution and include details about the population models and statistical framework.

POPULATION MODELS
Because our aim is to learn whether the primary and secondary masses are consistent with being drawn from the same distribution, we describe the primary and secondary mass distributions separately, but according to the same functional form.We model each of the onedimensional mass distributions as a mixture model between a smoothed power law component and a Gaussian component in order to make direct comparisons to the Power Law + Peak model used by the LIGO-Virgo-KAGRA collaboration (LVK) to describe the distribution of primary masses (Talbot & Thrane 2018;Abbott et al. 2021bAbbott et al. , 2023a)).
There are several ways to construct a two-dimensional mass model for the binary from this one-dimensional mass model for the component masses.However, as we discuss in Appendix A, in order to explicitly compare primary and secondary mass distributions, it is necessary to use the pairing function framework first introduced for GW population modeling in Fishbach & Holz (2020).Explicitly, where p 1 (m 1 |Λ 1 ) is the underlying distribution of primary masses, p 2 (m 2 |Λ 2 ) is the underlying distribution of secondary masses, f (q) is a pairing function that depends on the mass ratio of the system1 , Λ 1 and Λ 2 are the hyper-parameters2 describing the underlying primary and secondary mass distributions, respectively, and Λ = {Λ 1 , Λ 2 , β q } is the set of all mass model hyperparameters.In this work, we use a pairing function of the form f (q; β q ) = q βq Θ(q ≤ 1), though other forms may provide a better fit to the data (e.g.Farah et al. 2022).
As illustrated in Figure 1, different choices for the relationship between Λ 1 and Λ 2 result in distinct morphologies for the two-dimensional mass distribution.Below, we list each of the variations we consider in this work, along with the panels in Figure 1 to which they correspond.A table describing these variations in terms of choices for the hyper-prior is given in Appendix B.3.
• Pairing:Symmetric: This model sets Λ 1 = Λ 2 , making the distribution symmetric under the transformation m 1 ↔ m 2 .It corresponds to the first row of Figure 1: any feature in one of the distributions has to be present in both, so it always makes two bands that connect on the diagonal.The different columns in Figure 1 correspond to different power law spectral indices, β q , for the mass ratiodependent pairing function.The leftmost panels show models where components in the binary are allowed to pair up independently of mass ratio, the middle column shows a model where components have a slight preference to pair up with partners that are equal mass, and the rightmost panel shows the case where components are very "picky": they almost always pair up with equal mass partners (Fishbach & Holz 2020).When components pair nearly independently of mass ratio, the asymmetric models produce noticeably different distributions than the symmetric models.However, in the case of very picky binaries, the two scenarios become difficult to tell apart.There is therefore a degeneracy between the steepness of the pairing function and the existence of distinct features in the two mass distributions (see Tiwari 2023, for a discussion of this phenomenon in terms of Jacobian transformations).
In all models considered in this work, the redshift distribution is modeled as a power law with spectral index κ (Fishbach et al. 2018).We use the Default spin model from Abbott et al. (2021bAbbott et al. ( , 2023a) ) to describe the spin magnitudes and tilts of each component.These are the same redshift and spin distributions used with LVK 2023 in Abbott et al. (2023a).The explicit form of the full population model is given in Appendix B.
Using these parameterized models, we construct a hierarchical Bayesian inference (described in Appendix C) to determine the appropriate population-level parameters for the mass distribution, Λ, given the observed set of data {D j } for N observed events (Loredo 2004;Mandel et al. 2019).We model the data as an inhomogeneous Poisson process with the rate density (number of events per unit time per single-event-parameter hypervolume) given by where R acts as a normalizing constant that sets the overall magnitude of the rate.

RESULTS
We fit each model described in Section 2 to the BBHs in GWTC-3.The resulting two-dimensional mass distributions for the Pairing:Symmetric and Pairing:Generic models are shown in Figure 2.These plots represent an average over the hyper-posterior for each model; this average is sometimes referred to as a posterior population distribution (PPD).The contours differ in morphology to those illustrated in Figure 1 because the actual distribution of BBH component masses exhibits two peaks: one at ∼ 10 M ⊙ and another at ∼ 35 M ⊙ (Abbott et al. 2021b;Tiwari 2021;Edelman et al. 2022a;Sadiq et al. 2022;Abbott et al. 2023a;Edelman et al. 2022b;Farah et al. 2023;Ray et al. 2023), whereas we only place one peak in the models shown in Figure 1.The peak at ∼ 10 M ⊙ creates bands in all panels that have very little vertical extent because the peak is proximal to the minimum BH mass.Nonetheless, fits using both models exhibit vertical and horizontal bands, indicating peaks in both p 1 (m 1 ) and p 2 (m 2 ).While the Pairing:Symmetric model requires this behavior, the Pairing:Generic model does not, meaning that the secondary mass distribution appears to exhibit its own feature at ∼ 35 M ⊙ .Furthermore, the fact that Pairing:Generic produces a PPD similar to that of Pairing:Symmetric indicates that the data support consistent primary and secondary mass distributions.
The bands in all panels do not go the full extent of parameter space but rather taper off away from the diagonal, indicating a preference for equal-mass binaries either through a pairing function or a mass ratio distribution that favors m 1 ≈ m 2 .

Primary and secondary masses are consistent with having the same underlying distribution
The underlying distributions (i.e.before a pairing function is applied) of the primary and secondary masses are shown in Figure 3 for the two pairing function models.There is a region of overlap between the primary and secondary mass distributions under the Pairing:Generic model, indicating that the primary and secondary mass distributions are consistent with one another.As expected, this region of overlap also coincides with p(m), the distribution describing both component masses in the Pairing:Symmetric model.
Models not explicitly parameterized in terms of a pairing function are unable to produce underlying distributions such as those shown in Figure 3, so for the sake of comparison to previous work, we turn to con- Pairing:Generic Figure 2. Two-dimensional posterior population distributions (PPDs) for Pairing:Symmetric (left, orange), Pairing:Generic (middle, purple), and Pairing:No Peak in p2(m2) (right, pink ), all fit to the BBHs in GWTC-3.Darker colors indicate a higher rate of sources, and each panel is individually normalized to it maximum value.All three models find a higher rate of events near the m1 = m2 diagonal, as well as peaks at m1 ∼ 10 M⊙ and m1 ∼ 35 M⊙.Pairing:Symmetric and Pairing:Generic both find peaks in m2 as well, as indicated by horizontal bands in the two leftmost panels, whereas Pairing:No Peak in p2( m2) is unable to model features in the m2 direction.Additionally, Pairing:Symmetric and Pairing:Generic display peaks in similar locations despite the fact that Pairing:Symmetric requires that both m1 and m2 follow the same underlying distribution but Pairing:Generic is able to model each component separately.This suggests that both components may be drawn from the same underlying distribution, up to a pairing function.Under the Pairing:Symmetric model, p1(m1) = p2(m2) ≡ p(m), so only p(m) is shown.p(m) under Pairing:Symmetric is more tightly constrained than p1(m1) (purple filled band ) or p2(m2) (dot-dashed lines) under Pairing:Generic as it has twice as many observations per free parameter.These distributions are constructed to describe the population of black holes before the function that pairs components into merging binaries is applied.All three underlying distributions are consistent with one another, though p2(m2) appears to have a large peak at ∼ 35 M⊙ while p1(m1) has some support for no peak in that region.p(m) does not have support for no peak, but its peak is constrained to be small, while the peak in p2(m2) is less well-constrained and may be larger in amplitude.This hints at the possibility that the peak identified in the primary mass distribution by the LVK 2023 formalism may have been driven by a peak in p2(m2) rather than p1(m1), though hyper-parameter uncertainties within Pairing:Generic are too large to definitively determine this, and the data is consistent with p1(m1) = p2(m2).Lines denote the mean posterior probability, marginalized over hyper-parameter uncertainty, and credible intervals are omitted for clarity.When m1 is below the peak, all models agree.When m1 is above or within the peak, Pairing:Symmetric and Pairing:Generic exhibit a larger preference for m2 to be in the peak than the LVK 2023 model does because the latter is constructed to behave like a power law in the range [mmin, m1].The fact that Pairing:Generic and Pairing:Symmetric approximately agree on the peak location and height indicates that the primary and secondary masses may follow the same underlying distribution up to a pairing function.
butions.In particular, under the Pairing:Symmetric and Pairing:Generic models, we see a peak in the conditional secondary mass distribution p(m 2 |m 1 = C) when C ≥ 35 M ⊙ (solid and dotted lines).While m 2 generally prefers to be near m 1 , as indicated by an upward trend in all models, the two pairing function models have more support for m 2 being within the peak than being nearly equal to m 1 when m 1 is large (dotted lines).This behavior is in contrast to the LVK 2023 model, which forces the conditional secondary mass distribution p(m 2 |m 1 = C) to monotonically increase for all values of C because it does not allow for a peak in m 2 .While the Pairing:Generic model can replicate this behavior, it instead recovers similar distributions to the Pairing:Symmetric model.The Bayes factor of Pairing:Generic relative to Pairing:Symmetric is B(Pairing:Generic)/B(Pairing:Symmetric) = 0.38, indicating an inability to distinguish between the two models, with perhaps a modest preference for the Pairing:Symmetric model.

Improved constraints on peak location
If the primary and secondary mass distributions are identical, we can better measure the properties of features in that common distribution by using measurements of both m 1 and m 2 .Figure 5 shows the onedimensional hyper-posterior for the location, µ, of the Gaussian peak in the mass distribution under the Pairing:Symmetric and LVK 2023 models.We measure µ with a standard deviation of 1.50 M ⊙ under the Pairing:Symmetric model compared to 2.08 M ⊙ under the LVK 2023 model, representing a 28% improvement.This increase in precision is similar to that expected from using twice the number of independent events (1 − 1/ √ 2 = 0.29), because now both m 1 and m 2 contribute to the inference, rather than just m 1 .
Furthermore, the Pairing:Symmetric model recovers a lower value of µ = 31.4+2.3 −2.6 M ⊙ (median and 90% credible interval) compared to the LVK 2023 result of µ = 33.6 +2.6 −4.0 , because it refers to a feature in the underlying p(m) distribution rather than the marginal m 1 distribution.Features in p(m) appear at larger masses in the marginal m 1 distribution and lower masses in the marginal m 2 distribution due to the constraint that m 1 > m 2 .
If the ∼ 35 M ⊙ peak is believed to be a feature of the supernova remnant mass distribution, it is best to use its location in the underlying p(m) distribution rather than its location in the marginal distribution.For example, analyses wishing to compare this feature with expectations from a pair-instability supernova pileup (e.g.Stevenson et al. 2019;Farmer et al. 2020) should use the value presented here (µ = 31.4+2.3  −2.6 M ⊙ ).Interestingly, this lower peak location is in further tension with the latest theoretical predictions for a pair-instability pileup (Farag et al. 2022).
When possible, though, it is best to compare theoretical predictions directly to the two-dimensional mass distribution (such as those in Figure 2) rather than to the values of specific hyper-parameters, since hyperparameters have different meanings in different models.Correspondingly, it is important that models used to fit the data are intentionally constructed to allow for features in both the primary and secondary mass distribution.Many of the current parametric (Fishbach & Holz 2017a;Talbot & Thrane 2018;Baxter et al. 2021;Abbott et al. 2021bAbbott et al. , 2023a) ) and non-parametric (Tiwari 2021;Edelman et al. 2022a;Callister & Farr 2023;Godfrey et al. 2023) BBH mass distribution models enforce asymmetry between m 1 and m 2 , excluding the possibility that the primary and secondary mass distributions share the same features4 (with a few exceptions, e.g.Edelman et al. 2022b;Ray et al. 2023;Sadiq et al. 2023).It is possible to construct a mass distribution that allows for features in both the primary and secondary mass distribution without using a pairing function (e.g.Ray et al. 2023), but directly comparing the primary and secondary mass distributions using these parametrizations is less straightforward.
Another application of our improved measurement of the peak location is "spectral siren" cosmology, which uses such features in the mass distribution to infer the expansion history of the universe (Chernoff & Finn 1993;Messenger & Read 2012;Taylor et al. 2012;Farr et al. 2019;Ezquiaga & Holz 2021, 2022;Abbott et al. 2023b).Previous analyses found the location of this peak to be the parameter most correlated with the Hubble constant (see Figures 5 and 13 of Abbott et al. 2023b), so an im-proved constraint on µ should therefore have a relatively large effect on the constraints of cosmological parameters if all other mass distribution parameters remain equally-well constrained.

Evidence for a peak in the secondary mass distribution
Consistency between the primary and secondary mass distributions implies a peak in the secondary mass distribution at m 2 ∼ 30 M ⊙ , as this feature has already been robustly identified in the primary mass distribution (Abbott et al. 2021b;Tiwari 2021;Abbott et al. 2023a;Edelman et al. 2022b;Callister & Farr 2023;Farah et al. 2023).However, the data is also consistent with differing primary and secondary mass distributions, as shown by the regions in which the filled and dashed bands do not overlap in Figure 3.In this case, it is worthwhile to explicitly determine whether there exists a peak in the secondary mass distribution.
The left panel of Figure 6 shows the posterior distribution for the parameters governing the height of the peak in m 1 and m 2 under the Pairing:Generic model, λ 1 and λ 2 , respectively.Setting λ {1,2} = 0 means that no binaries in the astrophysical population are in the Gaussian peak5 , while λ {1,2} = 1 means all binaries are in the Gaussian peak.The lower left and upper right regions of the two-dimensional posterior are excluded, meaning that λ 1 and λ 2 cannot both be zero, nor can they both be large.This indicates that either one of the two distributions has a large peak while the other has none, or that both distributions have moderate peaks.
In fact, the posterior on λ 2 peaks at a higher value than λ 1 .Under the LVK 2023 model, the 1D posterior on λ peaks between the marginal λ 1 and λ 2 posteriors of the Pairing:Generic model.Additionally, in the Pairing:Generic model, λ 1 has more support at zero relative to λ, while λ 2 has reduced support at zero.This hints at the intriguing possibility that the secondary masses may be driving the nonzero value of λ found by Abbott et al. (2021b) and Abbott et al. (2023a).In other words, it is possible that the ∼ 30 M ⊙ peak exists the secondary mass distribution rather than the primary mass distribution.However, the data are still consistent with λ 1 = λ 2 : the dashed grey line in Figure 6 intersects the 1.5-σ contour of the hyper-posterior.
As shown in Appendix A.3, the LVK 2023 formalism behaves similarly to Pairing:No Peak in p 2 (m 2 ), which is nested within Pairing:Generic when λ 2 = 0. Therefore, ruling out λ 2 = 0 would indicate that the pairing function formalism is strongly preferred.We find that λ 2 = 0 is disfavored but not ruled out: there is 4.2 times more posterior density at λ 2 = 0.17 (its median a posteriori value) than at λ 2 = 0.It is therefore difficult to tell with the current number of detections whether the data prefer for the peak to exist in p 1 (m 1 ), p 2 (m 2 ), or both.
These three possibilities are somewhat degenerate because of the preference for equal mass BBHs.To illustrate this, the right panel of Figure 6 compares the power law spectral index, β q , of the pairing function under Pairing:Generic and Pairing:No Peak in p 2 (m 2 ).When there is no peak in p 2 (m 2 ), the β q posterior shifts to higher values, which corresponds to a larger preference for q ∼ 1.This is because the Pairing:No Peak in p 2 (m 2 ) model has a peak in p 1 (m 1 ), so it allows for a high fraction of secondary masses to lie within the peak by making more binaries equal-mass.For comparison, we also plot the posterior on β q under Pairing:Symmetric to show how a different set of assumptions about λ 2 changes β q .The shift in β q when λ 2 = 0 is larger than when λ 2 = λ 1 , suggesting that it is driven by an excess of events with m 2 ∼ 35 M ⊙ rather than by other model assumptions or a different realization of the inference.
In summary, we find modest evidence for a peak in p 2 (m 2 ), suggesting that pairing function models may be preferred over the LVK 2023 formalism, though current observations do not allow us to definitively rule out that the peak exists only in p 1 (m 1 ).The data is consistent with Λ 1 = Λ 2 , so there is not strong evidence against the possibility that m 1 and m 2 are drawn from the same underlying distribution, up to a pairing function.

Binary black holes are picky
Pairing functions provide an intuitive way to quantify whether the properties of one black hole in a binary influences those of its companion.Binaries that pair up independently of each component's masses are described by a pairing function with β q = 0, whereas a preference for equal-mass binaries is described by β q > 0. We find β q > 0 to > 99.99% for Pairing:Symmetric and to 99.95% for Pairing:Generic.This is consistent with an earlier study by Fishbach & Holz (2020) who use the pairing function formalism on GWTC-1 (Abbott et al. 2019a).However, BBHs are not maximally picky: the posterior on β q is not railing against the high end of the prior bounds. of BBHs have secondary masses in the Gaussian peak while 4% have primary masses in the peak, though this preference for a larger peak in p2(m2) may be due to the relatively poor constraint on λ2 rather than a true preference in the data.The lower left portion of the two-dimensional posterior is excluded, indicating that there must be a peak in either p1(m1) or p2(m2), and a slight preference exists for the peak to be in p2(m2) (upper left portion of plot) or both distributions (central portion of plot), though the peak being in p1(m1) only is not completely ruled out.For reference, the 1D posterior on λ under the LVK 2023 model is shown in green, with λ = 0.04 +0.03 −0.02 (Abbott et al. 2023a).The dashed grey line indicates where λ1 = λ2.Right: Marginal posteriors on βq, the hyper-parameter controlling the steepness of the pairing function under the Pairing:Symmetric, Pairing:Generic, and Pairing:No Peak in p2(m2) models.When we set λ2 = 0, βq increases, indicating a preference for equal-mass components.This behavior is likely caused to accommodate an excess of events with m2 ∼ 35 M⊙, which can either be caused by a peak in p2(m2) or with a peak in p1(m1) and a strong preference for equal-mass binaries (see discussion in Appendix A.1).
Our inferred two-dimensional mass distribution is generally consistent with that shown in Sadiq et al. (2023), who use a non-parametric approach and assume that the population is symmetric under m 1 ↔ m 2 .They show evidence for peaks in the mass ratio distribution at q ∼ 0.5 when m 1 ∼ 15 M ⊙ and m 1 ∼ 70 M ⊙ , which they interpret as a lack of preference for similar-mass pairings.However, we note that these peaks in the mass ratio distribution translate to peaks in the secondary mass distribution at ∼ 7 M ⊙ and ∼ 35 M ⊙ , the same locations at which the primary mass distribution exhibits peaks.Parametrizations that assume symmetry under m 1 ↔ m 2 but do not use a pairing function are unable to disentangle the effects of preference for similar-mass pairings from features in one or both component mass distributions.We therefore conclude that both our results and those presented in Sadiq et al. (2023) are con-sistent with a preference for similar-mass pairings as well as structure in the secondary mass distribution.
In this work, we have only examined power law forms for the pairing function, but more complex models are under investigation and will be presented in future work.If the pairing function is a more complex function of mass ratio, or if it depends on multiple parameters, such as primary mass or spin (e.g.Farah et al. 2022), it may become easier to distinguish between pairing function models and models parameterized similarly to LVK 2023, as complex pairing functions may find more support away from the m 1 = m 2 diagonal.The same is true if the marginal mass ratio distribution has multiple modes, or if mass ratio is allowed to correlate with other parameters (e.g.effective spin, Callister et al. 2021).
The LVK 2023 and pairing function models differ most in their predictions for the number of events in the region m 1 ∈ [35, 40] ∩ m 2 ∈ [m min , 30], because they disagree on whether to place m 2 in the Gaussian peak or somewhat evenly throughout the available parameter space.We can therefore determine which model will be preferred in the future by counting the number of detected events in that region.
It is expected that 260 +330 −150 BBHs will be detected by the end of the LVK's fourth observing run (O4), and 870 +1100 −480 BBHs will be detected by the end of the fifth observing run (O5) (Kiendrebeogo et al. 2023).With 260 (870) total BBH detections, we expect 23.6 ± 4.6 (79.0 ± 8.5) BBHs to fall in the region 30] under the LVK 2023 model and 12.6 ± 3.5 (42.3 ± 6.3) under the Pairing:Symmetric model.This means that by the end of O4, we will be able to distinguish the LVK 2023 and Pairing:Symmetric models to > 2σ, and we will be able to distinguish them to > 4σ by the end of O5.
Of course, it will be necessary to perform a full hierarchical analysis, as measurement uncertainties of detected systems will cause events to scatter into and out of this region.Additionally, in a hierarchical Bayesian context, posterior predictive checks performed on the event-level parameters (such as the true masses of individual events) are less sensitive than those performed on the observed data, such as the strain or the maximum likelihood event parameters (Sinharay & Stern 2003;Gelman 2006;Bayarri & Castellanos 2007;Loredo 2013).This is why, for example, Fishbach et al. (2020) performs posterior predictive checks on the maximum-likelihood values of the observed and predicted events rather than the posterior distributions of those events.Therefore, more sensitive posterior predictive checks may be able to distinguish between the two frameworks with fewer events than we project here.

SUMMARY AND IMPLICATIONS FOR ASTROPHYSICS OF MERGING BBHS
We present the first analysis of the secondary mass distribution of merging BBHs.This allows us to explore whether the primary and secondary component masses in merging BBHs are drawn from the same distribution, or whether the BBH mass distribution is asymmetric in m 1 ↔ m 2 as is commonly assumed in BBH population studies.We find the data to be consistent with two possibilities: either the primary and secondary mass distributions are similar, or a larger peak exists in the secondary mass distribution than the primary mass distribution.In either scenario, a peak likely exists in the secondary mass distribution.This possibility is not con-sidered in many previous analyses of the BBH population, which fix the secondary mass distribution to be a power law (e.g.Abbott et al. 2019bAbbott et al. , 2021bAbbott et al. , 2023a) ) or assume that m 1 and m 2 are interchangeable (e.g.Fishbach & Holz 2020;Sadiq et al. 2023).The existence of this secondary-mass peak has implications for the formation channels of merging BBHs and the origin of the ∼ 35 M ⊙ peak.
If the mass distribution is indeed symmetric under m 1 ↔ m 2 as in our preferred model Pairing:Symmetric, this may imply that a large fraction of merging BBHs are formed through dynamical assembly, in which the two component BHs are born through a similar process and then find each other in a dense stellar environment.In this case, the pairing function and its dependence on mass ratio and/or total mass may encode valuable information about dynamical processes like mass segregation and binary exchanges.
The appearance of a peak in both primary and secondary mass distributions is consistent with an origin in the BH remnant mass function.However, the peak location at 31.4 +2.3 −2.6 M ⊙ is in tension with predictions for pair-instability supernovae, so this feature is likely caused by another astrophysical process.As another application of our work, we suggest that future spectral siren measurements consider the Pairing:Symmetric model when inferring the Hubble constant, because it provides an improved measurement of the peak location relative to the LVK 2023 model and the peak location is strongly correlated with the inferred value of H 0 .
On the other hand, the data remain consistent with an asymmetric BBH mass distribution p(m 1 , m 2 ), in which p 1 (m 1 ) ̸ = p 2 (m 2 ), as in the Pairing:Generic model.This would imply that a significant fraction of merging BBHs originate from field binaries, in which "primary" and "secondary" are physically meaningful labels if they tend to correspond to the first-born or secondborn BH in a binary.The component mass distributions p 1 (m 1 ) and p 2 (m 2 ) would then encode the binary IMF and the highly uncertain physics of binary stellar evolution.Specifically, the location and prominence of features in each component mass distribution may provide insight into how mass loss versus mass accretion affect stellar evolution and collapse, including their effect on supernovae and the BH remnant mass function.For example, when mass transfer is unstable, a peak is expected at ∼ 15 M ⊙ in the primary mass distribution, but near the minimum black hole mass in the secondary mass distribution, though the relative locations of these peaks depends on common-envelope and supernova kick physics (van Son et al. 2022a).
Notably, if the BBH mass distribution is asymmetric, we find that it is possible that no peak exists in the primary mass distribution, and the previously-identified peak in primary mass is actually driven by an overdensity in the secondary mass distribution.A peak that is more prominent in p 2 (m 2 ) than p 1 (m 1 ) may imply that pulsational pair-instability supernovae are more frequent or are shifted to lower masses among second-born BHs.
The degree of asymmetry in the BBH mass distribution provides insight into the frequency of mass inversion in binary stars, providing a complementary probe to BH spins (Mould et al. 2022).If mass inversion never occurs, m 1 will typically be the remnant of the donor star.However, if mass inversion is the norm, m 1 will be the remnant of the accretor star.A perfectly symmetric BBH distribution under isolated binary evolution would imply that mass inversion happens exactly 50% of the time, causing the primary and secondary mass distributions to be indistinguishable even though binary physics imparts different distributions on the first-born versus second-born black hole, but this is statistically unlikely.Determining whether primary and secondary BH masses follow the same underlying distributions is therefore a novel and promising probe of formation channels.However, this test should be interpreted within the full BBH population inference, as theoretical models should make consistent predictions for the mass, spin and redshift distributions.
We have only explored parametric models in this work, and can therefore only comment on the features in the mass distribution explicitly parameterized by our model.Future investigations may find it beneficial to use a nonparametric approach that separately fits the primary and secondary mass distributions by employing a pairing function.
We additionally find that BHs in merging binaries have a strong preference to pair with similar-mass BHs for all forms of the secondary mass distribution considered in this work.This is consistent with the results presented in Fishbach & Holz (2020), though we more confidently exclude the scenario in which BBHs pair in-dependently of mass ratio since we now have many more detected BBHs.
By the end of the LVK's fifth observing run, we expect to confidently distinguish between the different scenarios presented here for the BBH mass distribution.In this work, we found that secondary masses in merging BBH systems likely display a peak at ∼ 35 M ⊙ , whereas previous results identified this peak exclusively among primary masses.With a few hundred additional BBH observations, we expect to determine whether both component mass distributions have a peak at the same location and if the peak is more prominent among secondary masses.
Directly incorporating the distribution of secondary masses can serve as an important tool to constrain formation mechanisms of BBHs.
The authors thank Sharan Banagiri, Thomas Dent, Zoheyr Doctor, Will Farr, Davide Gerosa, Tom Loredo, and Jam Sadiq for helpful discussions.A.M.F. is supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE-1746045.D.E.H is supported by NSF grants AST-2006645 and PHY-2110507, as well as by the Kavli Institute for Cosmological Physics through an endowment from the Kavli Foundation and its founder Fred Kavli.This research has made use of data or software obtained from the Gravitational Wave Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA.This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation.The authors are grateful for computational resources provided by the LIGO Laboratory and supported by National Science Foundation Grants PHY-0757058 and PHY-0823459.

. Differences between mass ratio distributions and pairing functions
Given a form for the primary mass distribution, there are several ways to construct a two-dimensional mass distribution.We discuss two possibilities here that are common in the literature, showing the different effects they have on the resulting two-dimensional mass distribution.
The main qualitative differences between the two parametrizations are illustrated in Figure 7.The top row has three examples of mass distributions that can be described by a model of the form for different values of β.Here, m 1 is the mass of the heavier component in the binary, m 2 is the mass of the lighter component, and q = m 2 /m 1 ≤ 1 is the mass ratio.This can be equivalently written as since m 2 = qm 1 .The parametrization described in Equation A1 is used by all parametric models presented in Abbott et al. (2021b) and Abbott et al. (2023a) that were used to model the primary mass distribution of BBHs, such as Broken Power Law and Power Law + Peak, and is also commonly used in other analyses (e.g., Fishbach & Holz 2017a;Kovetz et al. 2017;Talbot & Thrane 2018;Tiwari 2021;Edelman et al. 2022a;Callister & Farr 2023;Godfrey et al. 2023).For the remainder of this Appendix, we will refer to models parameterized by Equation A1 as "Conditioned-q," since they require the mass ratio distribution to be explicitly conditioned on primary mass.The bottom row of Figure 7 has three examples of mass distributions that can be described using a "pairing function", f (Fishbach & Holz 2020).Models with pairing functions have the form where f (q) is a pairing function that depends on the mass ratio of the system6 , and Λ = {Λ 1 , Λ 2 , β q } is the set of all model hyper-parameters.In this work, we use a pairing function of the form f (q; β q ) = q βq Θ(q ≤ 1), though other forms may provide a better fit to the data (e.g.Farah et al. 2022).In the examples illustrated in Figure 7, the primary and secondary mass distributions are equivalent, so p 1 (m) = p 2 (m) ≡ p(m).Alternatively, Λ 1 = Λ 2 .This describes a situation in which there is a single underlying mass distribution from which both components are drawn.The pairing function then describes how likely the two components are to be combined in a merging binary based on their mass ratio.A pairing function that prefers equal mass binaries will cause a marginal mass ratio distribution that has more support near q = 1, but the inverse is not necessarily true.The parametrization in Equation A3 factorizes the possibility that components in binaries prefer to be near-equal mass and the possibility that the primary and secondary masses have distinct probability distributions.In other words, pairing functions allow us to model the secondary mass separately from the primary mass, while also allowing for the possibility that component BHs prefer to pair with similar-mass BHs.
In Figure 7, a peak at 35 M ⊙ is placed in both models to show the effects of such features in both cases.One consequence of Conditioned-q models is that features such as Gaussian peaks can only appear in the the primary mass distribution.This is shown by the vertical bands in the top row of Figure 7 and lack of horizontal bands, since a band in the vertical (horizontal) direction are caused by a peak or dip in the primary (secondary) mass distribution for a range of secondary (primary) masses.For pairing function models, features can appear in p 1 (m 1 ), p 2 (m 2 ), or both.We have illustrated the case in which the same feature appears in both component mass distributions, and this appears as both the vertical and horizontal bands in the bottom row of Figure 7.If Λ 1 and Λ 2 are allowed to differ, features could appear in only one of these distributions.This would cause there to only be horizontal bands if features only existed in p 2 (m 2 ), and only vertical bands if features only existed in p 1 (m 1 ).Features are also able to appear in different locations in p 1 (m 1 ) vs p 2 (m 2 ) under the pairing function formalism.However, the Conditionedq formalism only allows for bands in the vertical direction, meaning that it is not flexible enough to capture true underlying distributions with features in p 2 (m 2 ).The behavior of the Conditioned-q formalism can in general be approximated by the pairing function formalism, while the opposite is not true.
The different columns in Figure 7 correspond to different power law spectral indices for the mass ratio distribution (top row, β) and mass ratio-dependent pairing function (bottom row, β q ).The top row's leftmost panel has a uniform mass ratio distribution, the middle panel has a mass ratio distribution that mildly favors equal-mass binaries, and the right panel's mass ratio distribution strongly favors equal mass binaries.Analogously, the bottom row's leftmost panel shows a model where components in the binary are allowed to pair up independently of mass ratio, the middle panel shows a model where components have a slight preference to pair up with partners that are equal mass, and the rightmost panel shows the case where components are "picky": they almost always pair up with equal mass partners (Fishbach & Holz 2020).When the mass ratio distribution is broad, or when components pair near-independently of mass ratio, the Conditioned-q models produce noticeably different distributions than the pairing function models.However, in the case of very picky binaries or a steeply rising mass ratio distribution, the Conditioned-q and pairing function models become difficult to tell apart, and likely explain the data equally well.There is therefore a degeneracy between the steepness of the pairing function and the existence of distinct features in the two mass distributions (see Tiwari 2023, for a discussion of this phenomenon in terms of Jacobian transformations).
Fortunately, as we will show in Section 3.4, we measure β ∼ 3.5 and β q ∼ 1, so the data lie somewhere between the middle and rightmost columns. 7This means that differentiating between the two scenarios will be difficult, but possible given enough data.
Mass distributions of compact objects are often visualised through a plot of the marginal component mass distributions.The marginal m 2 distribution is defined as where p(m 1 , m 2 |Λ) can be parameterized in the way described in Equation A1 or Equation A3.However, in all but the top left panel in Figure 7, the marginal secondary mass distribution exhibits a peak between 30 M ⊙ -40 M ⊙ , even though the secondary mass distribution is not able to have any features on its own.This is because features in the primary mass distribution induce features in the marginal secondary mass distribution if equal-mass binaries are preferred: when a binary's primary mass is within the peak, its secondary mass is also likely to be in that peak simply because m 2 ≃ m 1 is preferred.It is therefore difficult to distinguish between these different scenarios by looking at the marginal distributions alone.Instead, we analyze two-dimensional distributions such as the ones illustrated in Figure 7, as well as the secondary mass distributions conditioned at various values of m 1 .The latter can be thought of as one-dimensional slices of the former.It is our goal to determine whether the data prefer models described by pairing functions or by Conditioned-q functions.We also aim to determine if the primary and secondary masses are drawn from the same distribution, and whether we can draw physical insights from features that appear in the primary mass distribution, secondary mass distribution, or both.

A.2. "Pickiness"
Under the Conditioned-q formalism, we cannot determine how BHs in binaries choose their companions, though we can gain some insight from their marginal mass ratio distribution.Conditioned-q parametrizes the mass ratio distribution as a power law with spectral index β, where β > 0 corresponds to mass ratio distributions with more support for similar-mass binaries.Abbott et al. (2023a) find β = 1.1 +1.7 −1.3 , which similarly indicates a preference for equal-mass binaries.Note that the value for β under Conditioned-q is noticeably different than the value found for β q under the pairing function models because the two parameters cause different behaviors in the 2D mass distribution within their respective models.A low value for β does not imply that BBHs are not picky, just that the marginal mass ratio distribution rises slowly.As shown in Appendix B.3, the pairing function models and Conditioned-q model all produce near-identical marginal mass ratio distributions, despite different values for β q and β.Low values for β are partially due to the fact that the low-q end of the marginal mass ratio distribution will always be suppressed because of the existence of a minimum BH mass.This minimum mass makes it impossible to get extreme mass ratios unless m 1 is very large, and the m 1 distribution has very little support for m 1 ⪆ 60 M ⊙ .Therefore, β does not need to be large in order for the marginal mass ratio distribution to strongly disfavor unequal-mass binaries.In fact, the existence of a minimum mass plus a tapering at high m 1 means that even when β < 0, the marginal mass ratio distribution rises towards q = 1.Pairing function β q = 10 Figure 7. Illustration of some possible two-dimensional mass distributions under the commonly used "Conditioned-q" formalism described by Equation A1 (top row ) and the pairing function formalism described by Equation A3 (bottom row ).Overdensities/peaks in the mass distribution appear as darker filled contours in these figures.The Conditioned-q formalism is only able to produce models with features in the m1 distribution, as shown by the vertical bands in the top row, whereas the pairing function formalism can model features in either m1 or m2, or both.The different columns correspond to different power law spectral indices for the mass ratio distribution (top row, β) and mass ratio-dependent pairing function (bottom row, βq).In the case of a steeply rising mass ratio distribution, or if components strongly prefer to pair with nearly equal-mass partners, the Conditioned-q model and the pairing function model become difficult to tell apart and likely explain the data equally well, as shown by the two panels in the leftmost column.The diagonal contours in the middle and right columns are caused by a preference for equal-mass binaries and follow lines of constant mass ratio.

A.3. Mimicking Conditioned-q with a pairing function model
We show that the Pairing:No Peak in p 2 (m 2 ) model approximates the morphology of the LVK 2023 model, which uses the Conditioned-q formalism.
The right two panels of Figure 11 show the two-dimensional PPDs for the LVK 2023 and Pairing:No Peak in p 2 (m 2 ) models.Both exhibit vertical bands and no horizontal bands, and have a similar-magnitude drop in merger rate moving away from the diagonal.
The leftmost panel of Figure 11 shows the conditional m 2 distribution for both models.While the slopes of the two models differ slightly, a peak in m 2 appears in neither.We therefore find Pairing:No Peak in p 2 (m 2 ) to be an appropriate proxy for the behavior of the LVK 2023 model for the purposes of this work.A future measurement ruling out λ 2 = 0 with high confidence would therefore serve as evidence for the pairing function models over LVK 2023-like models.

B.1. Base model
We model the two-dimensional mass distribution by constructing separate one-dimensional distributions for the primary and secondary masses and combining draws from these two distributions according to a pairing function, as in Equation A3.Fishbach & Holz (2020) find the pairing function is most informative when parameterized by the mass ratio of the binary, so we adopt a pairing function described by a power law in mass ratio.
In all models considered in this work, the redshift distribution is modeled as a power law with spectral index κ (Fishbach et al. 2018) such that p(z) ∝ dV c dz We use the Default spin model from Abbott et al. (2021bAbbott et al. ( , 2023a) ) to describe the spin magnitudes and tilts of each component.These are the same redshift and spin distributions used with the Power Law + Peak mass model in the analysis presented in Abbott et al. (2023a).
A fit to the model described above is provided in Appendix B, though we focus on specific variations nested within this more general model for the remainder of this work.

B.2. Fit to general base model
We present results of a fit to the most general form of the pairing function model described in Section B.1 (Equation B5).Here, we do not fix any parameters to be equal between Λ 1 and Λ 2 , and instead infer them separately.
Figure 8 shows the posterior of all mass-related hyper-parameters within this model in the form of a corner plot.Because of the large number of free parameters, several are not well-constrained.Nonetheless, all parameters describing p 1 (m 1 ) are consistent with those describing p 2 (m 2 ).Additionally, strong correlations exist between α 1 , α 2 , and β q , making it difficult to meaningfully constrain all three parameters.We therefore set α 1 = α 2 in the Pairing:Generic model.
Figure 9 shows the underlying distributions, p 1 (m 1 ) and p 2 (m 2 ) inferred by the most general version of the base model.As expected, the constraints on these distributions weaken relative to that of the Pairing:Generic model.However it is still clear that p 1 (m 1 ) and p 2 (m 2 ) are consistent with one another and that it is possibile that p 2 (m 2 ) has a larger peak than p 1 (m 1 ) does.Notably, there now seems to be very little evidence for a peak in p 1 (m 1 ), whereas a clear peak still exists in p 2 (m 2 ).
These results are consistent with what has been presented using the Pairing:Generic model in the main body of this work.

B.3. Model comparison
In this Section, we compare the various models considered in this work, including the Conditioned-q model from LVK 2023.Table 1 lists the hyper-prior choices made to construct each model, along with descriptions of each hyper-parameter.
Figure 10 shows the m 1 , m 2 , and q PPDs, marginalized over all other parameters.These show a high level agreement between the different models, despite their two-dimensional PPDs differing in Figure 2. Notably, the marginal

Figure 1 .
Figure1.Illustration of some possible two-dimensional mass distributions under the models considered in this work.These distributions are not indicative of specific astrophysical predictions, but are instead meant to be illustrative of the morphologies accessible by current models.Overdensities/peaks in the mass distribution appear as darker filled contours in these figures.The different columns correspond to different power law spectral indices for the mass ratio-dependent pairing function, βq.In the case where components strongly prefer to pair with nearly equal-mass partners, it becomes difficult to determine whether a feature appears only in one component mass distribution (as in the Pairing:No Peak in p2(m2) model, bottom row ) or in both (Pairing:Symmetric, top row ).The diagonal contours in the middle and right columns are caused by a preference for equal-mass binaries and follow lines of constant mass ratio.The goal of this paper is to distinguish between the different scenarios illustrated in this figure.

Figure 3 .
Figure 3. Underlying distribution of primary and secondary masses for the Pairing:Symmetric (orange) and Pairing:Generic (purple) models.Under the Pairing:Symmetric model, p1(m1) = p2(m2) ≡ p(m), so only p(m) is shown.p(m) under Pairing:Symmetric is more tightly constrained than p1(m1) (purple filled band ) or p2(m2) (dot-dashed lines) under Pairing:Generic as it has twice as many observations per free parameter.These distributions are constructed to describe the population of black holes before the function that pairs components into merging binaries is applied.All three underlying distributions are consistent with one another, though p2(m2) appears to have a large peak at ∼ 35 M⊙ while p1(m1) has some support for no peak in that region.p(m) does not have support for no peak, but its peak is constrained to be small, while the peak in p2(m2) is less well-constrained and may be larger in amplitude.This hints at the possibility that the peak identified in the primary mass distribution by the LVK 2023 formalism may have been driven by a peak in p2(m2) rather than p1(m1), though hyper-parameter uncertainties within Pairing:Generic are too large to definitively determine this, and the data is consistent with p1(m1) = p2(m2).

Figure 4 .
Figure4.Conditional PPDs for the models considered in this work.We show the secondary mass distributions conditioned on m1 = 20 M⊙ (dashed ), m1 = 35 M⊙ (solid ) and m1 = 70 M⊙ (dotted ) to exemplify when the primary mass is below, inside, and above the Gaussian peak, respectively.Orange lines correspond to Pairing:Symmetric, purple lines to Pairing:Generic, and green lines to LVK 2023.Lines denote the mean posterior probability, marginalized over hyper-parameter uncertainty, and credible intervals are omitted for clarity.When m1 is below the peak, all models agree.When m1 is above or within the peak, Pairing:Symmetric and Pairing:Generic exhibit a larger preference for m2 to be in the peak than the LVK 2023 model does because the latter is constructed to behave like a power law in the range[mmin, m1].The fact that Pairing:Generic and Pairing:Symmetric approximately agree on the peak location and height indicates that the primary and secondary masses may follow the same underlying distribution up to a pairing function.

Figure 5 .
Figure 5. Marginal posteriors on µ, the location of the Gaussian peak, for the LVK 2023 and Pairing:Symmetric models.This parameter is more precisely measured under the Pairing:Symmetric formalism by ∼ 1/√ 2, indicating that m1 and m2 independently contribute to this measurement.µ is the most correlated parameter with the Hubble constant when using "spectral sirens"(Ezquiaga & Holz  2022;Abbott et al. 2023b), so using the Pairing:Symmetric model will presumably improve cosmological measurements.The central values of the two distributions differ because the hyper-parameter has slightly different effects on the resulting 2D mass distribution each model.

Figure 6 .
Figure6.Hyperposteriors under pairing function models.Left: Corner plot of hyper-parameters governing the height of the Gaussian peak for the primary (λ1) and secondary (λ2) mass distributions under the Pairing:Generic model.The medians and 90% credible intervals are shown above their respective marginal distributions.The medians indicate that roughly 17% of BBHs have secondary masses in the Gaussian peak while 4% have primary masses in the peak, though this preference for a larger peak in p2(m2) may be due to the relatively poor constraint on λ2 rather than a true preference in the data.The lower left portion of the two-dimensional posterior is excluded, indicating that there must be a peak in either p1(m1) or p2(m2), and a slight preference exists for the peak to be in p2(m2) (upper left portion of plot) or both distributions (central portion of plot), though the peak being in p1(m1) only is not completely ruled out.For reference, the 1D posterior on λ under the LVK 2023 model is shown in green, with λ = 0.04+0.03−0.02(Abbott et al. 2023a).The dashed grey line indicates where λ1 = λ2.Right: Marginal posteriors on βq, the hyper-parameter controlling the steepness of the pairing function under the Pairing:Symmetric, Pairing:Generic, and Pairing:No Peak in p2(m2) models.When we set λ2 = 0, βq increases, indicating a preference for equal-mass components.This behavior is likely caused to accommodate an excess of events with m2 ∼ 35 M⊙, which can either be caused by a peak in p2(m2) or with a peak in p1(m1) and a strong preference for equal-mass binaries (see discussion in Appendix A.1).

Figure 8 .
Figure 8. Corner plot of all mass-related hyper-posteriors in the most general form of the base model.Power law spectral indecies governing the slope of the primary and secondary mass distributions are degenerate with one another and with the pairing function power law, so only two of the three parameters can be meaningfully constrained at a time.All other parameters that perform the same function in p1(m1) and p2(m2) are consistent between the two underlying distributions.