# Dropping Anchor: Understanding the Populations of Binary Black Holes with
Random and Aligned Spin Orientations

Vishal Baibhav<EMAIL_ADDRESS>Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
1800 Sherman Ave, Evanston, IL 60201, USA  Zoheyr Doctor Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA),
Northwestern University, 1800 Sherman Ave, Evanston, IL 60201, USA  Vicky
Kalogera  Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 1800 Sherman Ave, Evanston, IL
60201, USA  Department of Physics and Astronomy, Northwestern University, 2145
Sheridan Road, Evanston, IL 60208, USA

###### Abstract

The relative spin orientations of black holes (BHs) in binaries encode their
evolutionary history: BHs assembled dynamically should have isotropically
distributed spins, while spins of the BHs originating in the field should be
aligned with the orbital angular momentum. In this article, we introduce a
simple population model for these dynamical and field binaries that uses spin
orientations as an anchor to disentangle these two evolutionary channels. We
then analyze binary BH mergers in the Third Gravitational-Wave Transient
Catalog (GWTC-3) and ask whether BHs from the isotropic-spin population
possess different distributions of mass-ratio, spin magnitudes, or redshifts
from the preferentially-aligned-spin population. We find no compelling
evidence that binary BHs in GWTC-3 have different source-property
distributions depending on their spin alignment, but we do find that the
dynamical and field channels cannot both have mass-ratio distributions that
strongly favor equal masses. We give an example of how this can be used to
provide insights into the various processes that drive these BHs to merge. We
also find that the current detections are insufficient in extracting
differences in spin magnitude or redshift distributions of isotropic and
aligned spin populations.

##  1 Introduction

Gravitational waves offer an unprecedented look into some of the most rare
phenomena in the universe, in particular the inspirals and mergers of pairs of
compact objects. The properties of these merging black holes (BHs) and neutron
stars such as their masses, spins, and redshifts provide clues on how these
mergers come to occur. So far, nearly 100 gravitational waves have been
detected by the Advanced Laser Interferometer Gravitational-Wave Observatory
(LIGO) (Aasi et al., 2015) and Advanced Virgo (Acernese et al., 2015)
detectors. These detections have been reported on by the LIGO-Virgo-KAGRA
Scientific Collaboration (LVK) (Abbott et al., 2021a, b, 2019a) as well as by
other teams (e.g. Nitz et al., 2021; Olsen et al., 2022).

So far, gravitational waves (GWs) have been seen from mergers of two BHs, two
neutron stars, and neutron stars with BHs. The individual GW signals have
shown a range of progenitor masses, mass ratios, spins, and redshifts (Abbott
et al., 2016, 2020a, 2020b), and one even produced a display of
electromagentic counterparts that were observed by other facilities (Abbott et
al., 2017a, b). Beyond inferences on individual systems, the full catalog of
gravitational waves can be leveraged to measure the population distributions
of underlying merger properties in the universe (e.g. Mandel et al., 2019).
These population-level inferences offer even more clues into the histories of
merging compact objects. For example, the distribution of black-hole masses of
the more-massive merger components $m_{1}$ peaks around 5-10 $M_{\odot}$ and
has an additional bump around 35 $M_{\odot}$, and there is a preference for
the secondary masses $m_{2}$ to be near their respective primary’s masses
(Abbott et al., 2021c; Tiwari & Fairhurst, 2021; Tiwari, 2022; Edelman et al.,
2022a; Abbott et al., 2021d; Sadiq et al., 2022). In terms of spins, the
distribution of effective inspiral spin $\chi_{\rm eff}$, the mass-weighted
projection of component spin vectors onto the orbital angular momentum, peaks
near zero, but is definitively biased towards positive values (Abbott et al.,
2021d, c). Furthermore, $\chi_{\rm eff}$ appears to be anti-correlated with
the binary mass ratio $q\equiv m_{2}/m_{1}$ (Callister et al., 2021; Abbott et
al., 2021d; Adamcewicz & Thrane, 2022) with a possible correlation with
redshift as well (Biscoveanu et al., 2022). Rather than modeling effective
inspiral spins, one can also model the spins of individual BHs and their spin
directions. The individual BH dimensionless spins tend to be small $\sim 0.2$,
but the tilt angle distributions are not particularly well-measured when
accounting for systematic uncertainties (Galaudage et al., 2021; Callister et
al., 2022; Vitale et al., 2022; Edelman et al., 2022b; Golomb & Talbot, 2022).
Ultimately, any proposed formation pathways of the detected BH mergers must
agree with these empirical findings.

The two most well-studied formation channels are isolated binary evolution in
the field and dynamical interactions in clusters (see e.g. (Mandel & Farmer,
2022; Mapelli, 2020) for reviews). For isolated binaries, binary interactions
such as stable mass transfer and CE phases may aid in the formation of two BHs
that can merge within a Hubble time, though the exact efficacy of each of
these processes in producing BH mergers is still under debate (see Gallegos-
Garcia et al. (2021) and references therein). Alternatively, dynamical
channels predict that binary BHs (BBHs) form and harden through three-body
encounters in dense stellar clusters. Other scenarios for the formation and
merger of BBHs include chemically homogenous evolution of isolated binaries
(Marchant et al., 2016; de Mink & Mandel, 2016), AGN disks (Leigh et al.,
2018; Stone et al., 2017; Bartos et al., 2017), secular interactions in
triples (Silsbee & Tremaine, 2017; Hoang et al., 2018; Fragione et al., 2019),
and primordial BHs (Bird et al., 2016). Different formation pathways leave
different imprints on the properties of the BBH population, including the
binary masses, spins, eccentricities, and redshift evolution. Measuring these
distributions informs us on the environment in which BBHs form and evolve
(Zevin et al., 2017; Taylor & Gerosa, 2018; Wysocki et al., 2019; Roulet &
Zaldarriaga, 2019; Abbott et al., 2019b).

One of the most promising signatures is the distribution of BH spin
orientations: systems formed through dynamical interactions are expected to
have isotropic spin orientations, whereas binaries born in the field are more
likely to have spins aligned with the orbital angular momentum (Gerosa et al.,
2013; Rodriguez et al., 2016; Farr et al., 2018; Gerosa et al., 2018;
Stevenson et al., 2017; Vitale et al., 2017b). Using GW observations, we can
separate populations of isotropic and aligned binaries. In this work, we ask
the question: Can we observe differences between the population of BHs with
isotropic spins and the population with aligned spins? If the properties of
systems from the isotropically-distributed spin directions are different from
aligned-spin systems, we could use that as an anchor to study the evolutionary
histories of BHs coming from cluster and field formation channels.

This work is organized as follows: we first introduce a generic mixture-model
framework for separating isolated and cluster binary BHs, using spin tilts as
an anchor. We then apply this framework to investigate whether the mass ratio
distributions are different between these two formation channels. Section 3.1
illustrates how these results can be compared back to population synthesis
studies. We conclude in Section 4 with some remarks on how the framework
described herein can be employed.

##  2 Mixture models

We employ a hierarchical Bayesian inference framework to measure distributions
of properties of individual BBH mergers (such as masses, spins, redshifts,
etc.) using only gravitational-wave data. The distribution of the individual
binary properties, $\vec{\theta}$ can be parameterized in terms of unknown
population-level hyper-parameters $\vec{\Lambda}$. We wish to infer this
$\vec{\Lambda}$ given a catalog of GW detections. We elaborate on the
hierarchical Bayesian inference framework in Sec. A.

To model potential differences in populations of field and cluster BH-merger
properties, we divide the BBH population into two sub-populations, one
requiring isotropic spins and another with preferentially aligned spins,
respectively. We further separate the binary parameters into two sets:

$\vec{\theta}=\\{\vec{\theta}_{\rm mix},\vec{\theta}_{\rm pure}\\}\,,$

where $\vec{\theta}_{\rm mix}$ come from a mixture of the two subpopulations,
and $\vec{\theta}_{\rm pure}$ are drawn from a common distribution for both
subpopulations. Similarly, the hyperparameters corresponding to
$\\{\vec{\theta}_{\rm mix}$ and $\vec{\theta}_{\rm pure}\\}$ can also be
divided into

$\vec{\Lambda}=\\{\vec{\Lambda}_{\rm mix},\vec{\Lambda}_{\rm pure}\\}\,.$

$\vec{\Lambda}_{\rm mix}$ further consists of

$\vec{\Lambda}_{\rm mix}=\\{\vec{\Lambda}_{I},\vec{\Lambda}_{A}\\}$ (1)

where $\vec{\Lambda}_{I}$ is the set of hyperparameters that parameterize the
binary parameters $\vec{\theta}_{\rm mix}$ of the isotropic subpopulation,
while $\vec{\Lambda}_{A}$ parameterizes $\vec{\theta}_{\rm mix}$ of the
aligned subpopulation.

One approach to differentiate the aligned and isotropic binaries is to look at
their spin directions directly:

$\cos\theta_{i}=\hat{\mathbf{S_{i}}}\cdot\hat{\mathbf{L}}$ (2)

where $\theta_{i}$ ($=1,2$) is the angle between the BH spin $\mathbf{S_{i}}$,
and the orbital angular momentum $\mathbf{L}$. Extending the model introduced
in Talbot & Thrane (2017), we assume these tilts are identically distributed
as a mixture between an isotropic component and a preferentially-aligned
component:

$\displaystyle p_{\rm mix}(\vec{\theta}_{\rm
mix},\cos\theta_{i}|\vec{\Lambda}_{\rm
mix},\zeta,\sigma_{t})=\frac{(1-\zeta)}{4}\;p(\vec{\theta}_{\rm
mix}|\vec{\Lambda}_{I})$ (3) $\displaystyle\hskip
56.9055pt+\zeta\;p(\vec{\theta}_{\rm
mix}|\vec{\Lambda}_{A})\;p(\cos\theta|\sigma_{t})\,.$

Here, $\zeta$ is the mixing fraction of events comprising the aligned
subpopulation, and the remaining fraction $1-\zeta$ is assigned to the
isotropic-spin population. The spin directions of aligned binaries are modeled
as a truncated-Gaussian centered at $\cos\theta=1$

$p(\cos\theta|\sigma_{t})=\prod_{i=1,2}{\cal{N}}_{[-1,1]}(\cos\theta_{i}|1,\sigma_{t})$
(4)

Unfortunately, individual properties of BH spin (such as magnitudes and
directions) can not be measured very accurately with GWs. Instead the most
well-measured spin parameter is the effective inspiral spin

$\chi_{\rm eff}=\frac{\chi_{1}\cos\theta_{1}+q\chi_{2}\cos\theta_{2}}{1+q}$
(5)

which is the leading order spin contribution and a constant at 2PN level
(Damour, 2001; Racine, 2008). Here $q=m_{2}/m_{1}$ is the ratio of BH masses
(where $m_{1}$ is the heavier BH while $m_{2}$ is the lighter BH); $\chi_{1}$
and $\chi_{2}$ are the spins of primary and secondary BHs respectively.
$\chi_{\rm eff}$ has long been proposed as a tool to differentiate between
binaries formed in the isolated channel and those assembled dynamically. Since
orientations of BHs formed in stellar clusters are isotropically distributed,
this leads to a symmetric $\chi_{\rm eff}$ distribution centered at $\chi_{\rm
eff}=0$, while the field binaries should be preferentially aligned (Kalogera,
2000), leading to a distribution skewed towards positive values of $\chi_{\rm
eff}$ (Rodriguez et al., 2016; Farr et al., 2018; Wysocki et al., 2018; Gerosa
et al., 2018; Ng et al., 2018). For this reason, another approach to
differentiate population parameters the aligned and isotropic binaries is to
model this mixture as

$\displaystyle p_{\rm mix}(\vec{\theta}_{\rm mix},\chi_{\rm
eff}|\vec{\Lambda}_{\rm mix},\zeta)=$ (6) $\displaystyle\hskip
28.45274pt\zeta\;p(\vec{\theta}_{\rm
mix}|\vec{\Lambda}_{A})\;N_{[-1,1]}(\chi_{\rm eff}|\mu_{\rm eff,A},\sigma_{\rm
eff,A})$ $\displaystyle\hskip 28.45274pt+(1-\zeta)\;p(\vec{\theta}_{\rm
mix}|\vec{\Lambda}_{I})N_{[-1,1]}(\chi_{\rm eff}|\mu_{\rm eff,I}=0,\sigma_{\rm
eff,I})$

where the first term models the distribution of $\chi_{\rm eff}$ for aligned
spin population as a gaussian truncated between $[-1,1]$ with $\mu_{\rm
eff,A}>0$, while the second term represents the isotropic population with a
truncated gaussian centered around $\mu_{\rm eff,I}=0$.

In this study, we will focus on using the spin directions $\cos\theta_{1,2}$
as an indicator of isotropic and aligned populations. We will assume that
priors on $\zeta$ are uniform between $0$ and $1$, while priors on
$\sigma_{t}$ are uniform between $0.1$ and $4$. Our full population prior,
$p(\vec{\theta}|\vec{\Lambda})$, is the product of the mixture [Eqs. 3] and
pure population models:

$p(\vec{\theta}|\vec{\Lambda})=p_{\rm mix}(\vec{\theta}_{\rm
mix},\cos\theta_{i}|\vec{\Lambda}_{\rm mix},\zeta,\sigma_{t})\times
p(\vec{\theta}_{\rm pure}|\vec{\Lambda}_{\rm pure})$ (7)

In this study, we will often parameterize some of the binary properties
$\vec{\theta}$ (either pure $\vec{\theta}_{\rm pure}$ or in the mixture of
aligned and isotropic components $\vec{\theta}_{\rm mix}$). The parameters not
studied will be assumed to come from distributions with hyperparameters fixed
to their median values as inferred from Abbott et al. (2021d).

##  3 Do isotropic and aligned binaries pair differently?

Figure 1: Left panel: posteriors of all the hyperparameters in the q-mix
model: $\beta_{q,I}$, $\beta_{q,A}$, $\zeta$ and $\sigma_{t}$. We plot the
$68\%$ and $90\%$ intervals. Right: we plot the posteriors of
$(\beta_{q,I},\,\beta_{q,A})$. Grey shaded region marks the $90\%$ credible
regions on $\beta_{q}$ recovered in GWTC-3 (Abbott et al., 2021d). We also
show the $\beta_{q,I}=\beta_{q*}$ extracted using the globular cluster models
in Zevin et al. (2021) for different values of BH natal spins $\chi_{b}$
(shown by markers of different shape). Similarly, we use the CE models to
represent aligned binaries with the $\beta_{q,A}=\beta_{q*}$ extracted for
different values of CE efficiencies $\alpha_{\text{\tiny CE}}$(shown by
markers of different colors).

We use the mixture model describe in Sec. 2 to discern if isotropic and
aligned populations have mass pairing. For simplicity, we assume that the
primary-mass distribution of both isotropic and aligned binaries follow the
PowerLaw+Peak model described in (Abbott et al., 2021d) that is the
distribution of primary masses $m_{1}$ is parameterized as the mixture model
of a power law and a gaussian distribution (Talbot & Thrane, 2018). The
primary mass model consists of $7$ parameters: power law slope $\alpha$,
mixture fraction $\lambda_{\rm peak}$, minimum ($m_{\min}$) and maximum BH
mass ($m_{\max}$), mean ($\mu_{m}$) and standard deviation ($\sigma_{m}$) of
the gaussian component and a smoothing factor ($\delta_{m}$). With
$\vec{\Lambda}_{m}=\\{\alpha,m_{\min},m_{\max},\delta_{m},\mu_{m},\sigma_{m},\lambda_{\rm
peak}\\}$

$p(m_{1}|\vec{\Lambda}_{m})=(1-\lambda_{\rm peak}){\mathcal{P}}+\lambda_{\rm
peak}\,{\mathcal{G}}$ (8)

where

${\mathcal{P}}\propto
m_{1}^{-\alpha}S(m_{1},m_{\min},\delta_{m})\Theta(m_{\max}-m_{1})$ (9)

is the power-law distribution and

${\mathcal{G}}\propto\exp\left(-\frac{(m_{1}-\mu_{m})^{2}}{2\sigma_{m}^{2}}\right)S(m_{1},m_{\min},\delta
m).$ (10)

is a Gaussian peak. In addition, this model employs a smoothing function at
the low-masses, $S(m,m_{\min},\delta m)$, which rises from zero to one as mass
increases from $m_{\min}$ to $m_{\min}+\delta_{m}$,

$\displaystyle S(m,m_{\min},\delta_{m})$
$\displaystyle=\frac{1}{1+e^{f(m-m_{\min},\delta_{m})}},$ $\displaystyle
f(m,\delta_{m})$
$\displaystyle=\frac{\delta_{m}}{m}-\frac{\delta_{m}}{m-\delta_{m}}.$ (11)

Since we are only focusing on discerning whether binaries with aligned and
isotropic spins have different mass-ratio distributions, we fix the
hyperparameters governing the primary mass distribution (that is
$\vec{\Lambda}_{m}=\\{\alpha,m_{\min},m_{\max},\delta_{m},\mu_{m},\sigma_{m},\lambda_{\rm
peak}\\}$) to the median values obtained in Abbott et al. (2021d). We allow
the aligned and isotropic binaries to have different mass ratio distributions,
that is, $\vec{\theta}_{\rm mix}={q}$ governed by hyperparameters
$\vec{\Lambda}_{\rm mix}=\\{\beta_{q,I},\beta_{q,A}\\}$. Here $\beta_{q,I}$ is
the power-law slope of the $q$ distribution of isotropic-spin binaries while
$\beta_{q,A}$ is the corresponding power-law slope of aligned-spin binaries
(allowing for a smooth turn-on for low-mass secondary BHs),

$p(q|m_{1},\beta_{q,i},\delta_{m},m_{\min})\propto
q^{\beta_{q,i}}S(m_{2},m_{\min},\delta m)\Theta(m_{1}-m_{2}).$ (12)

where $i=I,A$. We call this model q-mix. We assume priors on $\beta_{q,I}$ and
$\beta_{q,A}$ are uniform between $-2$ and $7$.

Fig. 1 shows results when analyzing GWTC-3 binary BHs with the q-mix model.
The left panel shows the $68\%$ and $90\%$ credible regions for the recovered
power-law slopes of $q$ distribution for isotropic and aligned populations
($\beta_{q,I},\beta_{q,A}$), the fraction of aligned binaries ($\zeta$) and
standard deviation of aligned-spin tilts $\sigma_{t}$. Most notably, part of
$\beta_{q,I}=\beta_{q,A}$ line is included in the $90\%$ credible region,
implying that there is no strong evidence for different $\beta_{q}$ between
the isotropic-spin and aligned-spin populations. Furthermore, the Bayes factor
in favor of different $\beta_{q}$’s is only $1.62$. Nevertheless, we can probe
the rest of the posterior to understand what alternative hypotheses may still
be viable given the data. Firstly, we observe that $\beta_{q,I}$ and
$\beta_{q,A}$ are possibly anti-correlated. Consequently, two regions in
$(\beta_{q,I},\,\beta_{q,A})$ parameter space are disfavored: i) both
$\beta_{q,I}$ and $\beta_{q,A}$ being small ($\lesssim 1$), ii) both
$\beta_{q,I}$ and $\beta_{q,A}$ being large ($\gtrsim 2$). This implies that
both isotropic and aligned binaries cannot have an extremely selective pairing
or random pairing, which comports with the inferred $\beta_{q}\sim 1$ assuming
a single mass-ratio distribution for all BBHs (Abbott et al., 2021d).
Interestingly, there is some support for $\beta_{q,I}$ and $\beta_{q,A}$
having opposite signs. In that case, equal mass-ratio events would be
dominated by one channel, while unequal mass events would be dominated by the
other. Most of the posterior samples have dissimilar $\beta_{q,I}$ and
$\beta_{q,A}$. In addition, the posteriors on $\beta_{q,I}$ are least
informative when $\beta_{q,A}$ is between $0$–$2$. A similar statement can be
made for $\beta_{q,A}$ when $\beta_{q,I}\simeq 1$.

Our results also have implications for the fraction of aligned systems
$\zeta$. If $\beta_{q,I}\simeq 1$, then the data slightly prefer that the
majority of binaries have random spin orientations, but a range of mixing
fractions are possible. On the other hand, if $\beta_{q,A}\simeq 1$, there is
likely a larger contribution from the aligned-spin population. In addition, we
can exclude $\zeta=0$ if $\beta_{q,I}<0$ or $\beta_{q,I}\gtrsim 3$. This is
because if isotropic-spin binaries have unequal masses (equal masses), the
aligned spin channel must be invoked to explain the mergers involving similar
masses (unequal masses). For a similar reason, $\zeta=1$ is disfavored for
$\beta_{q,A}<0$ or $\beta_{q,A}\gtrsim 3$.

There is a positive correlation between $\zeta$ and $\sigma_{t}$ as also
reported in Abbott et al. (2021c): if the tilts of aligned binaries are small,
then their fraction is also small, and most of the BBHs must be isotropic to
explain the observations. If we increase $\sigma_{t}$, the aligned binaries
are allowed to have larger tilts, and binaries that were earlier considered
isotropic are now considered aligned. A majority of the samples identified as
aligned have very large tilts.

As mentioned earlier, BBHs detected during O3 show evidence of anti-
correlation between the mass ratios and spins, with equal masses possessing
smaller $\chi_{\rm eff}$ and unequal-mass mergers exhibiting larger $\chi_{\rm
eff}$ (Callister et al., 2021; Abbott et al., 2021d). There is a possibility
that this behavior can be explained using by assuming isotropic and aligned
binaries have different $q$ distributions. In particular, if isotropic
binaries have a distribution that strongly favors $q=1$ (that is large
$\beta_{q,I}$) while the aligned binaries dominate at smaller mass ratios
(either small $\beta_{q,A}$ or $\beta_{q,A}\leq 0$). In this case, the $q=1$
region would be populated with isotropic binaries with $\chi_{\rm eff}$
symmetric around $0$, while the aligned population will possess smaller mass
ratios and larger values of $\chi_{\rm eff}$. So as we move from $q=1$ to
smaller $q$, the mean of the $\chi_{\rm eff}$ distribution increases as the
presence of aligned binary increases, thereby explaining the $q-\chi_{\rm
eff}$ correlation. Another component that controls the $q-\chi_{\rm eff}$
corrleation is the spin magnitude of individual BHs. We defer the discussion
on spin magntudes of isotropic and aligned systems to Appendix B.

In this section, we allowed the isotropic and aligned populations to have
different $q$ distributions, however, this could also be extended to include
primary masses ($m_{1}$) as well, i.e., $\vec{\theta}_{\rm
mix}=\\{m_{1},q\\}$. In this study, we do not evaluate the differences in the
primary-mass distributions of the two subpopulations because the model becomes
complicated with $16$ hyperparameters ($8$ for each subpopulation), and defer
such investigations to future work with a larger catalog of sources.

Figure 2: Cumulative density function of $q$ distribution of BBHs in the
dynamical (top) and isolated (bottom) scenario. In the top panel (bottom) we
show the CDFs for varying $\chi_{b}$ ($\alpha_{\text{\tiny CE}}$). The dashed
lines show the power-law fit ($p(q)\propto q^{\beta_{q*}}$) to these
astrophysical distributions. Figure 3: The distribution of $\beta_{q}$ for the
LVK PowerLaw+Peak model. Vertical lines mark the $\beta_{q*}$ obtained by
fitting astrophysical models: dot-dashed lines represent the pairing functions
of GC models for $\chi_{b}\in\\{0,0.1,0.2,0.5\\}$ with smaller $\chi_{b}$
yielding smaller $\beta_{q*}$, while the dotted lines represent the pairing
functions of CE models for $\alpha_{\text{\tiny
CE}}\in\\{0.2,0.5,1.0,2.0,5.0\\}$ with smaller $\alpha_{\text{\tiny CE}}$
yielding larger $\beta_{q*}$. The triangle marker represents the $90\%$ upper
limit of the $\beta_{q}$ distribution.

Figure 4: Left panel: posteriors of all the hyperparameters in the q-mix model
when $\beta_{q,I}=3,4$ and $5$: $\beta_{q,A}$, $\zeta$ and $\sigma_{t}$. We
plot the $68\%$ and $90\%$ intervals. Right: the distribution of $\beta_{q,A}$
at the same fixed $\beta_{q,I}$ as the left panel. Vertical lines mark the
$\beta_{q*}$ obtained by fitting CE models for $\alpha_{\text{\tiny
CE}}\in\\{0.2,0.5,1.0,2.0,5.0\\}$ with smaller $\alpha_{\text{\tiny CE}}$
yielding larger $\beta_{q*}$. The triangle marker represents the $90\%$ upper
limit of the $\beta_{q,A}$ distribution.

###  3.1 Astrophysical implications - an illustration

One of the long-standing problems in GW astrophysics is the origin of merging
compact-object binaries and processes that drive mergers. Multiple channels
have been proposed to explain such mergers, some of which could contribute to
the detection of BBHs by LVK. However, categorizing these detections by the
formation pathways is a challenging task, and to complicate matters, these
pathways are plagued by theoretical uncertainties. The task’s difficulty is
greatly alleviated by the tell-tale signatures exhibited by only a subset of
the proposed formation channels. Since the broad differences in spin
directions are immune to astrophysical uncertainties, we use a framework to
employ the spin directions as one such signature to separate multiple sub-
populations. Studying other properties (for example, the mass-pairing function
discussed above) of these sub-populations can further allow us to get a handle
on different unknowns that plague their formation mechanisms.

This section gives an example of how the framework developed here can be
applied to constrain and compare various astrophysical processes across
different formation channels. For illustration purposes, we will use one of
the isolated and dynamical channels from Zevin et al. (2021). Among the three
field formation scenarios in Zevin et al. (2021), we use the late-phase common
envelope (CE) to represent the isolated channel111Note that Gallegos-Garcia et
al. (2021) finds that common-envelope phases are unlikely to lead to BBH
mergers, but we nevertheless make use of common-envelope models for
illustration, since detailed astrophysical modeling in not in the scope of
this article.. For this channel, the POSYDON framework (Fragos et al., 2022)
was used to combine detailed MESA binary simulations with the COSMIC
population synthesis code (Breivik et al., 2020). The mass ratio distribution
of this channel is only governed by the CE efficiency $\alpha_{\text{\tiny
CE}}\in[0.2,0.5,1.0,2.0,5.0]$, with large $\alpha_{\text{\tiny CE}}$ leading
to efficient CE evolution. Small $\alpha_{\text{\tiny CE}}$($\sim 0.2$) tends
to show a preference for BBHs with similar masses. For $\alpha_{\text{\tiny
CE}}>0.2$, we find that the $q$-distribution peaks around $\sim 0.8$, with
larger $\alpha_{\text{\tiny CE}}$ producing more unequal mass mergers. We also
plot the cumulative density function for different $\alpha_{\text{\tiny CE}}$
in Fig. 2. The mass ratio distribution of the common-envelope channel used in
Zevin et al. (2021) does not depend on the natal spin of BHs $\chi_{b}$.

For the dynamical channel, we consider BBH mergers in old, metal-poor globular
clusters (GC). The GC models are taken from a grid of 96 simulated clusters
using the cluster Monte Carlo code CMC (Rodriguez et al., 2019): 24 models
(with a range of initial cluster masses, metallicities, and half-mass radii)
with 4 different natal BH spins ($\chi_{b}\in[0,0.1,0.2,0.5]$). The mass ratio
distribution of BBHs originating in the GCs depends on $\chi_{b}$. This
dependence is primarily due to hierarchical mergers. BHs born from stellar
collapse (henceforth 1g BHs) have a similar $q$-distribution for all
$\chi_{b}$. BHs with small 1g spins impart very small GW recoils to the
remnant BHs (henceforth 2g BHs), which can be retained inside the cluster and
likely merge again with 1g BHs. These 1g+2g mergers have mass ratios peaking
at $\sim 0.5$ (Rodriguez et al., 2019; Kimball et al., 2021, 2020; Kritos et
al., 2022). As the BH natal spin $\chi_{b}$ increases, GW recoils received by
the remnants increase drastically, making their retention difficult. This
significantly reduces the number of BBHs mergers near $q\sim 0.5$ in GCs. We
also plot the cumulative density function for different $\chi_{b}$ in Fig. 2.

Since, in this study, we assume that the mass ratios of both isotropic and
aligned sub-populations are distributed as a power law, we fit the cumulative
density of mass ratios of astrophysical populations used in Zevin et al.
(2021) as

$p(q^{*}<q)=\frac{q^{\beta_{q*}+1}-q_{\rm min}^{\beta_{q*}+1}}{1-q_{\rm
min}^{\beta_{q*}+1}}$ (13)

where $q_{\rm min}$ is the minimum mass ratio in the simulation and
$\beta_{q*}$ is the power law slope of the probability density function. We
plot their cumulative in Fig. 2, and we use the $\beta_{q*}$ as a proxy for
mass-pairing resulting from the CE phase in the field and dynamics in globular
clusters. Here, we stress that, in practice, power laws are not ideal for
fitting most of the astrophysical distributions. This is especially true if a
particular channel also contains multiple sub-populations. For example, in
GCs, the power law fails to appropriately fit the 1g+2g mergers, which
dominate at $q\sim 0.5$. This is important when natal spins are small
($\chi_{b}\simeq 0$) as up to $15\%$ of all GC mergers contain a 2g BH. In
addition, the power laws also fail to fit appropriately when the
q-distribution peaks away from $1$. This is the case for all common-envelope
models (except $\alpha_{\text{\tiny CE}}=0.2$) used in this study.
Nonetheless, for illustration, we assume that $q$-distributions in these
astrophysical models follow a power law.

We can first compare the mass pairing in the astrophysical models with the
results from PowerLaw+Peak model when $q$ distributions are not separated
based on the distribution of their spin directions Abbott et al. (2021d) also
measured the pairing of BBHs $\beta_{q}=$ with a 90% upper limit of $2.38$.
This excludes astrophysical scenarios where all contributing formation
channels produce extremely selective pairing with $q=1$. In Fig. 3, we plot
the recovered $\beta_{q}$ from the LVK analysis, along with $\beta_{q*}$
extracted from the aforementioned individual astrophysical models. Under the
strong assumption that only one channel contributes to the BBH population, we
find that the only models that are not disfavored are common-envelope models
with $\alpha_{\text{\tiny CE}}\geq 0.5$. On the other hand, all the globular
cluster models and CEs models with $\alpha_{\text{\tiny CE}}=0.2$ are unable
to fit the data with as the sole formation channel.

Of course, there is no reason for only one formation channel to contribute to
the BBH population. If BBHs with different spin-direction behavior have
different $q$ distributions, we can allow for both large $\beta_{q}$, which is
typical for dynamical environments, and small $\beta_{q}$, which is typical
for highly efficient CE phases, for example. We observe this in Fig. 1 where
isotropic binaries could mostly prefer equal mass mergers, while aligned
binaries also allow for smaller mass ratios. However, Fig. 1 shows that the
data disfavor the parameter space where both isotropic and aligned channels
strongly prefer $q=1$.

To further study this, we fix $\beta_{q,I}=3,4,5$ as a proxy for the mass
ratios expected in the dynamical scenario. We plot the recovered $\beta_{q,A}$
in Fig. 4. We find the distribution $\beta_{q,A}$ is not very sensitive to the
given values of fixed $\beta_{q,I}$. For $\beta_{q,I}=3$, $4$ and $5$ the
$90\%$ upper limit on $\beta_{q,A}$ is $2.25$, $1.85$ and $1.70$ respectively.
In other words, in these scenarios with a strong mass pairing as expected for
GCs, the aligned channel must have a weaker mass pairing. The results of
detailed modeling can be mapped onto these findings, as exemplified in Figs. 1
and 4 using the CE and GC models of Zevin et al. (2021).

To conclude, here we illustrate how differences in isotropic and isolated
subpopulations can be leveraged to gain insights into astrophysical processes
that govern the merger of BBHs. We caution that in this section’s illustration
we only employed detailed modeling from only two formation channels (one each
for isotropic and aligned spin binaries) out of a slew of scenarios proposed
for synthesizing BBHs. These proposed formation channels are plagued by many
astrophysical uncertainties, affecting population properties or merging
binaries in highly degenerate ways. Consequently, including only a subset of
formation channels can lead to biased inferences on astrophysical
uncertainties (Zevin et al., 2021). While it is not easy to disentangle
formation scenarios with similar spin-direction predictions, we can at least
break degeneracies between populations that predict isotropic and aligned
spins. Using spin directions as an anchor allows us to make more robust claims
than those based on astrophysical models alone.

###  3.2 Caveats

In this study, we used BH spin tilts to indicate the formation environment:
isotropic spins pointing towards the dynamical assembly of BHs and aligned
spins implying an isolated origin. However, GWs contain little information
about the individual BH spins; hence, they are very poorly measured.
Gravitational-wave signals instead depend primarily on the effective spin
parameter $\chi_{\rm eff}$ (defined as the BH spin contribution onto the
binary’s orbital angular momentum) at the leading 2PN order (Damour, 2001;
Racine, 2008). The only scenarios where individual spin components have been
well measured are when the binary is perfectly aligned, when mass ratios are
small or when the binary is observed close to edge-on (Vitale et al., 2017a)
or when the binary is detected with large signal-to-noise ratios (Pürrer et
al., 2016). Since the spin tilts are poorly measured, the estimation of
hyperparameters that govern the properties of isotropic and aligned
populations are also not well measured and depend significantly on the model
employed to describe them (Vitale et al., 2022; Callister et al., 2022; Tong
et al., 2022; Galaudage et al., 2021). Hence, the distribution of
$\beta_{q,I}$ and $\beta_{q,A}$ recovered in this study are not well measured.

In addition, the mixture model used in this study (as well as the Default spin
model in Abbott et al. (2021d)) are plagued by the correlation between $\zeta$
and $\sigma_{t}$. If the distribution of aligned-spin tilts (parameterized by
$\sigma_{t}$) is constrained to smaller angles, then the inferred branching
ratio $\zeta$ is also small, requiring a larger population of isotropic
binaries. On the other hand, if $\sigma_{t}$ is left unconstrained, it shows a
tendency to identify even isotropic binaries as aligned by possessing large
$\sigma_{t}$.

Since $\chi_{\rm eff}$ is measured more accurately than the spin directions
$\theta_{1,2}$, it can also be used to differentiate between binaries formed
in the isolated channel and those assembled dynamically. In this case, the
isotropic binaries will possess a symmetric $\chi_{\rm eff}$ distribution
centered at $\chi_{\rm eff}=0$, while the preferentially-aligned binaries will
have a distribution skewed towards positive values of $\chi_{\rm eff}$
(Rodriguez et al., 2016; Farr et al., 2018; Wysocki et al., 2018; Gerosa et
al., 2018; Ng et al., 2018). If we use the mixture model described in Eq. (6),
we find that there is no evidence of multiple populations. We find that the
observations can be explained with a single distribution consistent with the
Gaussian model employed in Abbott et al. (2021d). This is consistent with
results by Callister et al. (2021). However, it is possible that as the number
of detections increases with technological improvement in GW detectors, one
could discern the differences in isotropic and aligned spin populations using
the $\chi_{\rm eff}$ distribution as well.

Finally, the model used in this work (and the Default spin model in Abbott et
al. (2021d)) assumes that these are the only two categories of spin-direction
distributions: isotropic and aligned. However, formation scenarios – such as
BBHs synthesized in AGN disks – might predict binary where one or both BH
spins are antialigned with orbital angular momentum (Tagawa et al., 2020;
McKernan et al., 2020). We have not considered such populations with
preferentially antialigned spins. Since current observations indicate a dearth
of binaries with $\cos\theta_{i}=-1$ (Abbott et al., 2021d; Galaudage et al.,
2021; Callister et al., 2022; Vitale et al., 2022), this should not
significantly affect our results. But if future observations uncover more
features in the spin distribution of BBHs, our model can be extended to
include an anti-aligned spin component.

##  4 Conclusions

The spin orientations in BH mergers are possibly the cleanest observables to
shed light on BH binary formation: BHs assembled during dynamical encounters
are expected to be isotropic spins, while those formed in isolation should
have spins preferentially aligned with the orbital angular momentum. The
approach to distinguish isotropic and preferentially-aligned binaries has
already been implemented by Abbott et al. (2021d). However, the conventional
model does not provide any astrophysical insights into the two subpopulations
other than their relative abundance or the tilts distribution of the aligned
population (which can shine a light on the supernova kicks that misalign the
BH spins or other processes that realign them (Gerosa et al., 2018)).

In this work, we extend the Default spin model of Abbott et al. (2021d), and
use spin tilts as an anchor to extract more information about the distribution
of binary properties in the isotropic-spin and aligned-spin populations. We
find evidence that BHs coming from the two subpopulations have opposite
tendencies when forming a pair: if BHs with isotropic spins strongly prefer
partners with similar masses, then BHs with aligned spins should be less
picky, and vice versa. We discuss what implications this has on the relative
abundance of isotropic and aligned binaries and their various correlations
with the tilt distribution of the aligned population. We also demonstrate how
differentiating binaries by spin alignment on the population level can provide
insights into the unknown physics and various processes that drive BHs to
merge. For illustration, we compare our results with the mass-pairing function
of BH mergers in globular clusters and those driven by common-envelope (as
presented in Zevin et al. (2021)). We find that the mass-pairing for
isotropic-spin BBHs are consistent with the extremely selective pairing
expected from globular clusters. In addition, our model allows us to put
constraints on astrophysical parameters, such as the common-envelope
efficiency. However, such inferences that employ various astrophysical models
significantly depend on the underlying uncertainties and could be degenerate
with other formation channels not considered in this study. Hence, it is hard
to disentangle individual subpopulations that make up the isotropic/aligned
populations. However, using the prescription outlined in this study, we can
still make broader claims about the relationship between the mass-ratio
distribution of overall isotropic or aligned populations.

We extend our analysis to distributions of spin magnitudes and redshifts in
isotropic and aligned subpopulations in Appendix B. However, the current
observations are insufficient in discerning the differences in spin magnitude
or redshift distribution of isotropic and aligned spin binaries. This could be
due to the fact that individual spin tilts are poorly measured through GWs. As
the current GW detectors undergo further improvements, they will be able to
observe a larger number of sources, enabling us to put tighter constraints on
the contributions of field and dynamical formation channels to the binary
black hole population.

## Acknowledgements

We are grateful to Salvatore Vitale and Sylvia Biscoveanu for the helpful
discussions. This material is based upon work supported by NSF’s LIGO
Laboratory which is a major facility fully funded by the National Science
Foundation. Z.D. acknowledges support from the CIERA Board of Visitors
Research Professorship. V.K. was partially supported through a CIFAR Senior
Fellowship, a Guggenheim Fellowship, the Gordon and Betty Moore Foundation
(grant award GBMF8477), and from Northwestern University, including the Daniel
I. Linzer Distinguished University Professorship fund. This research was
supported in part through the computational resources and staff contributions
provided for the Quest high performance computing facility at Northwestern
University which is jointly supported by the Office of the Provost, the Office
for Research, and Northwestern University Information Technology.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quant. Grav., 32, 074001, doi: 10.1088/0264-9381/32/7/074001
  * Abbott et al. (2016) Abbott, B. P., et al. 2016, Phys. Rev. D, 93, 122003, doi: 10.1103/PhysRevD.93.122003
  * Abbott et al. (2017a) —. 2017a, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Abbott et al. (2017b) —. 2017b, Astrophys. J. Lett., 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2019a) —. 2019a, Phys. Rev. X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2019b) —. 2019b, Astrophys. J. Lett., 882, L24, doi: 10.3847/2041-8213/ab3800
  * Abbott et al. (2020a) Abbott, R., et al. 2020a, Phys. Rev. D, 102, 043015, doi: 10.1103/PhysRevD.102.043015
  * Abbott et al. (2020b) —. 2020b, Astrophys. J. Lett., 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2021a) —. 2021a.  https://arxiv.org/abs/2111.03606
  * Abbott et al. (2021b) —. 2021b, Phys. Rev. X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Abbott et al. (2021c) —. 2021c, Astrophys. J. Lett., 913, L7, doi: 10.3847/2041-8213/abe949
  * Abbott et al. (2021d) —. 2021d.  https://arxiv.org/abs/2111.03634
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quant. Grav., 32, 024001, doi: 10.1088/0264-9381/32/2/024001
  * Adamcewicz & Thrane (2022) Adamcewicz, C., & Thrane, E. 2022, doi: 10.1093/mnras/stac2961
  * Ashton et al. (2019) Ashton, G., et al. 2019, Astrophys. J. Suppl., 241, 27, doi: 10.3847/1538-4365/ab06fc
  * Baibhav et al. (2019) Baibhav, V., Berti, E., Gerosa, D., et al. 2019, Phys. Rev. D, 100, 064060, doi: 10.1103/PhysRevD.100.064060
  * Baibhav et al. (2020) Baibhav, V., Gerosa, D., Berti, E., et al. 2020, Phys. Rev. D, 102, 043002, doi: 10.1103/PhysRevD.102.043002
  * Bartos et al. (2017) Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, Astrophys. J., 835, 165, doi: 10.3847/1538-4357/835/2/165
  * Bavera et al. (2021) Bavera, S. S., et al. 2021, Astron. Astrophys., 647, A153, doi: 10.1051/0004-6361/202039804
  * Bird et al. (2016) Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, Phys. Rev. Lett., 116, 201301, doi: 10.1103/PhysRevLett.116.201301
  * Biscoveanu et al. (2022) Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, Astrophys. J. Lett., 932, L19, doi: 10.3847/2041-8213/ac71a8
  * Breivik et al. (2020) Breivik, K., et al. 2020, Astrophys. J., 898, 71, doi: 10.3847/1538-4357/ab9d85
  * Callister et al. (2021) Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, Astrophys. J. Lett., 922, L5, doi: 10.3847/2041-8213/ac2ccc
  * Callister et al. (2022) Callister, T. A., Miller, S. J., Chatziioannou, K., & Farr, W. M. 2022, Astrophys. J. Lett., 937, L13, doi: 10.3847/2041-8213/ac847e
  * Damour (2001) Damour, T. 2001, Phys. Rev. D, 64, 124013, doi: 10.1103/PhysRevD.64.124013
  * de Mink & Mandel (2016) de Mink, S. E., & Mandel, I. 2016, Mon. Not. Roy. Astron. Soc., 460, 3545, doi: 10.1093/mnras/stw1219
  * du Buisson et al. (2020) du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, Mon. Not. Roy. Astron. Soc., 499, 5941, doi: 10.1093/mnras/staa3225
  * Edelman et al. (2022a) Edelman, B., Doctor, Z., Godfrey, J., & Farr, B. 2022a, Astrophys. J., 924, 101, doi: 10.3847/1538-4357/ac3667
  * Edelman et al. (2022b) Edelman, B., Farr, B., & Doctor, Z. 2022b.  https://arxiv.org/abs/2210.12834
  * Farr et al. (2018) Farr, B., Holz, D. E., & Farr, W. M. 2018, Astrophys. J. Lett., 854, L9, doi: 10.3847/2041-8213/aaaa64
  * Fishbach et al. (2017) Fishbach, M., Holz, D. E., & Farr, B. 2017, Astrophys. J. Lett., 840, L24, doi: 10.3847/2041-8213/aa7045
  * Fishbach et al. (2018) Fishbach, M., Holz, D. E., & Farr, W. M. 2018, Astrophys. J. Lett., 863, L41, doi: 10.3847/2041-8213/aad800
  * Fishbach & Kalogera (2021) Fishbach, M., & Kalogera, V. 2021, Astrophys. J. Lett., 914, L30, doi: 10.3847/2041-8213/ac05c4
  * Fragione et al. (2019) Fragione, G., Grishin, E., Leigh, N. W. C., Perets, H. B., & Perna, R. 2019, Mon. Not. Roy. Astron. Soc., 488, 47, doi: 10.1093/mnras/stz1651
  * Fragos et al. (2022) Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2022, arXiv e-prints, arXiv:2202.05892.  https://arxiv.org/abs/2202.05892
  * Galaudage et al. (2021) Galaudage, S., et al. 2021, Astrophys. J. Lett., 921, L15, doi: 10.3847/2041-8213/ac2f3c
  * Gallegos-Garcia et al. (2021) Gallegos-Garcia, M., Berry, C. P. L., Marchant, P., & Kalogera, V. 2021, Astrophys. J., 922, 110, doi: 10.3847/1538-4357/ac2610
  * Gerosa & Berti (2017) Gerosa, D., & Berti, E. 2017, Phys. Rev. D, 95, 124046, doi: 10.1103/PhysRevD.95.124046
  * Gerosa et al. (2018) Gerosa, D., Berti, E., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 98, 084036, doi: 10.1103/PhysRevD.98.084036
  * Gerosa et al. (2013) Gerosa, D., Kesden, M., Berti, E., O’Shaughnessy, R., & Sperhake, U. 2013, Phys. Rev. D, 87, 104028, doi: 10.1103/PhysRevD.87.104028
  * Golomb & Talbot (2022) Golomb, J., & Talbot, C. 2022.  https://arxiv.org/abs/2210.12287
  * Hoang et al. (2018) Hoang, B.-M., Naoz, S., Kocsis, B., Rasio, F. A., & Dosopoulou, F. 2018, Astrophys. J., 856, 140, doi: 10.3847/1538-4357/aaafce
  * Kalogera (2000) Kalogera, V. 2000, Astrophys. J., 541, 319, doi: 10.1086/309400
  * Kimball et al. (2020) Kimball, C., Talbot, C., L. Berry, C. P., et al. 2020, Astrophys. J., 900, 177, doi: 10.3847/1538-4357/aba518
  * Kimball et al. (2021) Kimball, C., et al. 2021, Astrophys. J. Lett., 915, L35, doi: 10.3847/2041-8213/ac0aef
  * Kovetz et al. (2018) Kovetz, E. D., Cholis, I., Kamionkowski, M., & Silk, J. 2018, Phys. Rev. D, 97, 123003, doi: 10.1103/PhysRevD.97.123003
  * Kritos et al. (2022) Kritos, K., Strokov, V., Baibhav, V., & Berti, E. 2022.  https://arxiv.org/abs/2210.10055
  * Leigh et al. (2018) Leigh, N. W. C., et al. 2018, Mon. Not. Roy. Astron. Soc., 474, 5672, doi: 10.1093/mnras/stx3134
  * Mandel & Farmer (2022) Mandel, I., & Farmer, A. 2022, Phys. Rept., 955, 1, doi: 10.1016/j.physrep.2022.01.003
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, Mon. Not. Roy. Astron. Soc., 486, 1086, doi: 10.1093/mnras/stz896
  * Mapelli (2020) Mapelli, M. 2020, Proc. Int. Sch. Phys. Fermi, 200, 87, doi: 10.3254/ENFI200005
  * Mapelli & Giacobbo (2018) Mapelli, M., & Giacobbo, N. 2018, Mon. Not. Roy. Astron. Soc., 479, 4391, doi: 10.1093/mnras/sty1613
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, Astron. Astrophys., 588, A50, doi: 10.1051/0004-6361/201628133
  * McKernan et al. (2020) McKernan, B., Ford, K. E. S., O’Shaughnessy, R., & Wysocki, D. 2020, Mon. Not. Roy. Astron. Soc., 494, 1203, doi: 10.1093/mnras/staa740
  * Ng et al. (2018) Ng, K. K. Y., Vitale, S., Zimmerman, A., et al. 2018, Phys. Rev. D, 98, 083007, doi: 10.1103/PhysRevD.98.083007
  * Nitz et al. (2021) Nitz, A. H., Kumar, S., Wang, Y.-F., et al. 2021.  https://arxiv.org/abs/2112.06878
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, Astrophys. J., 730, 70, doi: 10.1088/0004-637X/730/2/70
  * Olsen et al. (2022) Olsen, S., Venumadhav, T., Mushkin, J., et al. 2022, Phys. Rev. D, 106, 043009, doi: 10.1103/PhysRevD.106.043009
  * Pürrer et al. (2016) Pürrer, M., Hannam, M., & Ohme, F. 2016, Phys. Rev. D, 93, 084042, doi: 10.1103/PhysRevD.93.084042
  * Racine (2008) Racine, E. 2008, Phys. Rev. D, 78, 044021, doi: 10.1103/PhysRevD.78.044021
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027, doi: 10.1103/PhysRevD.100.043027
  * Rodriguez et al. (2016) Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016, Astrophys. J. Lett., 832, L2, doi: 10.3847/2041-8205/832/1/L2
  * Roulet & Zaldarriaga (2019) Roulet, J., & Zaldarriaga, M. 2019, Mon. Not. Roy. Astron. Soc., 484, 4216, doi: 10.1093/mnras/stz226
  * Sadiq et al. (2022) Sadiq, J., Dent, T., & Wysocki, D. 2022, Phys. Rev. D, 105, 123014, doi: 10.1103/PhysRevD.105.123014
  * Silsbee & Tremaine (2017) Silsbee, K., & Tremaine, S. 2017, Astrophys. J., 836, 39, doi: 10.3847/1538-4357/aa5729
  * Stevenson et al. (2017) Stevenson, S., Berry, C. P. L., & Mandel, I. 2017, Mon. Not. Roy. Astron. Soc., 471, 2801, doi: 10.1093/mnras/stx1764
  * Stone et al. (2017) Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, Mon. Not. Roy. Astron. Soc., 464, 946, doi: 10.1093/mnras/stw2260
  * Tagawa et al. (2020) Tagawa, H., Haiman, Z., Bartos, I., & Kocsis, B. 2020, Astrophys. J., 899, 26, doi: 10.3847/1538-4357/aba2cc
  * Talbot et al. (2019) Talbot, C., Smith, R., Thrane, E., & Poole, G. B. 2019, Phys. Rev. D, 100, 043030, doi: 10.1103/PhysRevD.100.043030
  * Talbot & Thrane (2017) Talbot, C., & Thrane, E. 2017, Phys. Rev. D, 96, 023012, doi: 10.1103/PhysRevD.96.023012
  * Talbot & Thrane (2018) —. 2018, Astrophys. J., 856, 173, doi: 10.3847/1538-4357/aab34c
  * Taylor & Gerosa (2018) Taylor, S. R., & Gerosa, D. 2018, Phys. Rev. D, 98, 083017, doi: 10.1103/PhysRevD.98.083017
  * Tiwari (2022) Tiwari, V. 2022, Astrophys. J., 928, 155, doi: 10.3847/1538-4357/ac589a
  * Tiwari & Fairhurst (2021) Tiwari, V., & Fairhurst, S. 2021, Astrophys. J. Lett., 913, L19, doi: 10.3847/2041-8213/abfbe7
  * Tong et al. (2022) Tong, H., Galaudage, S., & Thrane, E. 2022, Phys. Rev. D, 106, 103019, doi: 10.1103/PhysRevD.106.103019
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Callister, T., et al. 2022, Astrophys. J., 931, 17, doi: 10.3847/1538-4357/ac64a3
  * Vitale et al. (2022) Vitale, S., Biscoveanu, S., & Talbot, C. 2022.  https://arxiv.org/abs/2209.06978
  * Vitale et al. (2022) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2022, in Handbook of Gravitational Wave Astronomy. Edited by C. Bambi, 45, doi: 10.1007/978-981-15-4702-7_45-1
  * Vitale et al. (2017a) Vitale, S., Lynch, R., Raymond, V., et al. 2017a, Phys. Rev. D, 95, 064053, doi: 10.1103/PhysRevD.95.064053
  * Vitale et al. (2017b) Vitale, S., Lynch, R., Sturani, R., & Graff, P. 2017b, Class. Quant. Grav., 34, 03LT01, doi: 10.1088/1361-6382/aa552e
  * Wysocki et al. (2018) Wysocki, D., Gerosa, D., O’Shaughnessy, R., et al. 2018, Phys. Rev. D, 97, 043014, doi: 10.1103/PhysRevD.97.043014
  * Wysocki et al. (2019) Wysocki, D., Lange, J., & O’Shaughnessy, R. 2019, Phys. Rev. D, 100, 043012, doi: 10.1103/PhysRevD.100.043012
  * Zevin et al. (2017) Zevin, M., Pankow, C., Rodriguez, C. L., et al. 2017, Astrophys. J., 846, 82, doi: 10.3847/1538-4357/aa8408
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, Astrophys. J., 910, 152, doi: 10.3847/1538-4357/abe40e

##  Appendix A Hierarchical Bayesian Analysis

We employ a hierarchical Bayesian inference framework to measure the mass,
spin and redshift distributions of BBH mergers using only gravitational-wave
data. We can parameterize the distribution of the individual binary
properties, $\vec{\theta}$ (like masses, spins, redshifts) in terms of unknown
population-level hyper-parameters $\vec{\Lambda}$. We wish to infer this
$\vec{\Lambda}$ given the catalog $D$ ($\equiv\\{d_{i}\\}$) consisting of the
$N_{\mathrm{obs}}=69$ BBHs reported in GWTC-3 with false alarm ratio smaller
than 1 per year (Abbott et al., 2021d). The posterior of the hyperparameters
$\vec{\Lambda}$ governing the distributions of $\vec{\theta}$, (Mandel et al.,
2019; Fishbach et al., 2018; Vitale et al., 2022) is

$p(\vec{\Lambda}|D)\propto{\pi(\vec{\Lambda})}\prod_{i=1}^{N_{\rm
obs}}\frac{p(d_{i}|\vec{\Lambda})}{\xi(\vec{\Lambda})}\,.$ (A1)

where $\pi(\vec{\Lambda})$ is the population prior and
$p(d_{i}|\vec{\Lambda})$ is the likelihood of individual BBHs. Since we are
only interested in the shape of $\vec{\theta}$ distribution, we have
marginalized the overall merger rate density assuming a logarithmically-
uniform prior $p(N_{\mathrm{obs}})\propto N_{\mathrm{obs}}^{-1}$.

For an individual event, the $p(d_{i}|\vec{\Lambda})$ can be expressed as

$\displaystyle p(d_{i}|\vec{\Lambda})=$
$\displaystyle\int\differential\vec{\theta}\ p(d_{i}|\vec{\theta})\
p(\vec{\theta}|\vec{\Lambda})=\int\differential\vec{\theta}\;\frac{p_{\rm
pe}(\vec{\theta}|d_{i})p(\vec{\theta}|\vec{\Lambda})}{p_{\rm
pe}(\vec{\theta})},$ $\displaystyle\simeq$
$\displaystyle\Bigl{<}\frac{p(\vec{\theta}_{i}|\vec{\Lambda})}{p_{\rm
pe}(\vec{\theta}_{i})}\Bigr{>}$

where $p_{\rm pe}(\vec{\theta}|d_{i})$ is the posterior distribution for the
event; $\pi_{\rm pe}(\vec{\theta})$ is the original prior adopted during the
individual-event parameter estimation; and $p(\vec{\theta}|\vec{\Lambda})$ is
the population model for how the hyperparameters $\vec{\Lambda}$ govern the
distribution of binary parameters $\vec{\theta}$. Since we do not have $p_{\rm
pe}(\vec{\theta}|d_{i})$ for each event, but discrete samples drawn from
$p_{\rm pe}(\vec{\theta}|d_{i})$, in the last line, we recast the integral as
Monte Carlo average over the posterior samples.

Detection efficiency $\xi(\vec{\Lambda})$ in the denominator of Eq. (A1) is
the fraction of detectable BBHs given the proposed population $\vec{\Lambda}$.
This correction accounts for observational selection bias and is given by,

$\xi(\vec{\Lambda})=\int
p_{\mathrm{det}}(\vec{\theta})p(\vec{\theta}|\vec{\Lambda})d\vec{\theta},$
(A2)

where $p_{\mathrm{det}}(\lambda)$ is the probability that an event with
properties $\vec{\theta}$ can be recovered by detection pipelines with an
FAR$<1\,\mathrm{yr}^{-1}$.

We can also calculate $\xi(\vec{\Lambda})$ by summing over $N_{\mathrm{det}}$
detectable injections out of $N_{\mathrm{inj}}$ signals drawn from some
reference distribution $p_{\mathrm{inj}}(\vec{\theta})$,

$\xi(\vec{\Lambda})=\frac{1}{N_{\mathrm{inj}}}\sum_{i=1}^{N_{\mathrm{det}}}\frac{p(\vec{\theta}|\vec{\Lambda})}{p_{\mathrm{inj}}(\vec{\theta})}\,,$
(A3)

We use the GWPopulation (Talbot et al., 2019) for applying the hierarchical
bayesian framework, and we use the emcee sampler (implemented in bilby (Ashton
et al., 2019)) to draw samples from the posterior on $\vec{\Lambda}$.

We use parameter estimation samples released LVK: Overall_posterior parameter
estimation samples for BBHs detected in GWTC-1; PrecessingSpinIMRHM samples
for events published in GWTC-2 and GWTC-2.1 (Abbott et al., 2021b); and
C01:Mixed for events in GWTC-3 (Abbott et al., 2021a).222 Parameter estimation
samples are available at
GWTC-1: https://dcc.ligo.org/LIGO-P1800370/public
GWTC-2: https://dcc.ligo.org/LIGO-P2000223/public
GWTC-2.1: https://zenodo.org/record/5117703
GWTC-3: https://zenodo.org/record/5546663 .

We evaluate the detection efficiency $\xi(\vec{\Lambda})$ using successfully
recovered BBH injections, released by LVK in
o1+o2+o3_bbhpop_real+semianalytic-
LIGO-T2100377-v2.hdf5333https://zenodo.org/record/5636816.

Figure 5:  Posteriors of all the hyperparameters in the s-mix model where
($\mu_{\chi,I}$, $\mu_{\chi,A}$, $\sigma_{\chi,I}^{2}$, $\sigma_{\chi,A}^{2}$
fit the $\vec{\theta}_{\rm mix}=\\{\chi_{1},\chi_{2}\\}$, while $\zeta$ and
$\sigma_{t}$ fit mixing fraction and tilts distribution of aligned systems. We
show the $68\%$ and $90\%$ intervals. The grey region marks the parameter
space in $(\mu_{\chi,I},\sigma_{\chi,I}^{2})$ and
$(\mu_{\chi,A},\sigma_{\chi,A}^{2})$ that is excluded to avoid the singularity
of the Beta distribution [Eq. (B5)]. The grey dotted lines are the resulting
effective priors $p_{\rm eff}(\mu_{\chi})$ and $p_{\rm
eff}(\sigma_{\chi}^{2})$ [Eq. (B6) and Eq. (B7) respectively].

##  Appendix B Spin and Redshift distributions

In the previous sections, we analyzed how the pairing function differs between
the isotropic and aligned systems. This section will focus on the differences
in their spin magnitude and redshift distributions.

Spin magnitudes could be similar for both field binaries and cluster binaries
if they are only controlled by stellar collapse physics (O’Connor & Ott,
2011). However, during the isolated evolution of binaries, processes like
stable mass transfer and tidal interactions can spin up the BHs (Bavera et
al., 2021; du Buisson et al., 2020). In dense stellar clusters, BHs born from
previous mergers (instead of stellar collapse) have a characteristic spin
$\sim 0.7$ (Gerosa & Berti, 2017; Fishbach et al., 2017; Kovetz et al., 2018).
If these BHs merge repeatedly, a small fraction of BBHs originating from
clusters could have larger spins. These differences in spin magnitudes of
field and cluster binaries, could be observed in the population of detected
BBHs. To analyse if this is the case, we will allow the aligned and isotropic
binaries to have different spin-magnitude distributions i.e.
$\vec{\theta}_{\rm mix}=\\{\chi_{1},\chi_{2}\\}$ governed by hyperparameters
$\vec{\Lambda}_{\rm
mix}=\\{\mu_{\chi,I},\sigma^{2}_{\chi,I},\mu_{\chi,A},\sigma^{2}_{\chi,A}\\}$.
We call this model s-mix. Here, we assume component spin magnitudes
$\chi_{1,2}$ are identically distributed according to a Beta distribution
(Talbot & Thrane, 2017) (we have dropped subscripts “I” and “A” for brevity),

$p(\chi_{i}|\alpha,\beta)=\frac{\chi_{i}^{1-\alpha}\,(1-\chi_{i})^{1-\beta}}{c(\alpha,\beta)}\,,$
(B1)

where
$c(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}$ is
a normalization constant, where $\Gamma$ is the gamma function. We sample the
the mean and variance ($\mu_{\chi},\sigma^{2}_{\chi}$) of the beta
distribution (and not in $(\alpha,\,\beta)$ directly)

$\displaystyle\mu_{\chi}$
$\displaystyle=\frac{\alpha_{\chi}}{\alpha_{\chi}+\beta_{\chi}}$ (B2)
$\displaystyle\sigma_{\chi}^{2}$
$\displaystyle=\frac{\alpha_{\chi}\beta_{\chi}}{(\alpha_{\chi}+\beta_{\chi})^{2}(\alpha_{\chi}+\beta_{\chi}+1)}.$
(B3)

The shape parameters $\alpha$ and $\beta$ can be recovered from $\mu_{\chi}$
and $\sigma^{2}_{\chi}$ through

$\displaystyle\alpha_{\chi}=$
$\displaystyle\frac{\mu_{\chi}^{2}(1-\mu_{\chi})}{\sigma_{\chi}^{2}}-\mu_{\chi},$
$\displaystyle\beta_{\chi}=$
$\displaystyle\frac{\mu_{\chi}(1-\mu_{\chi})^{2}}{\sigma_{\chi}^{2}}-(1-\mu_{\chi})$
(B4)

It also implies that there are certain constraints on the distribution’s mean
and variance $(\mu,\variation)$ (which is what we actually sample). To ensure
$\alpha_{\chi}>1,\,\beta_{\chi}>1$

$\sigma_{\chi}^{2}<\begin{dcases}\frac{\mu_{\chi}^{2}(a_{\rm
max}-\mu_{\chi})}{(a_{\rm max}+\mu_{\chi})}&0<\mu_{\chi}<\frac{a_{\rm
max}}{2}\\\ \frac{\mu_{\chi}(a_{\rm max}-\mu_{\chi})^{2}}{(2a_{\rm
max}-\mu_{\chi})}&\frac{a_{\rm max}}{2}<\mu_{\chi}<a_{\rm max}\end{dcases}$
(B5)

and the maximum variance that can be achieved is $a_{\rm max}^{2}/12$ at
$\mu=a_{\rm max}/2$. We will assume $a_{\rm max}=1$ and that the priors on
$\mu_{\chi,I}$ and $\mu_{\chi,A}$ are uniform between $0$ and $1$, while the
prior on $\sigma^{2}_{\chi,I}$ and $\sigma^{2}_{\chi,A}$ are uniform between
$0$ and $1/12$. However, since we have imposed that the Beta distribution is
not singular, a vast region in the prior space is restricted. In this case,
the effective one-dimensional prior on $\mu_{\chi}$ is

$p_{\rm
eff}(\mu_{\chi})=\frac{1}{c}\begin{dcases}\frac{\mu_{\chi}^{2}(1-\mu_{\chi})}{(1+\mu_{\chi})}&0<\mu_{\chi}<\frac{1}{2}\\\
\frac{\mu_{\chi}(1-\mu_{\chi})^{2}}{(2-\mu_{\chi})}&\frac{1}{2}<\mu_{\chi}<1\end{dcases}$
(B6)

where $c=4\log({3}/{2})-{19}/{12}$ is the normalization constant, while the
effective one dimensional prior on $\sigma_{\chi}^{2}$ is

$p_{\rm eff}(\sigma_{\chi}^{2})\propto
1+2\real\left((1+i\sqrt{3})\mathcal{S}_{\chi}^{1/3}\right)$ (B7)

where

$\mathcal{S}_{\chi}=1+3\sigma_{\chi}\left(\sqrt{3(\sigma_{\chi}^{4}+11\sigma_{\chi}^{2}-1)}-6\sigma_{\chi}\right)\,.$

Fig. 5 shows results when analyzing GWTC-3 binary BHs with the s-mix. We plot
the $68\%$ and $90\%$ intervals for the mean and variance of $\chi_{1,2}$
distribution for isotropic and aligned populations
($\mu_{\chi,I},\sigma^{2}_{\chi,I}$, $\mu_{\chi,A},\sigma^{2}_{\chi,A}$), the
fraction of aligned binaries ($\zeta$) and standard deviation of aligned-spin
tilts $\sigma_{t}$. We also show the effective priors $p_{\rm
eff}(\mu_{\chi})$ and $p_{\rm eff}(\sigma_{\chi}^{2})$ [Eq. (B6) and Eq. (B7)
respectively] as well as the region excluded to avoid the singularity of the
Beta distribution [Eq. (B5)]. We see the presence of multiple modes in Fig. 5.
However, these results are consistent with the distribution of spin magnitudes
obtained by Default spin model in Abbott et al. (2021d). We observe two modes
in the mixing fraction, $\zeta\simeq 0.1$ (implying that majority of systems
are isotropic) and $\zeta\simeq 0.9$ (implying the majority of systems are
aligned). When aligned systems are favored, we find that their spin-magnitude
distribution peaks at $\mu_{\chi,A}\simeq 0.24$. This is consistent with the
Default spin model in Abbott et al. (2021d). Also, the spin magnitudes of
isotropic binaries are consistent with the prior when $\zeta$ is large. On the
contrary, when isotropic systems are favored (small $\zeta$), isotropic
binaries reproduce the Default spin model while the aligned binaries extract
the prior. In both cases, we find that the overall model is effectively trying
to reproduce the Default model.

For simplicity, in the above model (s-mix), we only fit $\vec{\Lambda}_{\rm
mix}=\\{\mu_{\chi,I},\sigma^{2}_{\chi,I},\mu_{\chi,A},\sigma^{2}_{\chi,A}\\}$
and $\zeta,\sigma_{t}$. We keep all the other hyperparameters associated with
describing the mass and redshift distributions fixed to their median values
obtained in Abbott et al. (2021d). However, Abbott et al. (2021d); Callister
et al. (2021) have shown that the mass ratios and spins of detected BBHs show
evidence of anti-correlation. Since $\chi_{\rm eff}$ depends on the spin
tilts, this would require that distribution of spins be fitted alongside the
mass-ratio distribution. We use two approaches to test if fitting $q$
distribution alongside spins affects our results. In the first approach, we
assume isotropic and aligned binaries have different spin-magnitude
distributions (governed by $\vec{\Lambda}_{\rm
mix}=\\{\mu_{\chi,I},\sigma^{2}_{\chi,I},\mu_{\chi,A},\sigma^{2}_{\chi,A}\\}$)
but the same mass ratio distribution (governed by $\beta_{q}$). This model
yields a result similar to the s-mix model with $\beta_{q}$ consistent with
Powerlaw+Peak model in Abbott et al. (2021d). In the second approach, we
assume isotropic and aligned binaries have both different spin-magnitude
distributions (governed by
$\\{\mu_{\chi,I},\sigma^{2}_{\chi,I},\mu_{\chi,A},\sigma^{2}_{\chi,A}\\}$) and
different mass ratio distribution (governed by
$\\{\beta_{q,I},\beta_{q,A}\\}$). In this case, we find that the spin
magnitudes and tilts reproduce the s-mix model with bimodal features. However,
$\\{\beta_{q,I},\beta_{q,A}\\}$ do not follow the anti-correlation observed in
the q-mix model. Instead, $\\{\beta_{q,I},\beta_{q,A}\\}$ also show bimodality
similar to spin magnitudes i.e. dominant channel recovers the LVK result while
the subdominant channel is consistent with prior. This is because assuming
different spin-magnitude distributions for aligned and isotropic binaries
enforces that either channel dominates and reproduces the LVK result. In Sec.
3, we discussed that features in $\beta_{q,I}$ and $\beta_{q,A}$, observed
with q-mix model, could help explain the $q-\chi_{\rm eff}$ correlation
observed in GWTC-3 (Abbott et al., 2021d). If that were the case, one would
also expect that isotropic and aligned subpopulations have different spin
magnitudes. This could be understood through following toy model. Let’s assume
that BBHs with aligned (isotropic) binaries have spin magnitudes $\chi_{A}$
($\chi_{I}$). Also the q-mix model is consistent with aligned (isotropic)
binaries having $q\simeq 1$ (small mass ratios). In that case, the $\chi_{\rm
eff}$ distribution of isotropic binaries should be

$\chi_{\rm eff,I}=\chi_{I}\frac{\cos\theta_{1}+\cos\theta_{2}}{2}$ (B8)

This yields a triangle distribution (see Baibhav et al. (2020) for derivation)
of $\chi_{\rm eff,I}$ between $-\chi_{I}$ and $\chi_{I}$. On the other hand,
$\chi_{\rm eff}$ distribution of aligned binaries (assuming small mass ratios)
should be

$\chi_{\rm eff,A}\simeq\chi_{A}\cos\theta_{1}$ (B9)

This would yield $\chi_{\rm eff}$ distribution symmetric around $0$ at
$q\simeq 1$, and positive values with peak at $\chi_{\rm eff}=\chi_{A}$ at
small $q$. However, since the current data can not discern differences in the
spin-magnitude distribution of isotropic and aligned binaries, we can not test
our hypothesis that $q-\chi_{\rm eff}$ correlation is caused by isotropic and
aligned binaries having different properties. However, as more information is
available with the next observing runs of LVK, it might be possible to gain
insight into the cause behind such features.

The formation history of BBHs in most channels is predicted to be set by the
star formation history and hence, could be similar across binaries originating
in galactic fields and stellar clusters. However, different channels have
different processes that drive the merger. These processes are not equally
efficient and might have different time delays between binary formation and
merger (van Son et al., 2022; Mapelli & Giacobbo, 2018; Baibhav et al., 2019).
For example, even in the field, binaries mergers driven by stable mass
transfer take significantly longer than those driven by CE. These differences
might be observable in the redshift distribution of BBHs detected by LVK
(Fishbach & Kalogera, 2021). To analyze these differences, we allow the
aligned and isotropic binaries to have redshift distributions, i.e.,
$\vec{\theta}_{\rm mix}=\\{z\\}$ governed by hyperparameters
$\vec{\Lambda}_{\rm mix}=\\{\kappa_{I},\kappa_{A}\\}$. We call this model
z-mix. Here $\kappa_{I,A}$ is the power-law slope that governs the evolution
of the source-frame merger rate,

$p(z|\kappa_{I,A})\propto\frac{1}{1+z}\frac{dV_{c}}{dz}\left(1+z\right)^{\kappa_{I,A}}\,,$
(B10)

where $\frac{dV_{c}}{dz}$ is the differential comoving volume per unit
redshift. We assume an uniform prior on $\kappa_{I}$ and $\kappa_{A}$ between
$-10$ and $10$. We find that z-mix model can not distinguish between the
redshift distribution of isotropic and aligned systems. We find that the
$\zeta$ distribution is consistent with flat prior (unlike s-mix where $\zeta$
was bimodal). In addition, we find that either $\kappa_{I}$ or $\kappa_{A}$
recover $\kappa$ calculated in Abbott et al. (2021d).