# Gravitational waves carry information beyond effective spin parameters but
it is hard to extract

Simona J. Miller<EMAIL_ADDRESS>Department of Physics, California
Institute of Technology, Pasadena, California 91125, USA  LIGO Laboratory,
California Institute of Technology, Pasadena, California 91125, USA  Zoe Ko
<EMAIL_ADDRESS>Department of Physics, University of California Berkeley,
Berkeley, California 94720, USA  Tom Callister<EMAIL_ADDRESS>Kavli
Institute for Cosmological Physics, The University of Chicago, Chicago,
Illinois 60637, USA  Katerina Chatziioannou<EMAIL_ADDRESS>Department of Physics, California Institute of Technology, Pasadena,
California 91125, USA  LIGO Laboratory, California Institute of Technology,
Pasadena, California 91125, USA

###### Abstract

Gravitational wave observations of binary black hole mergers probe their
astrophysical origins via the binary spin, namely the spin magnitudes and
directions of each component black hole, together described by six degrees of
freedom. However, the emitted signals primarily depend on two effective spin
parameters that condense the spin degrees of freedom to those parallel and
those perpendicular to the orbital plane. Given this reduction in
dimensionality between the physically relevant problem and what is typically
measurable, we revisit the question of whether information about the component
spin magnitudes and directions can successfully be recovered via
gravitational-wave observations, or if we simply extrapolate information about
the distributions of effective spin parameters. To this end, we simulate three
astrophysical populations with the same underlying effective-spin distribution
but different spin magnitude and tilt distributions, on which we conduct full
individual-event and population-level parameter estimation. We find that
parameterized population models can indeed qualitatively distinguish between
populations with different spin magnitude and tilt distributions at current
sensitivity. However, it remains challenging to either accurately recover the
true distribution or to diagnose biases due to model misspecification. We
attribute the former to practical challenges of dealing with high-dimensional
posterior distributions, and the latter to the fact that each individual event
carries very little information about the full six spin degrees of freedom.

##  I Introduction

The spins of black holes (BHs) in binaries (BBHs) are a unique probe of
physics on multiple scales, from fundamental BH properties to stellar
interiors and the astrophysical environments in which compact binaries form.
Each binary possesses six spin degrees of freedom: the spin magnitudes, polar
angles (tilts) and azimuthal angles of each binary component [1]. BH spins are
encoded in the gravitational waves (GWs) the binary emits and can, at least in
principle, be constrained from observation [2, 3] by the LIGO [4] and Virgo
[5] detectors. The magnitudes and directions of the spins at merger are
determined by the spin each BH has upon formation as well as the binary’s
evolutionary history, e.g. [6, 7, 8]. Spin measurements are therefore a
promising way to determine whether BBHs form dynamically or in the field, e.g.
[9, 10, 11, 12, 13, 14], and answer questions such as the role of angular
momentum transfer in stars, tidal interactions, and mass transfer, e.g. [15,
16, 17, 18, 19].

Despite their astrophysical importance, spins remain poorly constrained in GW
data. Their imprint on the signal is typically subdominant to other intrinsic
effects such as the BH masses, e.g. [20, 21, 22, 23, 24]. Furthermore, not all
six spin degrees of freedom affect the signal equally. Though waveform models
formally depend on the full spin vectors [25, 26, 27, 28], analytical post-
Newtonian calculations indicate that the dominant spin effect is captured by
two effective parameters: the effective aligned spin $\chi_{\mathrm{eff}}$
that includes the spin components parallel to the Newtonian orbital angular
momentum [29], and the effective precessing parameter $\chi_{\mathrm{p}}$ that
includes the perpendicular components [30]. The former primarily affects the
length of the signal while the latter describes spin-precession, the change in
binary orientation due to spin-orbit and spin-spin interactions [31].
Unsurprisingly, then, constraints on the astrophysical distributions of
$\chi_{\mathrm{eff}}$ and $\chi_{\mathrm{p}}$ can be typically obtained with
fewer observations and are less prone to population model systematics than the
spin components [32, 33, 34].

Although less well measurable, it is instead the underlying spin components
that are of prime astrophysical interest. GW signals contain some information
about component spins. However, unlike $\chi_{\mathrm{eff}}$ and
$\chi_{\mathrm{p}}$ that appear prominently in the GW phase and amplitude and
whose measurability can be predicted with analytic arguments [35, 36],
individual spin components have a significantly subdominant effect on the
waveform. The resulting constraints on the astrophysical distribution of spin
components are correspondingly weaker and in many cases subject to
uncertainties about the role of population models [34, 37]. Indeed, even
though it is widely accepted that BBHs have a range of $\chi_{\mathrm{eff}}$
values that are not symmetric about zero [38, 2, 3] and that not all BBHs have
a vanishing $\chi_{\mathrm{p}}$ [2, 3], the exact shape of the inferred
distribution for spin magnitudes and directions depends on the
parameterization of the corresponding population model. For example, different
parametrizations for the angle between the spins and the Newtonian orbital
angular momentum lead to varied conclusions about where the distribution peaks
and the degree of spin-orbit misalignment [3, 34, 39, 37, 40, 41].

Central to this discussion are the questions of how much information GW
signals actually contain about the BH spin components versus
$\chi_{\mathrm{eff}}$ alone, how feasible it practically is to reliably
extract this information, and the extent to which conclusions are driven by
informative data or simply by overly restrictive models. In this paper, we
approach these issues by posing three questions, from which we conclude:

  1. 1.

Do GWs carry information about spin components, or are we just extrapolating
the effective aligned spin $\chi_{\mathrm{eff}}$? (Sec. IV)

_Yes, we can distinguish between populations with low, moderate, and high
spins even when they have identical effective spin distributions._

  2. 2.

Can component spin distributions be accurately measured? (Sec. V)

_Even though we can qualitatively tell apart BH populations with different
spin distributions, characterizing them accurately is practically
challenging._

  3. 3.

Can we tell when measurements of component spin distributions are biased?
(Sec. VI)

_Common tests based on posterior predictive checks cannot identify modeling
biases in component spin distributions due to the fact that individual-event
posteriors are extremely weakly informative about spin components._

The remainder of this paper presents our analysis in support of these
conclusions. We discuss spin degrees of freedom, effective spin parameters,
and the notation used throughout in Sec. II. Our Methods are briefly described
in Sec. III, and are expanded upon in the Appendices. Results about measuring
the component spin distributions are presented in Sec. IV. Section V
introduces the extensive series of verification methods – both population and
individual-event level – we use to ensure the robustness of our results, all
of which are further elaborated upon in the Appendices. In Sec. VI, we
identify limitations of the traditional method of using posterior predictive
checks to assess biased population measurements, and identify the sources of
this bias. We compare our findings to those of past work in Sec. VII, and then
conclude in Sec. VIII.

##  II Spin magnitudes and tilts versus effective spin parameters

Each BH in the binary is described by a dimensionless spin vector
$\vec{\chi}_{i}$, $i\in\\{1,2\\}$. In a coordinate system where the $z$-axis
is aligned with the binary Newtonian orbital angular momentum $\vec{L}$, the
spin vector is characterized by a magnitude $\chi_{i}\in[0,1]$, polar angle
$\theta_{i}\in[0,\pi]$, and azimuthal angle $\phi_{i}\in[0,2\pi]$. Modulo
horizon absorption effects [42], the spin magnitude is constant throughout the
binary evolution, while the spin angles evolve due to spin-orbit and spin-spin
interactions causing the spin vector to precess [31, 43].

The full six spin degrees of freedom remain relatively poorly constrained by
GW signals. Rather, the dominant spin effects are expressed by two effective
parameters. The mass-weighted average spin projected onto $\vec{L}$

$\chi_{\mathrm{eff}}=\frac{\chi_{1}\cos\theta_{1}+q\chi_{2}\cos\theta_{2}}{1+q}\in(-1,1)\,,$
(1)

is referred to as _effective aligned spin_ , where we have defined the binary
mass ratio $q\equiv m_{2}/m_{1}$, where $m_{1}\geq m_{2}$ are the BH masses.
The effective aligned spin is, in general, better constrained as it is related
(in the equal mass limit) to the leading-order spin contribution in the post-
Newtonian expansion for the GW inspiral phase. Additionally,
$\chi_{\mathrm{eff}}$ is conserved under spin-precession and radiation
reaction to at least the second post-Newtonian order [29].

Spin-precession effects are captured with the _effective precessing parameter_

$\chi_{\mathrm{p}}=\mathrm{max}\left[\chi_{1}\sin\theta_{1},\left(\frac{3+4q}{4+3q}\right)q\,\chi_{2}\sin\theta_{2}\right]\in[0,1)\,\,.$
(2)

This parameter and its extensions [44, 45] are motivated by the fact that
spin-orbit precession (and the GW amplitude and phase modulations it induces)
are driven by in-plane spin components [46, 30]. Constraints on
$\chi_{\mathrm{p}}$ are typically much weaker than $\chi_{\mathrm{eff}}$
especially given the observed absence of large spin-precession in BBHs [2, 3].
In what follows, we therefore focus on $\chi_{\mathrm{eff}}$.

##  III Methodology

Figure 1: Spin distributions for the three simulated BBH populations we use to
assess the amount of recoverable information GW signals contain about
component spin distributions. The three populations share the same
$\chi_{\mathrm{eff}}$ distribution (chosen for consistency with current data),
but differ in their underlying component spin and $\chi_{\mathrm{p}}$
distributions. From left to right, panels show the spin magnitude $\chi$, spin
tilt $\cos\theta$, the effective spin parameter $\chi_{\mathrm{eff}}$, and the
effective precession parameter $\chi_{\mathrm{p}}$. (Navy) HighSpinPrecessing:
a BBH population with few vanishing spins that are preferentially oriented
close to the orbital plane. (Pink) MediumSpin: a BBH population with moderate
spin magnitudes peaking at $\chi=0.25$ and preferentially aligned with the
Newtonian orbital angular momentum. (Orange) LowSpinAligned: a BBH population
with low spin magnitudes peaking at $\chi=0.10$ with both strongly aligned and
anti-aligned sub-populations.

In order to isolate the amount of information included in GW signals about the
spin components relative to the effective spin, we simulate astrophysical
populations with identical $\chi_{\mathrm{eff}}$ distributions but different
underlying component spin distributions. We choose a $\chi_{\mathrm{eff}}$
distribution that is qualitatively similar to current constraints [3, 34, 47]
and decompose it into three populations with distinct spin magnitudes and tilt
angle distributions. The latter are not astrophysically motivated, but rather
selected as distinct test cases of potential distributions.

The three simulated astrophysical distributions are shown in Fig. 1, with
further details given in App. A.:

  * •

The HighSpinPrecessing population contains BHs with the most extremal spins
and tilts: the majority of the population has $\chi>0.5$ and tilts nearly in-
plane, corresponding to significant spin-precession.

  * •

The MediumSpin population is most similar to current constraints:
preferentially small to moderate spin magnitudes peaking at $0.25$, and a wide
range of tilts with a preference for alignment compared to anti-alignment.

  * •

The LowSpinAligned population has the smallest spin magnitudes, with nearly
all BHs having $\chi<0.5$. Uniquely, this population has a bimodal spin tilt
angle distribution, with a larger peak at $\cos\theta=1$ (perfect alignment)
and a smaller peak at $\cos\theta=-1$ (perfect anti-alignment). It is
therefore a test case of sensitivity to mixture models.

With these three populations, we conduct a full end-to-end injection/recovery
campaign. We draw parameters describing individual GW events from each
distribution, restrict to detectable events with a network optimal signal-to-
noise ratio (SNR) above $10$ in the LIGO Livingston, LIGO Hanford, and Virgo
detectors, simulate data assuming O3 sensitivity [48], and obtain samples from
the multidimensional posterior distribution of the binary parameters for each
event individually. We then hierarchically model the population distribution
of the simulated posteriors with parametrized population models.

The individual-event posterior sampling is conducted with the nested sampler
Dynesty [49] as implemented in Bilby [50, 51]. We use the IMRPhenomXPHM
waveform model [27] both for simulation and recovery as it models all six spin
degrees of freedom, contains higher order radiation modes, and is the least
computationally expensive option available. Although more computationally
expensive than approximate parameter estimation [52, 53, 54], it is essential
to use full stochastic sampling for this work. As we are trying to discern
subtle effects in the signals, we must properly characterize the individual-
event likelihoods. Full details about parameter estimation settings are given
in App. B. For hierarchical inference, we primarily use the Markov Chain Monte
Carlo sampler Emcee [55], with some follow-up studies run with Numpyro [56,
57]. The full hierarchical inference procedure is outlined in App. C, with the
parameterized population models detailed in App. D.

##  IV Different spin magnitude and tilt distributions can be distinguished

Figure 2:  Inferred distributions for various spin parameters (from left to
right: spin magnitude $\chi$, spin tilt $\cos\theta$, effective
$\chi_{\mathrm{eff}}$ spin, effective precessing spin $\chi_{\mathrm{p}}$) and
the three simulated populations (top to bottom). All results are obtained with
the Beta+DoubleGaussian model using 70 (pink) and 300 (navy) events. Traces
correspond to draws from the population posterior and solid lines enclose 90%
of the probability. The black solid line corresponds to the true, underlying
population. The dashed lines show the 90% credible intervals inferred by
sampling the prior on the population-level parameters, including the effective
sample cut as defined in Eq. (9).

Using the results of the signal injection and parameter estimation campaign,
we perform hierarchical inference [58, 59, 60] on events drawn from each
population shown in Fig. 1 in order to reconstruct their underlying spin
distributions. Population inference requires the adoption of a model for the
component spin and tilt distributions. We select an analytic model in which
spin magnitudes follow a non-singular Beta distribution and the spin tilt
angles follow a bimodal Gaussian distribution, which we hereby refer to as the
Beta+DoubleGaussian; see App. D.2 for a full description. While the true,
underlying distributions shown in Fig. 1 were not explicitly drawn from the
Beta+DoubleGaussian model, this model is expected to agree with each simulated
population to within statistical uncertainties.111We confirm that the
Beta+DoubleGaussian model is a good fit to the underlying populations through
a least-squares fit. The Kullback-Leibler divergences [61] between the best
fit and true underlying distributions are $<10^{-4}$ for all spin magnitude
distributions and $<0.08$ for all cosine tilt distributions. As a further
check, we also run hierarchical inference on catalogs of simulated, Gaussian
individual-event spin posteriors, with results discussed in App. G.2. We do
not model the spin azimuthal angle, effectively (and correctly) assuming that
it is distributed according to its uniform prior.

Fig. 2 shows the inferred distributions for various spin parameters for the
three simulated populations under the Beta+DoubleGaussian population model.
The results from a 70-event catalog are plotted in pink, and from a 300-event
catalog in navy, chosen to mimic O3 and projected O4 catalog size
respectively. Black, bold traces show the true underlying populations for
comparison. The $\chi$ (first column) and $\cos\theta$ (second column)
distributions are generated by random draws from the posteriors on the
population parameters, shown in Figs. 6, 7, and 8 in App. E.

The $\chi_{\mathrm{eff}}$ (third column) and $\chi_{\mathrm{p}}$ (fourth
column) distributions are generated by (i) randomly drawing from the inferred
$\chi$, $\cos\theta$, and $q$ distributions, (ii) calculating the effective
spin parameters from these draws, and (iii) generating a Gaussian kernel
density estimate. We find that the $\chi_{\mathrm{eff}}$ distributions are
reconstructed accurately across all three populations and for both catalog
sizes.

Switching to the component spins, we can qualitatively distinguish between
their distributions for each population. The spin magnitude inferred for the
HighSpinPrecessing population is the widest and has the largest mean. From
HighSpinPrecessing to MediumSpin to LowSpinAligned, the means and widths of
the inferred distributions get progressively smaller, as is the case for the
true, underlying populations. The mass ratio distributions are also
successfully recovered for all three populations and two catalog sizes, as
shown in App. E. We, therefore, _can distinguish between populations with low,
moderate, and high spins when they have identical $\chi_{\mathrm{eff}}$
distributions_.

Although we can qualitatively characterize the spin magnitude and tilt
distributions among these three populations, in some cases we cannot reliably
characterize their properties accurately. Specifically, the true underlying
distribution the LowSpinAligned population does not lie within the 90%
credible reconstructed region as measured with the Beta+DoubleGaussian model.
The bimodality of this population’s tilt angle distribution is not recovered
and the inferred spin magnitude distribution has a higher mean than truth. We
confirm that this mismatch between the true and inferred distributions is not
solely driven by overly restrictive priors on the population-level parameters,
plotted with dashed lines in Fig. 2. For example, the bias in the inferred
spin magnitude distribution for the HighSpinPrecessing population occurs over
a region where the prior generates a flat distribution.

Population measurements of the effective spin are more robust against bias
than component spins. Even a notable mismatch between the true and recovered
spin magnitude and tilt distributions for the LowSpinAligned population
results in a precisely and accurately constrained $\chi_{\mathrm{eff}}$
distribution. The $\chi_{\mathrm{p}}$ distributions are more susceptible to
inaccurate recovery in the component spins, effectively inheriting their
biases. For example, for the $\chi_{\mathrm{p}}$ distribution of the
LowSpinAligned population: the inference of more in-plane spin
($\cos\theta\sim 0$) combined with the slight over-estimation of the mean of
the spin magnitude distribution, leads to a corresponding over-estimate of the
bulk of the $\chi_{\mathrm{p}}$ distribution. This means that
$\chi_{\mathrm{eff}}$, but not necessarily $\chi_{\mathrm{p}}$ can be reliably
characterized on a population level by component spin measurements.

##  V Difficulties of Measuring Component Spin Distributions

The biased reconstruction of the LowSpinAligned population is unexpected. The
injection and recovery campaign was performed using the same waveform model,
and selection effects were self-consistently handled in both signal selection
and parameter estimation. Under these conditions, there is no a priori reason
why population recovery should fail. As such, our results instead suggest a
shortcoming in either the parameter estimation or population recovery stages
of the analysis.

To diagnose this shortcoming, we employ a slew of checks, all of which are
further elaborated upon in App. G. To ensure that the problem is not our
hierarchical inference framework and implementation we do the following:

  * •

Simulate Gaussian individual-event spin posteriors (App. G.2) – For each of
the 300 events per population, we generate a series of simulated Gaussian
individual-event spin posteriors with a range of measurement errors and
underlying correlations, and use these as input to hierarchical inference. In
these cases, the Beta+DoubleGaussian population model is able to recover the
underlying populations, as seen in Fig. 10.

Implications: The hierarchical inference and selection effect framework is
algorithmically robust, and the Beta+DoubleGaussian model is able to recover
the true population distributions to within statistical uncertainties.

  * •

Fix either the spin magnitude or tilt angle distribution to the truth (App.
G.3) – When only fitting for the $\chi$ or $\cos\theta$ population and not the
other, we are still unable to recover the correct distribution for the
LowSpinAligned population; see Fig. 12.

Implications: The observed bias in the spin magnitude and tilt angle
distributions is not related to correlations between the two distributions.

  * •

Use a different sampler for the hierarchical likelihood (App. G.4) – We repeat
the analysis of Fig. 2 with an independently-implemented hierarchical
inference code that is based on Numpyro instead of emcee. We obtain
essentially identical results, shown in Fig. 13.

Implications: The hierarchical inference and selection effect framework is
algorithmically robust.

  * •

Fit for the mass and redshift distributions instead of fixing it to truth
(App. G.5) – Our main results fit for the spin magnitude, spin tilt, and mass
ratio distributions, while fixing the distributions of the primary mass and
redshift to their true population values, given in App. A. Figure 13 extends
these results to also fit for the mass and redshift distributions and shows
the corresponding spin population posteriors, which remain unchanged. The mass
and redshift distributions are recovered with no bias.

Implications: We have not misspecified the mass or redshift distributions when
fixing them to truth during hierarchical inference, nor biased results of our
spin inference by neglecting to simultaneously fit for the mass and redshift
distributions.

  * •

Plot rates instead of probability distributions (App. G.6) – Apparent
disagreement between injected and recovered probability distributions can
sometimes be caused by neglecting to simultaneously measure and plot the
absolute rate of mergers. In the main text we do not infer the overall rates
of black hole mergers, but only the shapes of their spin distributions. To
check if neglecting the merger rate contributes to apparent disagreement
between injected and recovered populations, we repeat our hierarchical
inference while also fitting for the rate of black hole mergers as function of
spin magnitude and tilt. This yields the results shown in Fig. 14, which
remain qualitatively similar to those in Fig. 2.

Implications: We have not biased results of our spin measurements by failing
to measure or plot absolute merger rates, rather than probability
distributions.

  * •

Exclude spin selection effects (App. G.7) – The selection function only
negligibly affects spin magnitudes and tilt angles. To ensure that we are not
incorrectly implementing the selection function in the hierarchical
likelihood, we conduct hierarchical inference without including selection
effects in spin. This does not impact our results, as can be seen in Fig. 15.

Implications: The implementation of the selection function is algorithmically
robust.

  * •

Employ different methods of breaking the degeneracy in the bimodal Gaussian
model (App. G.7) – For a bimodal distribution, some method must be imposed to
break the degeneracy between the two components of the model. For the
Beta+DoubleGaussian model, this can be done in one of three ways: imposing an
ordering of the means, the widths, or limiting the mixing fraction be $\leq
0.5$. Sometimes one method of breaking the degeneracy converges better than
another. As shown in Fig. 15, we find that this is not the case here and
different methods perform comparably.

Implications: Our choice of degeneracy-breaking between the two Gaussian
components in our population model is not causing convergence issues.

  * •

Run hierarchical inference on different 70-event catalog instantiations (App.
G.7) – Finally, to get a sense of how much the specific 70 events we select
from the underlying population affect hierarchical inference, we repeat the
procedure with several different catalog instantiations. While there is
expected variance in the results –see Fig. 16– it cannot account for the
degree of mismatch seen in the bottom row of Fig. 2. Additionally, each
catalog instantiation leads to a different number of per-event effective
samples, which we find are not correlated to the goodness of fit.

Implications: The observed bias does not arise from an insufficient number of
per-event effective samples.

We then move towards investigating the underlying individual-event parameter
estimation with the following checks:

  * •

Sampler settings – We run parameter estimation with a large variety of sampler
settings in Bilby, and eventually adopt the standard, reviewed settings for
our headline results of Fig. 2.

Implications: Running with more aggressive sampler settings in Bilby may fix
convergence problems, but this was not the case for any configurations we
employed.

  * •

Probability-probability plots (App. G.1) – We generate probability-probability
(P-P) plots [62, 63] for reweighted individual-event Bilby posteriors. As seen
in Fig. 9, the test passes.

Implications: Either the Bilby individual-event posterior samples are
unbiased, or the biases are subtle enough to not be detectable by a P-P test,
as warned against in [64].

  * •

Use a different waveform model (App. G.3) – We re-run individual-event
inference on the same sets of events with Bilby using the IMRPhenomXP waveform
model instead of IMRPhenomXPHM both for injection and recovery. Results with
this waveform model are comparable or worse to that presented in the main text
with IMRPhenomXPHM, although the bimodality of the LowSpinAligned population
is slightly better constrained; see Fig. 11.

Implications: The existence of bias in the measured spin magnitude and tilt
distributions is not driven by our choice of waveform model, although the
specific details of how that bias manifests appear to be, i.e. different
waveforms yield different population-level results. This indicates that the
bias may be due to individual-event sampling issues.

  * •

Fix non-spin parameters to truth in individual-event sampling (App. G.3) –
Finally, we conduct individual-event inference with IMRPhenomXPHM fixing all
parameters aside from the spin magnitudes to tilt angles to truth (i.e. use
delta function priors at their injected values). In this case, the
Beta+DoubleGaussian population model is able to successfully recover the truth
for all three populations; see Fig. 11.

Implications: The added complexity going from sampling just spins to all
fifteen binary parameters is a likely culprit for the biased spin magnitude
and tilt angle distributions.

Although we can qualitatively tell apart the different populations in Fig. 1,
our results indicate that the spin distribution of all possible BBH
populations cannot necessarily be accurately measured. Full parameter
estimation with spin-precession is a technically challenging analysis. Despite
conducting tens of model checking procedures, we cannot fully identify the
driving source of the bias observed in Fig. 2. We hypothesize that the error
is due to issues related to sampling from the high-dimensional posterior for
individual events, as suggested by the final bulletpoint above.

##  VI Identifying bias is difficult: limitations of posterior predictive
checks

Figure 3:  Posterior predictive checks (PCCs) for spin parameters (left to
right: component spin magnitudes $\chi$, component spin tilts $\cos\theta$,
effective spin $\chi_{\mathrm{eff}}$, and effective precessing spin
$\chi_{\mathrm{p}}$) of the LowSpinAligned population under the
Beta+DoubleGaussian population model. Each trace is one catalog from the
observed and the predicted populations. 100 catalogs are shown, for results
with 70 (pink) and 300 (navy) events. In the absence of discrepancy between
the inferred and true population distributions, the traces should on average
follow the diagonal. However, here there is such a discrepancy, as can be seen
in the bottom row of Fig. 2, and the traces do on average follow the diagonal,
meaning that this PPC is not necessarily a good diagnostic tool for component
spins.

While the Beta+DoubleGaussian model was able to produce qualitatively correct
results for each of the three distinct populations, it was not able to
successfully recover the true underlying populations. This leads us to the
question: can commonly-used modeling diagnostics successfully identify poorly
performing fits to component spin distributions? There are multiple avenues
through which a population model can fail: either the model is theoretically a
good fit and for any number of reasons (e.g. those discussed in Sec. V) cannot
find the truth, or the model is intrinsically a poor fit, i.e. it does not
have enough flexibility to find the shape of the true, underlying
distribution. Both cases induce mismatch between the true distribution and the
inferred distribution, which we hope to diagnose using only the information
available to us. In this section, we begin by discussing the first scenario
(Fig. 3) and then the second (Fig. 4).

In reality, given the complexities of astrophysical BH spin evolution, it is
almost certain that our measured distributions are in some other way
discrepant with the truth; phenomenological models likely cannot perfectly
reflect the underlying populations. Model checks on current data sets are then
used to motivate more complicated parametric models that do not suffer from
identifiable deficiencies. In parallel, nonparametric inference introduces
more flexible models that are based on a large number of parameters, however
those are also subject to model uncertainties and impose correlations across
the population parameter space [39, 41, 47]. Detailed model checking remains
an essential ingredient of population constraints.

For end-to-end event simulation and population recovery such as Fig. 2, we a
priori know what the “true” underlying astrophysical distribution is. However,
when dealing with real GW observations, this is, of course, not the case. We
therefore diagnose the bias seen in Fig. 2 using only information available to
us when dealing with real observations. To do so, we use posterior predictive
checks (PPCs) that examine the predictive accuracy of the inferred models via
its ability to predict future data that are consistent with current
observations. PPCs are ubiquitous in the field of GW population analyses [52,
65, 34, 2, 3].

We now look at the results from the Beta+DoubleGaussian model presented in
Sec. IV: a case in which a population model is theoretically a good fit, but
cannot find the underlying distibution accurately. A PPC for the
LowSpinAligned population222We highlight the LowSpinAligned population
throughout Sec. VI as it displays the largest bias under the
Beta+DoubleGaussian model. is plotted in Fig. 3. Specifically, we plot the
spin parameters predicted by the fitted model against those of the observed
events. The “predicted” (horizontal axis) and “observed” (vertical axis) draws
and are generated as follows:

  1. 1.

Draw one sample from the posterior for the Beta+DoubleGaussian hyper-
parameters.

  2. 2.

Draw one sample from the detectable [47, 66] $\chi_{i}$ and $\cos\theta_{i}$
distribution corresponding to this hyper-parameter. This is the predicted
draw.

  3. 3.

Draw one sample from one individual-event posterior in the catalog, reweighted
to the population from Step 1, as described in App. F. This is the observed
draw.

  4. 4.

Repeat 70 times for the O3-like catalog, or 300 times for the O4-like catalog.

The predicted and observed values are sorted and plotted against each other,
generating one trace in Fig. 3. We repeat this procedure 100 times to generate
a collection of traces. If we have perfectly measured the true underlying
distribution and in the limit of infinitely many observations, the traces
should be an exact diagonal. For a number of finite observations, the average
of the traces should be diagonal [67, 68, 52, 69, 32]. As the number of
observed events increases, the spread of the traces around the diagonal should
decrease.

For all spin parameters shown, the traces on average do follow the diagonal,
even though the measured population does not match the truth. The 300 event
case (navy) traces are more tightly clustered around the diagonal than the 70
event case (pink), as expected. That the traces average to the diagonal but we
know the fit is poor indicates that this class of PPC, although widely used in
GW population analyses, is not a sufficient diagnostic of model mismatch or
inaccurate population inference in this case.

Figure 4:  Results for diagnosing model misspecification for the
LowSpinAligned population’s $\cos\theta$ distribution under the Beta+Gaussian
model. Results from three different simulated measurement uncertainties are
shown: individual-event spin measurement error of $\sigma_{\mathrm{meas}}=$
0.1 (blue), 0.3 (pink), and 0.5 (yellow). The “true” measurement error from
the Bilby runs averages to $\sigma_{\mathrm{meas}}=0.48$. All results are
shown from runs done on 70-event catalogs. (Top Row) Traces corresponding to
draws from the population posterior, compared to the true population (black).
(Middle Row) PPCs from 100 catalogs, each with 70 events. In the absence of
model misspecification, the traces should on average follow the diagonal.
(Bottom Row) Fraction of events from the posterior predictive checks with
$\cos\theta$ under-predicted. The shaded regions indicate three-sigma
uncertainty on each average value, marked by the crosses. In the absence of
model misspecification, the error bars should encompass a horizontal line at
0.50.

We continue to investigate the conditions under which PPCs succeed or fail by
next turning to the second case discussed previously: that in which a
population model is an intrinsically bad fit to the underlying astrophysical
distribution. We again perform population inference of simulated data, now
deliberately using a model that cannot reproduce the injected spin
distributions. In particular, we will adopt a model in which cosine tilts are
described only as a single Gaussian, which we call the Beta+Gaussian model and
is given analytically in App. D.1. This population model is capable of
reproducing the $\chi_{\mathrm{eff}}$ distribution, but at the level of
component spins cannot capture the bimodality present in the MediumSpin and
LowSpinAligned populations.

We here wish to isolate the effect of a bad model, without having to worry
about the shortcomings of inference per Fig. 2. As such, this time we do not
perform full signal injection and parameter estimation, but instead produce
mock spin magnitude and cosine tilt posteriors, allowing us to better control
and understand the interplay between individual-event and population-level
measurement uncertainty. We assume these mock-posteriors to be Gaussian
distributed with width $\sigma_{\mathrm{meas}}$. Our procedure for generating
these mock-posteriors is detailed in App. G.2.

Results from conducting hierarchical inference using the Gaussian mock-
posteriors are shown in Fig. 4. The left-hand column (blue), shows results
from individual-event spin posteriors with $\sigma_{\mathrm{meas}}=0.1$,
between $\sim\,$1-5 times more informative than the Bilby-produced
$\cos\theta$ posteriors, which averages to $\sigma_{\mathrm{meas}}=0.48$. In
the right-hand column (orange) are plotted results from more-realistic
individual-event measurement error of $\sigma_{\mathrm{meas}}=0.5$; the center
column (pink) is an intermediate case of $\sigma_{\mathrm{meas}}=0.3$.

The inferred $\cos\theta$ distributions for the LowSpin- Aligned population
under the Beta+Gaussian model are plotted in the top row of Fig. 4. For each
of the three different individual-event measurement errors, the traces are
clustered tightly, meaning that the inferred population is precisely measured,
even though the model is not a good fit to the underlying population. The
Beta+Gaussian model is, in essence, doing its job: even with large individual
event uncertainty, it identifies the mean and the over-all width of the
distribution very well, even though it cannot capture the full underlying
bimodal structure.

We again ask: if we did not know the injected distribution, would we have been
able to tell that this model is insufficient? Going further, in the case that
we can tell a PPC fails, we are looking for an estimate of how we should
ammend our population model to better fit the truth. Beyond just inspecting
the diagonality of PPCs by eye, we can calculate the fraction of events
over/under-predicted by our model across parameter space using the slopes of
the PPC traces. If the slope of a PPC trace is steeper (shallower) than the
diagonal, then the model is predicting more (fewer) events in that region of
parameter space than are observed. To find the slopes of each trace as a
function of each parameter of interest, we perform linear regression in a
small region around each point on a grid spanning that parameter. The
fractions of each spin parameters under-predicted for each simulated
population is then the fraction of traces with slopes shallower than the
diagonal (i.e. $<1$). If the model is a good fit to the data, the fraction
under-predicted should be consistent with 0.5.

PPCs and the corresponding fraction of events under-predicted are shown in the
middle and bottom rows of Fig. 4 respectively. Errors on the fraction under-
predicted are calculated by repeating the PPC procedure ten times, and
calculating the mean (crosses) and variance (shaded region) of the results.
For the $\sigma_{\mathrm{meas}}=0.1$ case, the PPC is inconsistent with the
diagonal, meaning that here we can identify that inferred distribution under
the Beta+Gaussian model is not a good fit. The fraction of events under-
predicted is correspondingly inconsistent with 0.5. For $\cos\theta\lesssim
0.25$ and $\cos\theta\gtrsim 0.75$, the fraction under-predicted is greater
than 0.5, meaning that the model ubiquitously under-predicts the population in
this region of parameter space. Between $0.25\lesssim\cos\theta\lesssim 0.75$,
the fraction is less than 0.5, meaning here the model over-predicts. Looking
at the top left corner of Fig. 4, we can see that this is exactly the case.
These results hint at how we could improve the $\cos\theta$ model: to find the
truth we should allow the model to predict more events at alignment and anti-
alignment, i.e., include a bimodality.

As the individual-event measurement uncertainty increases (left to right), the
PPCs become more consistent with the diagonal, and correspondingly the
fractions become consistent with 0.5. By realistic measurement uncertainty, we
lose our ability to diagnose inconsistency between the underlying and measured
populations using PPCs. A crucial step in generating PPCs is the reweighting
of individual-event posteriors to the inferred population. If individual-event
posteriors are sufficiently uninformative, then this process yields reweighted
posteriors that are all essentially identical to the measured population.
Thus, the “observed” and “predicted” draws will be the same, and the PPCs will
be on average diagonal. _This type of PPC is therefore insufficient for weakly
informative parameters as the reweighted posterior is dominated by the
population rather than the individual-event likelihood._ We propose that
alternative model-checking procedures must, therefore, be developed and
utilized for diagnosing model bias and misspecification for poorly-measured
BBH parameters such as spin components.

##  VII Comparison to past work on Hierarchical Inference with Spin-Precession

To our knowledge, our study includes the first full, end-to-end individual-
event and population-level GW injection campaign for multiple distinct
populations of spin magnitudes and tilt angles of BBHs. Past studies
performing injection-recovery campaigns for spin populations either use
different waveform models and sampling implementations, and/or consider cases
with reduced complexity compared to ours. Our work is in consistent with past
findings, as described below.

Talbot and Thrane [70] investigated the measurability of the spin tilt angle
distributions alone using an astrophysically-motivated model assuming that
some fraction of BBHs form in isolated binaries, while the rest form
dynamically. They performed an injection and recovery campaign for spin tilts
where they measured the fraction of binary mergers with preferentially-aligned
versus isotropically distributed tilts, and the typical degree of spin
misalignment for each BH. In their study, all simulated binaries share the
same masses, distance, and spin magnitudes, chosen to be similiar to LIGO’s
first event GW150914 [71]. Using the waveform model IMRPhenomPv2 [72] and
nested sampling implemented in LALInference [73], they found that they are
able to constrain the parameters of the tilt-angle distributions for five
different populations. Although they do sample over all fifteen BBH parameters
in individual-event inference, Talbot and Thrane [70] is most similar to the
follow-up studies we present in App. G.3: we too are able to better recover
the underlying distributions for all three of our populations when the
complexity of the explored parameter space is reduced (see, e.g. Fig. 12), on
either an individual-event or population level.

In the context of searching for unresolved binary signals, Smith  _et al._
[74] also simulated and recovered BBH spin magnitude and tilt distributions.
They looked at a single population, consistent with the LIGO/Virgo O1 and O2
observations [38], and used the IMRPhenomPv2 waveform model [72] implemented
in Bilby for individual-event inference. A crucial difference between this
study and ours is the use of a selection function: as Smith  _et al._ [74] are
looking at resolved and unresolved binaries, they ignore selection effects
entirely. Under these conditions, Smith  _et al._ [74] find that the spin
magnitude and tilt angle distributions are both accurately measurable. Given
that their simulated populations are most similiar to our MediumSpin
population, this is in agreement with our findings, as we are able to well
constrain the MediumSpin population for both 70 and 300 events. It is only
when more complex distributions are introduced that our inference fails.

Another point of comparison between our work and others’ is on the subject of
biased measurements from population model misspecification. In particular,
other authors have also identified shortcomings of traditionally and widely
used model checking techniques such as probability-probability plots (App.
G.1) and posterior predictive checks (Sec. VI). Biscoveanu  _et al._ [64]
discussed population model bias in the mass distribution of binary neutron
star populations arising from misspecification of spin priors. As part of
their work, they show that a P-P check on individual-event posteriors can pass
but still lead to highly biased population inference. This is in agreement
with our findings.

##  VIII Conclusions

In this work, we investigated the measurability of the spin magnitude and tilt
angle distributions of BBH populations via GW observations. To see if
realistic GW populations contain information about spin components or just the
effective spin, we simulated three BBH populations that have the same
underlying effective-spin distributions, but deliberately distinct spin
magnitude and tilt distributions, and on them conducted individual-event and
population-level parameter estimation. We then turned to the question of
whether mismatch between the injected and recovered spin magnitude and tilt
distributions can be identified using only the individual-event and
population-level data available to us, without knowledge of the true
underlying population. Our work focuses on the three questions posed in Sec.
I, the answers to which we summarize below.

(1) There is information in gravitational-wave signals beyond the effective-
spin. As discussed in Sec. IV, we can tell that our three different
populations have different spin magnitude and tilt distributions despite their
having identical $\chi_{\mathrm{eff}}$ distributions.

(2) Measuring component spin distributions accurately is practically
challenging. Under standard, reviewed parameter estimation settings, we were
able to accurately measure the spin magnitude and/or tilt angle distributions
for some of the populations, but not all three. The bimodal tilt distribution
of the LowSpinAligned population proved especially resistant to being
accurately constrained, even when using a population model that was inherently
bimodal. We employed a suite of verification methods to ensure the robustness
of these results, which are enumerated in Sec. V. Notably, we find that the
effective spin distribution is, however, accurately measured no matter the
degree of mismatch in the component spin model constraints; this is not true
for the effective precessing spin $\chi_{\mathrm{p}}$ which remains
susceptible to biased spin magnitude and tilt inference.

Although we cannot say for certain, we hypothesize that the root of the
mismatch between the recovered and underlying distributions is related to a
lack of convergence in individual-event posteriors, the specifics of which are
subtle enough to not present themselves via a standard P-P check (App. G.1).
We do not claim that accurately recovering component spin population
distributions is impossible at current sensitivity, just challenging. Running
Bilby with more aggressive settings, while computationally costly, may very
well fix the problems presented in this work. However, unlike our injection
set, real observations do not come with an answer key. If manually tuning
sampler settings is a requirement to recover truth, we must be aware what
these same errors could manifest in real LIGO/Virgo events.

(3) At current sensitivity, we cannot tell when measurements of component spin
distributions are biased via the currently widely-used method of posterior
predictive checks. Due to the fact that individual-event posteriors are
extremely weakly informative about spin components, reweighting these
posteriors to the inferred population distribution–a crucial step in
conducting posterior predictive checks–yields individual-event measurements
that are all nearly identical to the inferred population itself. Nearly any
population model can seem like a good fit to poorly constrained data, as
discussed in Sec. VI.

Fishbach  _et al._ [52] detailed different categories of posterior predictive
checking for GW data. The most commonly used level is what we do in this work:
performing consistency checks on the true underlying parameters of the
observed data versus predicted by the model. However, one can also conduct
PPCs on the observed parameters (e.g. max likelihood parameters) of the data
versus those predicted by the model. While checks on the true parameters are
susceptible to the issues related to reweighting that we discuss in Sec. 3,
checks on observed parameters might be more constraining. However, they are
far more computationally expensive to perform, as one must generate maximum
likelihood values predicted by the model: this involves either running an
optimization routine or conducting mock-parameter estimation on thousands of
events. While trustworthy mock-parameter estimation exists for some parameters
(e.g. masses) [52, 53, 54], the imprint of spin magnitudes and tilt angles on
data is more subtle and remains unincorporated into these algorithms.
Developing different, more-informative methods of posterior predictive
checking for poorly-constrained parameters such as spin is an essential topic
of future work.

## Data and Code Availability

The code used to produce all results presented in this paper can be found at
https://github.com/simonajmiller/measuring-bbh-component-spin. Our individual-
event and hierarchical-inference posteriors samples can be shared upon
request.

###### Acknowledgements.

We thank Sylvia Biscoveneau, Jacob Golomb, Ethan Payne, and Colm Talbot for
their extensive help with Bilby parameter estimation and for helpful feedback
on this work, and Sophie Hourihane for essential insights about reweighting
and probability-probability plots. Furthermore, we extend thanks to Maya
Fishbach, Reed Essick, Matthew Mould, Zoheyr Doctor, and Will Farr for their
suggestions related to the range of verification methods performed in this
paper. This work is supported by National Science Foundation grant No.
PHY-2150027 as part of the LIGO Caltech REU Program. SM and KC were supported
by NSF Grant PHY-2110111 and NSF Grant PHY-2308770. The authors are grateful
for computational resources provided by the LIGO Laboratory and supported by
National Science Foundation Grants PHY-0757058 and PHY-0823459. Software:
emcee [55], bilby (version 2.2.2) [50, 51], dynesty (version 2.1.2) [49],
numpy [75], scipy [76], matplotlib [77], seaborn [78], astropy [79, 80], jax
[81], numpyro [56, 57]

##  Appendix A Simulated Populations

We simulate three populations with the same $\chi_{\mathrm{eff}}$ but
different spin magnitude $\chi$ and tilt angle $\theta$ distributions. To
generate the populations, we first choose distributions of the mass ratio $q$
and spin $z$-component $s_{i,z}\equiv\chi_{i}\cos\theta_{i}$ to be shared in
common across all three populations; this ensures the same
$\chi_{\mathrm{eff}}$ distribution. The mass ratio distribution corresponds to
the median posterior value inferred with the PowerLaw+Peak model in Ref. [3],
while for $s_{i,z}$ we select a Gaussian with mean $0.10$ and standard
deviation $0.15$. This yields a Gaussian-like $\chi_{\mathrm{eff}}$
distribution with mean $0.10$ and standard deviation $0.11$.

To decompose $\chi_{\mathrm{eff}}$ into component spins, we choose a different
spin magnitude $\chi$ distribution for each population and then numerically
calculate the resultant $\cos\theta$ distribution implied by $p(\chi)$ and
$p(s_{i,z})$. This procedure results in the three populations shown in Fig. 1.
The $\chi$ distribution for the HighSpinPrecessing population is uniform
between $s_{i,z}$ and $1$, i.e., $s_{i,z}$ values are drawn from the Gaussian
described above and then a $\chi_{i}$ value is conditionally drawn based on
each $s_{i,z}$. For the MediumSpin (LowSpinAligned) population, each $\chi$
value is drawn from a Gaussian distribution about $s_{i,z}$, truncated on
$0\leq\chi\leq 1$, with a standard deviation of $0.20$ ($0.05$). For each
population, we assume $\chi_{1}$ and $\chi_{2}$ are identically but
independently distributed, as are $\cos\theta_{1}$ and $\cos\theta_{2}$.
Finally, each spin vector’s azimuthal angle $\phi_{i}$ is drawn uniformly
between $0$ and $2\pi$. Due to the different $\chi$ and $\cos\theta$
distributions, the $\chi_{p}$ distributions of each population differ as
well.333Different component spin distributions given identical
$\chi_{\mathrm{eff}}$ _and_ $\chi_{\mathrm{p}}$ distribution can only be
achieved by relaxing the assumption of identically distributed component spin
magnitudes and angles.

The astrophysical distribution of the remaining binary parameters is the same
for all populations. We inject primary masses drawn from the PowerLaw+Peak
model [82] with all parameters, except for $m_{\mathrm{min}}$, fixed to their
one-dimensional median values as found in Ref. [3]: $\alpha=3.51$,
$m_{\mathrm{max}}=88.21$, $\lambda_{\mathrm{peak}}=0.033$, $\mu_{m}=33.61$,
$\sigma_{m}=4.72$, and $\delta_{m}=4.88$ in the notation used therein. The
injected mass ratio distributions for all populations are described by a power
law with slope $\beta_{q}=0.96$ (see Eq. 12), again the median inferred value
from [3]. In the parameterization of the PowerLaw+Peak model, we use a
population minimum mass of $m_{\mathrm{min}}=6\,M_{\odot}$ instead of
$5\,M_{\odot}$ to set the shape of the distribution. We additionally impose a
mass cut of $8\,M_{\odot}$, as restricting to higher-mass events ensures
shorter analysis times. This mass cut effectively becomes the minimum mass,
but we renormalize the distribution to keep the same shape above the cutoff
mass as it would with $m_{\mathrm{min}}=6$. Explicitly setting
$m_{\mathrm{min}}=8$ in the PowerLaw+Peak model would change the over-all
shape of the distribution to be inconsistent with the desired results in Ref.
[3].

Finally, the BBH merger density rate in the source frame evolves with respect
to redshift $z$ as

$R(z)\propto\frac{dV_{c}}{dz}\left(1+z\right)^{2.7}\,,$ (3)

where $V_{c}$ is the comoving volume. Distances are calculated from redshifts
assuming the cosmology reported by the Planck 2013 survey [83]. All other
parameters are drawn uniformly from their respective physical range. Mass and
redshift distributions are plotted in Fig. 5.

##  Appendix B Individual-event Parameter Estimation

Figure 5: The underlying (navy) and detected (orange) distributions for
source-frame primary mass $m_{1}$ (top), source-frame secondary mass $m_{2}$
(middle), and redshift $z$ (bottom) shared between all three simulated
populations.

From each of the three astrophysical distributions described in App. A, we
draw $10^{5}$ events. We apply a network SNR [84, 85, 86, 87] cut of 10 in the
LIGO Livingston, LIGO Hanford, and Virgo detector network using the “O3
actual” power spectral densities (PSDs) provided in Ref. [48], and select 300
detectable events. Histograms of events from the underlying (navy) versus
detectable (orange) mass and redshift distributions are shown in Fig. 5.

For each event, we simulate GW data with the IMRPhenomXPHM waveform model [27]
including a Gaussian noise realization and draw samples from the
15-dimensional posterior distribution for the binary parameters using the same
waveform. We employ standard priors for all binary parameters [51], although
we use a targeted chirp mass [88, 89] prior of $\pm 15M_{\odot}$ about the
injected value to reduce computational cost, which we verify does not affect
results. Simulated data assume a detector network of LIGO-Hanford, LIGO-
Livingston [4], and Virgo [5], each at their O3 sensitivity [48] with a
sampling rate of $2048~{}\mathrm{Hz}$. We analyze data in the
$15-921.6(=0.9\times 2048/2)~{}\mathrm{Hz}$ frequency range, assuming perfect
knowledge of the detector calibration.

We use the nested sampler dynesty [49] as implemented in Bilby [50, 51] under
reviewed settings to stochastically sample from the individual-event
posteriors. For sampler settings, we use nlive = 1000, naccept= 60, and sample
= "acceptance-walk". Time marginalizaton is turned on, while distance and
phase marginalization remain off. Post-facto, we apply an optimal SNR cut of
10 on the posterior samples for consistency with our selection criteria [66].

##  Appendix C Hierarchical Inference

The hierarchical inference framework used in our analysis to obtain posteriors
distributions on the population parameters is implemented using the Python
Markov-Chain Monte-Carlo package emcee [55]. The likelihood
$\mathcal{L}(\\{d\\}|\Lambda)$ that a catalog of $N_{\mathrm{obs}}$ GW events
with data $\\{d_{i}\\}_{i=1}^{N_{\mathrm{obs}}}$ arises from an underlying
population $\pi_{\mathrm{pop}}$ described by parameters $\Lambda$ is given by
[58, 59, 90, 60]

$\mathcal{L}(\\{d\\}|\Lambda)\propto\prod^{N_{\mathrm{obs}}}_{i}\frac{\int
d\lambda_{i}\,\mathcal{L}(d_{i}|\lambda)\pi_{\mathrm{pop}}(\lambda|\Lambda)}{\xi(\Lambda)}\,,$
(4)

where $\lambda_{i}$ are the parameters of the $i^{\mathrm{th}}$ event in the
catalog (i.e. spins, masses, etc). In practice, we have access to the
individual-event posteriors $p(\lambda_{i}|d_{i})$ obtained with a default
parameter-estimation prior $\pi_{\mathrm{pe}}(\lambda_{i})$, rather than the
event likelihood $\mathcal{L}(d_{i}|\lambda_{i})$. We thus write Eq. (4) as

$\mathcal{L}(\\{d\\}|\Lambda)\propto\xi(\Lambda)^{-N_{\mathrm{obs}}}\prod_{i}\int
d\lambda_{i}\,\frac{p(\lambda|d_{i})}{\pi_{\mathrm{pe}}(\lambda_{i})}\pi_{\mathrm{pop}}(\lambda_{i}|\Lambda)\,.$
(5)

Additionally, rather than $p(\lambda_{i}|d_{i})$ itself, we have a discrete
set of $N_{i}$ independent samples $\\{\lambda_{i,j}\\}_{j=1}^{N_{i}}$ drawn
from $p(\lambda_{i}|d_{i})$. Using the standard procedure, we approximate the
integral of Eq. (5) via a Monte Carlo average,

$\mathcal{L}(\\{d\\}|\Lambda)\propto\xi(\Lambda)^{-N_{\mathrm{obs}}}\prod_{i}\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\frac{\pi_{\mathrm{pop}}(\lambda_{i,j}|\Lambda)}{\pi_{\mathrm{pe}}(\lambda_{i,j})}\,.$
(6)

The detection efficiency

$\xi(\Lambda)=\int
d\lambda\,\pi_{\mathrm{pop}}(\lambda|\Lambda)P_{\mathrm{det}}(\lambda)\,,$ (7)

is the fraction of events that we would successfully detect if the population
with parameters $\Lambda$ is the true underlying population. Here,
$P_{\mathrm{det}}(\lambda)$ is the probability that an individual event with
parameters $\lambda$ is detected. As with the population likelihood, we
calculate the detection efficiency with a Monte Carlo average. Given
$N_{\mathrm{inj}}$ injected signals drawn from some reference distribution
$p_{\mathrm{inj}}(\lambda)$, the detection efficiency is

$\xi(\Lambda)=\frac{1}{N_{\mathrm{inj}}}\sum_{i=1}^{N_{\mathrm{fnd}}}\frac{\pi_{\mathrm{pop}}(\lambda_{i}|\Lambda)}{p_{\mathrm{inj}}(\lambda_{i})}\,,$
(8)

where the sum is over the $N_{\mathrm{fnd}}$ injections that pass the
detection criteria. We generate the set of “found” injections over which the
Monte Carlo average is calculated in the same way that we produced catalogs of
events in App. B. The reference $p_{\mathrm{inj}}(\lambda)$ follows the true
mass and redshift distribution (App. A; Fig. 5), but is uniform in spin
magnitudes and isotropic in spin tilts such that we can resolve features
across the full underlying spin distribution. As in App. B, our detection
criterion is an optimal SNR greater than 10 using the waveform IMRPhenomXPHM
[27] in the LIGO Livingston, LIGO Hanford, and Virgo network at O3 sensitivity
[48]. We acknowledge that the optimal SNR is not a strictly accurate estimate
of selection effects on real data as it is solely a function of source
parameters, not detector noise. However, this approach remains formally self-
consistent as long as we apply the same optimal SNR cut on posterior samples,
as explained in Essick and Fishbach [66].

Following Farr [91], we account for uncertainty in the Monte Carlo integral by
demanding that the effective number of independent samples

$N_{\mathrm{eff}}(\Lambda)\equiv\frac{\left[\sum_{i=1}^{N_{\mathrm{fnd}}}w_{i}(\Lambda)\right]^{2}}{\sum_{i=1}^{N_{\mathrm{fnd}}}\left[w_{i}(\Lambda)\right]^{2}}\geq
4\,N_{\mathrm{obs}}\,,$ (9)

where the weights $w_{i}$ between the population distribution and parameter
estimation prior are defined as

$w_{i}(\Lambda)=\frac{\pi_{\mathrm{pop}}(\lambda_{i}|\Lambda)}{p_{\mathrm{inj}}(\lambda_{i})}\,,$
(10)

evaluated on the parameters of the found injections. This procedure rejects
samples from regions of parameter space in which there are not sufficient
injections to accurately probe. We use 200,000 injections to calculate $\xi$.
Our results never rail against the $N_{\mathrm{eff}}$ cut of Eq. (9), so we do
not believe it affects our results. For further investigation, we perform a
set of analyses without including spins when calculating $N_{\mathrm{eff}}$
(see App. G.7), under which our conclusions do not change.

In addition to including a cut on effective samples from the selection
function, one can also impose a cut on the per-event effective samples of the
posteriors used in calculating the hierarchical likelihood. Here, instead of
evaluating Eq. 10 on the found injections, it is evaluated on the posterior
samples for every event. If any events in the catalog have an effective sample
number below some threshold, the corresponding $\Lambda$ sample is tossed. In
this work, we do not include any per-event $N_{\mathrm{eff}}$ cuts in the
sampling of $\mathcal{L}(\\{d\\}|\Lambda)$, but calculate them post-facto as a
check of Monte Carlo convergence (see e.g. Fig. 16). Other tests of Monte
Carlo convergence are discussed in [92].

##  Appendix D Spin population models

We recover the simulated populations with two models. Generically, we
factorize the population models as

$\begin{split}\pi_{\mathrm{pop}}(\lambda|\Lambda)&=p(m_{1}|\Lambda)\,p(m_{2}|m_{1},\Lambda)\,p(z|\Lambda)\times\\\
&\,p(\chi_{1}|\Lambda)\,p(\chi_{2}|\Lambda)\,p(\cos\theta_{1}|\Lambda)\,p(\cos\theta_{2}|\Lambda)\,,\end{split}$
(11)

meaning, aside from the masses $m_{1}$ and $m_{2}$, there are no correlations
in the population. Our two parameterized models for the spin magnitude $\chi$
and tilt angle $\theta$ are described below. More details about the two models
are provided in Table 1. Spin magnitudes $\chi_{i}$ and tilt angles
$\theta_{i}$ are assumed identically and independently distributed.

During hierarchical inference, we fix the distributions of primary mass
$m_{1}$ and redshift $z$ to truth, as described in App. A. To account for
possible correlations between spins and mass ratio, although no underlying
correlation was injected, we follow Callister  _et al._ [93] and
simultaneously infer the distribution of binary mass ratios and spins using a
secondary mass distribution of

$p(m_{2}|m_{1})\propto m_{2}^{\beta_{q}}\quad\left(m_{\mathrm{min}}\leq
m_{2}\leq m_{1}\right)\,,$ (12)

where the power-law index $\beta_{q}$ is a free parameter with a Gaussian
prior of $\mathcal{N}(0,3)$. The true underlying distribution has
$\beta_{q}=0.96$. For all other individual-event parameters, we take
population distributions identical to the priors used during the original
Bilby parameter estimation. Most notably for this analysis, azimuthal spins
are distributed uniformly $\phi_{i}\in[0,2\pi)$. All parameters aside from
masses, redshift, and spin magnitudes and tilt angles are thus excluded from
hierarchical inference.

Model Name |  $p(\chi)$ $p(\cos\theta)$ | Parameter | Prior | Comments
---|---|---|---|---
Beta+Gaussian |  | $\mu_{\chi}$ | U(0,1) | Component spin model that cannot reproduce the simulated populations; used to study model misspecification
$\sigma_{\chi}$ | U(0.07,0.5)
$\mu_{\theta}$ | U(-1,1)
$\sigma_{\theta}$ | U(0.16,0.8)
Beta+DoubleGaussian |  | $\mu_{\chi}$ | U(0,1) | Component spin model that can reproduce the simulated populations; used to study the amount of information available in component spins. See App. G.7 about methods of breaking the degeneracy between the two Gaussian components.
$\sigma_{\chi}$ | U(0.07,0.5)
$\mu_{\theta,1}$ | U(-1,1)
$\sigma_{\theta,1}$ | U(0.16,0.8)
$\mu_{\theta,2}$ | U(-1,1)
|  | $\sigma_{\theta,2}$ | U(0.16,8) |
|  | $f$ | U(0,5) |
Table 1: Details about the two component spin models we employ. Columns give
the model names, an example $\chi$ and $\cos\theta$ plot, the parameters they
depend on, the parameter priors, and some brief comments. The notation
U$(a,b)$ means a prior uniform between $a$ and $b$. For $p(\chi)$ in both
models, we impose an additional prior requirement that the beta distribution
shape parameters $\alpha,\beta>1$, see Eq. (14), making the distribution non-
singular at the boundaries. Full expressions for $p(\chi)$ and $p(\cos\theta)$
can be found in Eqs. (13), (16), and (17).

###  D.1 Beta+Gaussian

Following the Default model in Ref. [3], we assume that spin magnitudes
$\chi_{i}$ are identically and independently distributed according to a Beta
distribution

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

where $B(\alpha,\beta)$ is the Beta function which ensures that the
distribution is normalized to unity on $0\leq\chi\leq 1$. Instead of sampling
in the shape parameters $\alpha$ and $\beta$, we sample in the more familiar
mean $\mu_{\chi}$ and standard deviation $\sigma_{\chi}$ which are related to
$\alpha$ and $\beta$ by

$\alpha=\mu_{\chi}\nu~{},\quad\beta=(1-\mu_{\chi})\nu\,,$ (14)

where

$\nu=\frac{\mu_{\chi}(1-\mu_{\chi})}{\sigma_{\chi}^{2}}-1\,\,.$ (15)

We adopt uniform priors on $\mu_{\chi}$ and $\sigma_{\chi}$ and impose an
additional cut such that $\alpha,\beta\geq 1$ to keep the distribution
bounded. This cut enforces $p(\chi_{i}|\alpha,\beta)=0$ at $\chi_{i}=0$ and
$1$. Aside from the HighSpinPrecessing population at $\chi=1$, this assumption
is valid, and even there the spin model captures the over-all distribution’s
shape.

For the tilt-angle distribution, we adopt a truncated, normalized Gaussian
distribution

$p(\cos\theta_{i}|\mu_{\theta},\sigma_{\theta})=\mathcal{N}_{[-1,1]}(\cos\theta_{i}|\mu_{\theta}\,,\sigma_{\theta})$
(16)

on the interval $-1\leq\cos\theta\leq 1$, and fit for the mean $\mu_{\theta}$
and standard deviation $\sigma_{\theta}$.

###  D.2 Beta+DoubleGaussian

The Beta+DoubleGaussian model uses the same spin magnitude distribution as the
Beta+Gaussian, as given in Eq. (13) and explained thereafter. The tilt angle
distribution is here instead given by a mixture of two truncated normalized
Gaussians

$p(\cos\theta_{i}|\mu_{\theta,1},\sigma_{\theta,1},\mu_{\theta,2},\sigma_{\theta,2},f)=\\\
f\mathcal{N}_{[-1,1]}(\cos\theta_{i}|\mu_{\theta,1},\sigma_{\theta,1})\\\
+(1-f)\mathcal{N}_{[-1,1]}(\cos\theta_{i}|\mu_{\theta,2},\sigma_{\theta,2})\,,$
(17)

to capture the multimodality of the some of the underlying distributions. We
measure the means $\mu_{\theta,1}$, $\mu_{\theta,2}$ and standard deviations
$\sigma_{\theta,1}$, $\sigma_{\theta,2}$ of the two Gaussians, and the mixing
fraction $f$ between them. We impose $\mu_{\theta,1}\leq\mu_{\theta,2}$ to
distinguish between the two components.

##  Appendix E Detailed Hierarchical Inference Results

We here show full posteriors on the hyper-parameters for the
Beta+DoubleGaussian component spin population model, as given in Eqs. (13) and
(17), for the HighSpinPrecessing (Fig. 6), MediumSpin (Fig. 7), and
LowSpinAligned (Fig. 8) populations. Results for 70 (pink) and 300 (navy)
event catalogs are shown in each figure. These posteriors are compared against
non-linear least-squares fit parameters (black dashed) calculated from 50,000
draws per population, representing the best possible fit for the true
underlying distributions within the Beta+DoubleGaussian model. Population
distributions generated from draws from these posteriors are plotted in Fig.
2. For all three populations, as expected, including more events makes hyper-
parameter measurements more precise. However, adding more events does not
necessarily make the results more accurate.

The hyper-parameters of the HighSpinPrecessing population (Fig. 6) are
recovered with minimal bias. While the mean $\mu_{\chi}$ of the spin magnitude
distribution of the HighSpinPrecessing population is very well constrained,
its width $\sigma_{\chi}$ is slightly underestimated in the case of both the
70 and 300 event catalogs. The means $\mu_{i,\cos\theta}$ are also accurately
constrained. The widths of the tilt angle distributions also seem to be under-
estimated, but per Fig. 2, the actual shape of the distribution converges on
the truth. This is because – due to the allowed bimodality – different
combinations of hyper-parameters can lead to the same unimodal distribution.

The MediumSpin population (Fig. 7), on the other hand, is reconstructed very
accurately by the Beta+DoubleGaussian model for both the 70 and 300 event
catalogs. Each “true” hyper-parameter either falls within the 90% measured
credible region. This is also reflected in Fig. 2 – the black traces
representing truth are enclosed by the 90% credible envelopes for both $\chi$
and $\cos\theta$.

Finally, the width of the spin magnitude distribution for the LowSpinAligned
population (Fig. 8) is accurately constrained, but its mean $\mu_{\chi}$ is
over-estimated. The most striking failure of the Beta+DoubleGaussian model is
its inability to identify the bimodality of the LowSpinAligned population’s
tilt angle distribution, for either 70 or 300 events. Aside from the mixing
fraction $f$, none of the tilt-angle distribution hyper-parameters’ posteriors
are consistent with truth at the 90% level.

For all three populations, the power law slope for the mass distribution is
recovered within 90% credibility about the injected value of $\beta_{q}=0.96$.
Masses are recovered without bias by our hierarchical inference procedure; we
only encountering biases when fitting for the spin populations.

Figure 6:  The posterior distributions on the hyper-parameters of the spin
magnitude and tilt angle distributions under the Beta+DoubleGaussian model for
the HighSpinPrecessing population for 70 (pink) and 300 (navy) event catalogs.
The labels above each one-dimensional posterior give the medians and 90%
credible intervals on each hyper-parameter for the two different catalog
sizes, while the contours in each two-dimensional posterior denote the 50% and
90% credible regions. See Table 1 for descriptions the hyper-parameters and
their priors. Black dashed lines labeled “truth” represent the theoretical
best-fit parameters for the population under the Beta+DoubleGaussian model, as
calculated using a least-squared fit on 50,000 draws from the population.
Figure 7:  Same as Fig. 6 but for the MediumSpin population. This population
is recovered by the Beta+DoubleGaussian model without bias.  Figure 8:  Same
as Figs. 6 and 7 but for the LowSpinAligned population. This population is
recovered by the Beta+DoubleGaussian model with considerable bias.

##  Appendix F Reweighting individual event posteriors

Given a set of discrete sample from an individual-event posterior
$p(\lambda|d_{i})$ calculated with prior $\pi_{\mathrm{pe}}(\lambda)$ and
discrete samples of hyper-parameters describing a population distribution
$\pi_{\mathrm{pop}}(\lambda|\Lambda)$, we can reweight the individual-event
posterior to the inferred population using a two-step algorithm [32, 69].
First, randomly select a hyper-parameter sample from the population
distribution $\Lambda_{i}\in\\{\Lambda\\}$ and calculate the following weights
for each individual-event posterior sample $\lambda_{j}$:

$w_{j}\propto\frac{\pi_{\mathrm{pop}}(\lambda_{j}|\Lambda_{i})}{\pi_{\mathrm{pe}}(\lambda_{j})}\,\,.$
(18)

Second, select one sample $\lambda_{j}\in\\{\lambda_{j}\\}$ subject to the
weights $w_{j}$. Repeat this process to build up a set of samples from a
reweighted individual-event posterior. This procedures ensures that events are
not double-counted during weighting [69].

##  Appendix G Verification of methods

To explore the origin of the bias observed in Fig. 2, we perform a number of
explorations that we elaborate upon in subsequent subsections. In App. G.1, we
generate probability-probability (P-P) plots for reweighted individual event
Bilby posteriors, which return unbiased. In App. G.2, we turn to hierarchical
inference and explore a range of simulated Gaussian individual-event spins
posteriors, rather than ones generated via stochastic sampling with Bilby. We
also reduce the complexity of the parameter spaces explored on both an
individual-event and population level. In App. G.3, we present results where
we fix various combinations of parameters to their true values in both the
individual-event and hierarchical inference levels. We here also discuss using
a less complex waveform – IMRPhenomXP rather than IMRPhenomXPHM – which
excludes higher order modes, in individual-event sampling. App. G.4 uses an
alternate hierarchical inference code, implemented in Numpyro instead of
emcee, and App. G.5 shows results from simultaneously inferring for the mass
and redshift distributions along with the spins. In App. G.6, we look at
hierarchically-inferred rates across parameter space rather than probability
density functions to ensure that the normalization is not obscuring the
results. Finally, other miscellaneous checks for the hierarchical inference
framework and implementation are: excluding selection effects in spin; trying
different methods of breaking the degeneracy in the double-Gaussian tilt
distribution; and looking at different 70-event catalog instantiations. Plots
showing these results can be found in App. G.7.

###  G.1 Probability-probability (P-P) plots

A crucial assumption of hierarchical inference is that the input individual-
event posteriors are themselves reliable. To test this assumption and ensure
that the stochastically-sampled Bilby individual-event posteriors (see App. B)
are indeed unbiased, we perform the common diagnostic check of generating a
probability-probability (P-P) plot [62, 63].

A P-P plot is generated by performing parameter estimation on events with
parameters distributed according to their individual-event priors, in Gaussian
noise. The percentiles, or credible intervals (CI), at which the injections
fall in their resultant one-dimensional, marginalized posterior distribution
shall be uniformly distributed if parameter estimation is unbiased. In our
case, where the injected distribution does not match the priors used in
parameter estimation, reweighting (see App. F) to the injected distribution
must be performed as a post-processing step. Specifically, we apply an optimal
SNR cut of 10 to the posteriors, and then reweight to the underlying
population; this procedure is analogous to not applying any SNR cut and
reweighting to the detected distribution.

P-P plots for spin magnitudes, spin tilt angles, masses, and redshifts for the
300 injections per simulated population are shown in Fig. 9. On the horizontal
axis, the (sorted) CIs are plotted. The vertical axis shows the cumulative
density function (CDF) of these CIs, i.e. the frequency at which each CI
occurs. This should be a diagonal line with a slope of 1 in the case of
infinitely many injections: e.g. 20% of the time, the injection should fall
within the lower 20% CI of its posterior. In the case of finitely many
injections, these CDFs should roughly fall within a 3-$\sigma$ region around
the diagonal, the width of which is a function of the number of events
injected, as indicated by the gray lines in Fig. 9.

To check if the Bilby posteriors pass the P-P test, we look at the
$p$-values444We calculate $p$-values using a Kolmogorov-Smirnov test. that
each set of $y$-axis values shown in Fig. 9 is uniform. Then, we take the
$p$-values of these $p$-values, which should also be uniformly distributed if
sampling error is random. The $p$-values (listed in the titles of Fig. 9) for
each of the three simulated populations is above the threshold of randomness
expected from the 7 parameters plotted ($1/7\sim 0.143$), indicating that the
Bilby posteriors pass the P-P test.

Figure 9:  Probability-probability (P-P) plots for spin magnitudes (dark
blue), spin tilt angles (light blue), masses (pink) and redshifts (orange) for
each simulated population (from left to right: HighSpinPrecessing, MediumSpin,
LowSpinAligned). For 300 events per population, confidence intervals (CIs;
horizontal axis) are plotted against the fraction of events for which the
true, injected value is recovered in that CI, stochastically sampled using the
Bilby implementation of the nested sampler Dynesty. The 1-, 2-, and 3-$\sigma$
regions for 300 events are plotted in gray; all parameters stay within the
3-$\sigma$ region, corresponding to the outermost gray lines. The $p$-values
for each of the three populations are greater than the threshold for 7
parameters ($\sim 0.143$), indicating that the P-P test is passed.

Though a necessary check, diagonal P-P plots are not a sufficient condition
for reliable individual-event posteriors. As also seen in [64], a sampling
algorithm can pass a P-P test but still result in biased hierarchical
inference recovery beyond the 3-$\sigma$ level. In Fig. 8, for example, the
truth for the means of the spin magnitude distribution and both modes of the
spin tilt distribution all lie outside of the 90% credible interval for the
recovered values. It remains unclear in our case whether such discrepancies
between the true and recovered populations are due to individual-event
sampling issues that are not picked up by the P-P test, as was the case in
[64], or further unknown biases.

###  G.2 Simulated Gaussian individual-event spin posteriors

To better understand the relation between individual-event and population-
level measurement uncertainty, we generate a series of simulated individual-
event spin magnitude and tilt angle posteriors for each of the same 300 GW
events per population that we stochastically sample in Bilby. We take these
mock posteriors to be Gaussian distributed with width
$\sigma_{\mathrm{meas}}$.

First, we generate a series of mock posteriors without any underlying spin-
spin correlations with the following steps. For each of the 300 injections per
population,

  1. 1.

Take the true, injected value of each spin parameter

$\lambda_{\mathrm{true}}\in\\{\chi_{1},\chi_{2},\cos\theta_{1},\cos\theta_{2}\\}\,,$

and from it draw an observed maximum likelihood value $\lambda_{\mathrm{obs}}$
from the Gaussian distribution
$\mathcal{N}(\lambda_{\mathrm{true}},\sigma_{\mathrm{meas}})$ with mean
$\lambda_{\mathrm{true}}$ and width $\sigma_{\mathrm{meas}}$.

  2. 2.

Draw $N$ samples from
$\mathcal{N}_{[a,b]}(\lambda_{\mathrm{obs}},\sigma_{\mathrm{meas}})$ where $N$
is the number of samples in the Bilby posterior for the injection of interest
and $\mathcal{N}_{[a,b]}$ is a Gaussian distribution truncated on $[a,b]$. For
spin magnitude this truncation is between $[0,1]$, and for the cosine tilt
angle $[-1,1]$.

Specifically, we look at cases where $\sigma_{\mathrm{meas}}=0.1$, $0.3$, and
$0.5$, as shown in Fig. 3. In all cases, we keep the Bilby mass and redshift
posteriors. Moreover, the simulated and Bilby posteriors all have the same
number of samples per event.

To simulate a more realistic case, we also generate a set of mock Gaussian
posteriors that do include underlying inter-spin correlations, with the same
covariance as Bilby individual-event posteriors. For each injection, we first
find the covariance of the corresponding four-dimensional Bilby posterior for
$\\{\chi_{1},\chi_{2},\cos\theta_{1},\cos\theta_{2}\\}$. We then generate a
mock four-dimensional spin posterior with that same covariance using the
procedure enumerated above: from truth, draw an observed maximum likelihood
value; then generate a posterior by sampling a truncated Gaussian centered at
that observed value. The only difference is, instead of separately generating
each one-dimensional posterior for magnitudes and tilts, we generate a four-
dimensional posterior that includes correlations.

Figure 10:  Inferred distributions obtained with the Beta+DoubleGaussian model
for spin magnitude $\chi$, spin tilt $\cos\theta$, effective
$\chi_{\mathrm{eff}}$ spin, effective precessing spin $\chi_{\mathrm{p}}$, for
the three simulated populations with simulated Gaussian individual-event spin
posteriors. Results are shown from 70 event catalogs with per-event spin
magnitude and tilt measurement error $\sigma_{\mathrm{meas}}=0.5$ (light blue
dotted), $\sigma_{\mathrm{meas}}=0.1$ (blue dashed), and with realistic
$\sigma_{\mathrm{meas}}$ and inter-parameter correlations taken from the Bilby
posteriors (pink solid). The shaded region denotes 90% of the probability,
while the black solid line corresponds to the true population.

The Beta+DoubleGaussian population model is able to recover the underlying
populations when using the simulated Gaussian posteriors, as seen in Fig. 10.
This is true for the most-informative individual-event mock-posteriors
($\sigma_{\mathrm{meas}}=0.1$; light blue dotted), the least informative
($\sigma_{\mathrm{meas}}=0.5$; blue dashed), and the most realistic (the
posteriors with correlations, labeled “realistic $\sigma_{\mathrm{meas}}$”;
pink solid). In all three cases, the true population lies within the 90%
credible interval of the recovered region. As the measurement error decreases,
the constraints get tighter around truth.

###  G.3 Reducing complexity of the explored parameter spaces

In this section, we present two simplified scenarios for individual-event
sampling in Bilby, then one simplified scenario for population inference in
emcee.

On the side of individual-event inference, we first reduce complexity by
reducing the number of parameters sampled over. The results shown in the main
text use posteriors where all fifteen dimensions of parameter-space are
sampled over (see App. B). To simplify, we first conduct parameter estimation
on all the same injections, but fixing their extrinsic parameters (i.e.
everything aside from masses and spins) to the true, injected values. This
yields the population constraints shown in the pink solid lines in Fig. 11
(for 100 events), labeled “fixed extrinsic parameters.” We then simplify
further and additionally fix the masses and spin azimuthal angles to truth,
generating the orange solid lines in Fig. 11 and labeled “fixed extrinsic
parameters + masses.” Notably, the “fixed extrinsic + masses” individual-event
posteriors yield population constraints that are significantly improved from
those shown in Fig. 2. Since sampling convergence is more challenging as the
dimensionality of the explored parameter-space increases, the trend we observe
suggests that convergence might at least partially contribute to the bias.

Next, we return to sampling over all parameters (masses, spins, and
extrinsic), but this time with a simpler waveform model: IMRPhenomXP [27].
Coming from the same family as IMRPhenomXPHM, the IMRPhenomXP model does not
contain higher order modes, which help break degeneracies between BBH
parameters. The yellow (light blue) dashed lines in Fig. 11 show the
population constraints from the same 70 (300) events as Fig. 2 but with
individual-event posteriors sampled with IMRPhenomXP. The recovered
HighSpinPrecessing and MediumSpin populations have a worse mismatch with the
truth than in Fig, 2, but the LowSpinAligned population is recovered
marginally better. Higher order modes become more important to accurately
constrain BBH parameters as the degree of spin precession increases. Thus, the
fact that the HighSpinPrecessing population is the worst constrained by
IMRPhenomXP is consistent with our understanding of the utility of higher
order modes. We emphasize that these findings are unrelated to waveform
systematics: we always inject and recover with the same waveform model.

On the population level, to reduce the complexity of the sampling, we conduct
analyses where we fit for only the spin magnitude or the tilt angle
distribution, while fixing the other to it’s true injected value. In theory,
this could help identify if one or the other of these parameters was the
driving factor for the mismatch between the true and recovered populations
seen in Fig. 2. The inferred spin magnitude distribution for the
LowSpinAligned population under the Beta+DoubleGaussian model with the tilt
distribution fixed to truth is shown in blue in the top panel of Fig. 12; the
bottom panel shows the inverse. These recoveries are indeed better than those
in Fig. 2 (plotted in navy dashed lines for comparison), i.e. the mean of the
$\chi$ distribution and mean of the larger sub-population of the $\cos\theta$
distribution are both more accurate. However, even in this much simplified
version, the Beta+DoubleGaussian model again fails to recover the truth, and
in particular still shows no signs of bimodality in the tilt distribution.

Figure 11:  Inferred distributions obtained with the Beta+DoubleGaussian model
for spin magnitude $\chi$, spin tilt $\cos\theta$, effective
$\chi_{\mathrm{eff}}$ spin, effective precessing spin $\chi_{\mathrm{p}}$, for
the three simulated populations with some form of reduced complexity when the
sampling of individual-event posteriors. The dashed lines show population
results with the IMRPhenomXP waveform (used for both injection and recovery),
which excludes higher order modes, for 70 (yellow) and 300 (blue) event
catalogs. The solid lines shown populations results from 100-event catalogs
using individual-event posteriors calculated with IMRPhenomXPHM, but with
various parameters fixed to their true values rather than sampled over. In
pink (“fixed extrinsic”), we have fixed the extrinsic parameters (i.e.
everything aside from masses and spins) to their true values. In orange
(“fixed extrinsic + masses”), we further restrict by only sampling over spin
magnitudes and tilts. Of these variations, only the “fixed extrinsic + masses”
individual-event posteriors yield population constraints that are improved
from those shown in Fig. 2.  Figure 12:  (Top) Inferred spin magnitude $\chi$
distribution (blue shaded) for 70 events from the LowSpinAligned population
under the Beta+DoubleGaussian population model with the tilt-angle
distribution fixed to truth (black line in the bottom subplot). (Bottom)
Inferred cosine of the tilt angle $\cos\theta$ distribution for the same 70
events and model with the spin magnitude distribution fixed to truth (black
line in the top subplot). In both subplots, 90% credible intervals for the
distributions inferred by fitting for both the spin magnitude and tilt
distributions simultaneously, as shown in Fig. 2, are shown in navy dashed
lines.

###  G.4 Hierarchical analysis with independent codes

It is always possible that our poor recovery of component spin distributions
is simply due to an unidentified error in the code used to perform
hierarchical inference. As a safeguard against this possibility, we have
repeated the hierarchical analysis of the HighSpinPrecessing, MediumSpin, and
LowSpinAligned populations using a second, distinct body of code. This
alternate analysis code was developed entirely independently, and furthermore
relies on a different stochastic sampler: whereas our main hierarchical
inference results are obtained using emcee, this alternative performs
inference using numpyro [56, 57], a probabilistic programming library
implemented with jax [81]. The numpyro-based code produces results nearly
identical to those obtained with our emcee-based code. This implies that our
results are not attributable to an unidentified error, unless that same error
was independently introduced into two bodies of code created by two different
analysts.

###  G.5 Simultaneously fitting the mass and redshift distributions

Yet another source of potential bias we investigate is the choice to fix,
rather than fit, the binary black hole primary mass and redshift
distributions. In principle, we do not expect significant covariance between
the inferred primary mass, redshift, and component spin distributions; our
simulated astrophysical populations have no underlying correlations between
mass, redshift, and spin. At the same time, inferred component spins are
expected to correlate strongly with the mass ratio distribution, which in turn
can depend systematically on the choice of primary mass distribution [35, 64].
Furthermore, it is known that assumptions regarding spin magnitudes can at
times affect inference of the high-redshift rate of black hole mergers [2, 3].
Given these possibilities, it is possible that fixing the presumed mass and
redshift distributions (even fixing them to the correct values, as we have
done) introduces bias into our spin measurements.

Figure 13:  As in Fig. 2, but now when additionally inferring the mass and
redshift distributions of the three simulated populations in conjunction with
their spin distributions. These results are furthermore produced with an
entirely independent body of code, using numpyro rather than emcee to
stochastically sample the population likelihood. Despite these differences,
the results are nearly identical to those in Fig. 2, indicating that the
difficult recovery of injected component spin distributions is due neither to
our choice to fix the mass and redshift distributions in the main text, nor to
unidentified errors in our hierarchical inference code.

To check this, we repeat our inference but now hierarchically inferring the
black hole mass and redshift distributions alongside the component spin
distributions. We model primary masses following the PowerLaw+Peak model [2,
3] and assume that the merger rate density evolves with redshift as
$(1+z)^{\kappa}$ for some parameter $\kappa$. This model also uses a slightly
different spin magnitude model: a truncated Normal distribution instead of a
Beta distribution. We perform this inference using the alternative numpyro-
based code introduced in App. G.4 above. The spin distributions inferred in
this case are shown in Fig. 13. The results are extremely similar to those in
Fig. 2. As before, we recover the HighSpinPrecessing and MediumSpin component
spin distributions reasonably well, but do not successfully measure the
LowSpinAligned distributions. In this latter case, we miss (or misplace) the
bimodality inherent in $\cos\theta$ and, accordingly, systematically
overestimate component spin magnitudes. Once more, though, the
$\chi_{\mathrm{eff}}$ distribution is well-recovered in all three cases.

###  G.6 Recovering rates and spins simultaneously

Figure 14:  As in Fig. 2, but now fitting for the absolute rate of black hole
mergers as a function of spin, rather than component spin probability
distributions. When producing these results, we additionally infer the black
hole mass and redshift distributions, as in Fig. 13 using Numpyro rather than
Emcee. The left-most column shows the differential merger rate
$dR/dm_{1}\,d\chi_{1}\,\chi_{2}$ per unit primary and secondary spin
magnitude, evaluated along the $\chi_{1}=\chi_{2}$ line. The second column
analogously shows the merger rate
$dR/dm_{1}\,d\cos\theta_{1}\,d\cos\theta_{2}$ across the
$\cos\theta_{1}=\cos\theta_{2}$ line, and the last two columns give the rates
and $dR/dm_{1}\,d\chi_{\mathrm{p}}$ as a function of effective spins. In all
cases, rates are evaluated at $m_{1}=30\,M_{\odot}$ and $z=0.2$. Even when
plotting rates, we draw the same qualitative conclusions as originally
identified in Fig. 2: the HighSpinPrecessing and MediumSpin populations are
recovered well, while the rate density of the LowSpinAligned population is
recovered poorly.

GW data are generated according to a Poisson point process, in which
individual compact binaries stochastically trace an underlying rate density
$dR/d\lambda$ of mergers across the space of binary parameters $\lambda$. When
performing hierarchical inference over GW catalogs, we are formally
reconstructing this rate density: measuring the “counts” of events occurring
in different regions of parameter space. Often, however, we are concerned only
about the shape of $dR/d\lambda$, not its normalization. In this case, it is
common to instead study the normalized probability distribution $p(\lambda)$
of source parameters. This is achieved after the fact by fitting for the
merger rate but only presenting $p(\lambda)$, or from the very outset by
marginalizing over and subsequently ignoring the total rate, a procedure that
gives Eq. (4).

While this procedure is usually well-behaved, there do exist instances in
which the choice to present the normalized $p(\lambda)$, rather than a
reconstructed rate density, can yield inadvertently misleading conclusions. In
some cases, models that successfully recover the correct rate density can
appear to fail in recovering the correct probability density; see discussion
in Sec. 5B of Callister and Farr [47]. Thus, when evaluating the goodness-of-
fit of a given model, the most robust results are obtained by comparing
injected and recovered merger rates, rather than injected and recovered
probability densities.

Given this discussion, does the poor agreement between injected and recovered
spin probability distributions signify a true modeling and inference failure?
Or is this disagreement illusory, due to our choice to compare probability
distributions rather than reconstructed merger rate densities? To check this,
we repeat the hierarchical analyses of the three simulated populations but now
fitting for the overall merger rate alongside the hyperparameters governing
the component spin distributions. As in Appendix G.5, we simultaneously infer
the primary mass and redshift distributions. Figure 14 shows our inferred
merger rates as a function of spin for each injected population. Our initial
conclusions hold: when simultaneously fitting for and presenting differential
merger rates, rather than probability densities, we still find that the
HighSpinPrecessing and MediumSpin populations are recovered well, but we do
not successfully recover the LowSpinAligned population. Hence our poor
recovery of LowSpinAligned is a real effect, rather than a bias or misleading
visualization related to our choice to marginalize over the absolute merger
rate.

###  G.7 Other miscellaneous checks for hierarchical inference

Finally, we present results from other miscellaneous verification methods for
our hierarchical inference procedure. First, we investigate different methods
of breaking the degeneracy between the two Gaussian components in tilt-
distribution portion of the Beta+DoubleGaussian model, see Eq. (17). For any
model defined as a mixture of multiple components, some method must be imposed
to break the degeneracy these components. For a bimodal Gaussian, this can be
done in three ways:

  1. 1.

Imposing an ordering of the means – assign “distribution 1” to be that with
the smaller mean and “distribution 2” to be that with the larger mean

  2. 2.

Imposing an ordering of the widths – assign “distribution 1” to be that which
is narrower, and “distribution 2” to be wider

  3. 3.

Limiting the mixing fraction be $\leq 0.5$ – assign “distribution 1” to be
that which contains a smaller fraction of events, and “distribution 2” to be
that which contains a larger fraction.

Sometimes one method of breaking the degeneracy converges better when used in
a hierarchical inference procedure than another. We find that this is not the
case in this work: different methods perform identically (within sampling
error), as seen in Fig. 15. Using the means (mixing fraction) of the Gaussians
to break the degeneracy yields the distributions plotted in blue (orange). The
results are consistent with each other. We opt to use the mixing fraction to
break degeneracy throughout the bulk of this work because it is more
computationally efficient.

Next, to ensure there is no misspecification in the selection function for
spins, see Eq. (8), we run hierarchical inference without any spin selection
effects. Results are shown for the LowSpinAligned population pink in Fig. 15
for 70 (dashed) and 300 (solid) events. Selection effects in component spins
are not strong, and thus are not expected to effect population inference
significantly.555Selections effects are strong for masses and redshift, on the
other hand. This can be seen when comparing the underlying and detected
distributions in Fig. 5. This is indeed the case, as the distributions
inferred without spin selection effects are nearly identical to this inferred
with them (shown in orange). We also note that for all results shown in this
work, the number of effective samples does not rail against the cut given in
Eq. (9).

Figure 15:  Inferred distributions (blue) for spin magnitude $\chi$ (top
subplot) and cosine of the tilt angle $\cos\theta$ (bottom subplot) for 70
events from the LowSpinAligned population obtained with the
Beta+DoubleGaussian population model under various conditions. All results
from 70 (300) event catalogs are plotted with dashed (solid) lines. In blue
and orange we compared two methods of breaking the degeneracy between the two
Gaussian sub-populations in the Beta+DoubleGaussian tilt distribution: first,
restricting the mean of distribution 1 to be smaller than that of distribution
2 ($\mu_{1,\cos\theta}<\mu_{2,\cos\theta}$; blue), and second, restricting the
mixing fraction to be less than one half, effectively assigning distribution 1
to be that containing less events ($f<0.5$; orange). In orange and pink, we
compare results where we do (orange) versus do not (pink) include spin
selection effects (both using the $f<0.5$ method of breaking degeneracy). None
of these variations produce population constraints improved from those shown
in Fig. 2; the bimodality of the tilt-angle distribution remains unrecovered.

Our final check is to conduct hierarchical inference on several different
random 70-event catalog realizations from the 300 total events per population.
Just like different Gaussian noise instantiations of the data lead to variance
in individual-event posteriors, so too can random catalog instantiation lead
to variance in the recovered posteriors on the population parameters. Some
catalogs will yield a more accurately recovered population than others, just
by random chance from working with finite numbers, see, e.g. Callister  _et
al._ [34]. Fig. 16 shows some expected variance in results, but nothing
corresponding to the degree of mismatch between the true and inferred tilt
distributions of the LowSpinAligned population seen in Fig. 2. We therefore
conclude that we cannot attribute bias between the injected and recovered
LowSpinAligned population to be from an “unlucky” catalog realization.
Additionally, each catalog instantiation leads to a different number of per-
event effective samples, see Eq. (10) and corresponding discussion in App. C,
which we find are not correlated to the by-eye goodness of fit, as seen in the
rightmost row of Fig. 16. This leads us to believe that the absence of
inferred bimodality is not due to the issue of our events not having enough
effective samples, although the fact that the minimum event-level
$N_{\mathrm{eff}}$ over catalog instantiations is small could be another
source of imperfect recovery.

Figure 16:  Inferred distributions obtained with the Beta+DoubleGaussian model
for spin magnitude $\chi$, spin tilt $\cos\theta$, effective
$\chi_{\mathrm{eff}}$ spin, effective precessing spin $\chi_{\mathrm{p}}$, for
various 70-event catalog instantiations of the three simulated populations,
each plotted in a different color. The true, underlying populations are shown
in black for comparison. While there is a small amount of variation, none of
the catalog instantiations yields population constraints improved from those
shown in Fig. 2. The rightmost column shows the minimum per-event
$N_{\mathrm{eff}}$ over the catalog used for inference for that particular
analysis (See Eq. (10) and corresponding discussion in App. C). The minumum
$N_{\mathrm{eff}}$ does not correlate with the by-eye goodness of fit of the
distributions in the first four columns.

## References

  * Blanchet [2014] L. Blanchet, Gravitational radiation from post-newtonian sources and inspiralling compact binaries, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration, & KAGRA Collaboration), Population properties of compact objects from the second ligo–virgo gravitational-wave transient catalog, Astrophys. J. Lett. 913, L7 (2021a).
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific Collaboration, Virgo Collaboration, & KAGRA Collaboration), The population of merging compact binaries inferred using gravitational waves through GWTC-3,  (2021b), arXiv:2111.03634 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, D. Aisa, N. Allemandou, _et al._ , Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32, 024001 (2015).
  * Gerosa and Berti [2017] D. Gerosa and E. Berti, Are merging black holes born from stellar collapse or previous mergers?, Phys. Rev. D 95, 124046 (2017), arXiv:1703.06223 [gr-qc] .
  * Rodriguez  _et al._ [2019] C. L. Rodriguez, M. Zevin, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, and C. S. Ye, Black holes: The next generation—repeated mergers in dense star clusters and their gravitational-wave properties, Phys. Rev. D 100, 043027 (2019), arXiv:1906.10260 [astro-ph.HE] .
  * Bavera  _et al._ [2020] S. S. Bavera, T. Fragos, Y. Qin, E. Zapartas, C. J. Neijssel, I. Mandel, A. Batta, S. M. Gaebel, C. Kimball, and S. Stevenson, The origin of spin in binary black holes: Predicting the distributions of the main observables of Advanced LIGO, Astron. Astrophys. 635, A97 (2020), arXiv:1906.12257 [astro-ph.HE] .
  * Vitale  _et al._ [2017a] S. Vitale, R. Lynch, R. Sturani, and P. Graff, Use of gravitational waves to probe the formation channels of compact binaries, Class. Quant. Grav. 34, 03LT01 (2017a), arXiv:1503.04307 [gr-qc] .
  * Rodriguez  _et al._ [2016] C. L. Rodriguez, M. Zevin, C. Pankow, V. Kalogera, and F. A. Rasio, Illuminating Black Hole Binary Formation Channels with Spins in Advanced LIGO, Astrophys. J. Lett. 832, L2 (2016), arXiv:1609.05916 [astro-ph.HE] .
  * Farr  _et al._ [2017] W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, Distinguishing Spin-Aligned and Isotropic Black Hole Populations With Gravitational Waves, Nature 548, 426 (2017), arXiv:1706.01385 [astro-ph.HE] .
  * Gerosa  _et al._ [2018] D. Gerosa, E. Berti, R. O’Shaughnessy, K. Belczynski, M. Kesden, D. Wysocki, and W. Gladysz, Spin orientations of merging black holes formed from the evolution of stellar binaries, Phys. Rev. D 98, 084036 (2018), arXiv:1808.02491 [astro-ph.HE] .
  * Mandel and Farmer [2022] I. Mandel and A. Farmer, Merging stellar-mass binary black holes, Phys. Rept. 955, 1 (2022), arXiv:1806.05820 [astro-ph.HE] .
  * Zevin  _et al._ [2021] M. Zevin, S. S. Bavera, C. P. L. Berry, V. Kalogera, T. Fragos, P. Marchant, C. L. Rodriguez, F. Antonini, D. E. Holz, and C. Pankow, One Channel to Rule Them All? Constraining the Origins of Binary Black Holes Using Multiple Formation Pathways, Astrophys. J. 910, 152 (2021), arXiv:2011.10057 [astro-ph.HE] .
  * Fuller and Ma [2019] J. Fuller and L. Ma, Most Black Holes are Born Very Slowly Rotating, Astrophys. J. Lett. 881, L1 (2019), arXiv:1907.03714 [astro-ph.SR] .
  * Qin  _et al._ [2018] Y. Qin, T. Fragos, G. Meynet, J. Andrews, M. Sørensen, and H. F. Song, The spin of the second-born black hole in coalescing binary black holes, Astron. Astrophys. 616, A28 (2018), arXiv:1802.05738 [astro-ph.SR] .
  * Bavera  _et al._ [2021] S. S. Bavera _et al._ , The impact of mass-transfer physics on the observable properties of field binary black hole populations, Astron. Astrophys. 647, A153 (2021), arXiv:2010.16333 [astro-ph.HE] .
  * Steinle and Kesden [2021] N. Steinle and M. Kesden, Pathways for producing binary black holes with large misaligned spins in the isolated formation channel, Phys. Rev. D 103, 063032 (2021), arXiv:2010.00078 [astro-ph.HE] .
  * Zevin and Bavera [2022] M. Zevin and S. S. Bavera, Suspicious Siblings: The Distribution of Mass and Spin Across Component Black Holes in Isolated Binary Evolution 10.48550/arXiv.2203.02515 (2022), arXiv:2203.02515 [astro-ph.HE] .
  * Pürrer  _et al._ [2016] M. Pürrer, M. Hannam, and F. Ohme, Can we measure individual black-hole spins from gravitational-wave observations?, Phys. Rev. D 93, 084042 (2016), arXiv:1512.04955 [gr-qc] .
  * Vitale  _et al._ [2017b] S. Vitale, R. Lynch, V. Raymond, R. Sturani, J. Veitch, and P. Graff, Parameter estimation for heavy binary-black holes with networks of second-generation gravitational-wave detectors, Phys. Rev. D 95, 064053 (2017b), arXiv:1611.01122 [gr-qc] .
  * Chatziioannou  _et al._ [2018] K. Chatziioannou, G. Lovelace, M. Boyle, M. Giesler, D. A. Hemberger, R. Katebi, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, and B. Szilágyi, Measuring the properties of nearly extremal black holes with gravitational waves, Phys. Rev. D 98, 044028 (2018), arXiv:1804.03704 [gr-qc] .
  * Biscoveanu  _et al._ [2021] S. Biscoveanu, M. Isi, V. Varma, and S. Vitale, Measuring the spins of heavy binary black holes, Phys. Rev. D 104, 103018 (2021), arXiv:2106.06492 [gr-qc] .
  * Miller  _et al._ [2023] S. J. Miller, M. Isi, K. Chatziioannou, V. Varma, and I. Mandel, GW190521: tracing imprints of spin-precession on the most massive black hole binary, arXiv:2310.01544 [astro-ph.HE] (2023).
  * Khan  _et al._ [2020] S. Khan, F. Ohme, K. Chatziioannou, and M. Hannam, Including higher order multipoles in gravitational-wave models for precessing binary black holes, Phys. Rev. D 101, 024056 (2020), arXiv:1911.06050 [gr-qc] .
  * Varma  _et al._ [2019] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer, Surrogate models for precessing binary black hole simulations with unequal masses, Phys. Rev. Research. 1, 033015 (2019), arXiv:1905.09300 [gr-qc] .
  * Pratten  _et al._ [2021] G. Pratten _et al._ , Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Ossokine  _et al._ [2020] S. Ossokine _et al._ , Multipolar Effective-One-Body Waveforms for Precessing Binary Black Holes: Construction and Validation, Phys. Rev. D 102, 044055 (2020), arXiv:2004.09442 [gr-qc] .
  * Racine [2008] E. Racine, Analysis of spin precession in binary black hole systems including quadrupole-monopole interaction, Phys. Rev. D 78, 044021 (2008), arXiv:0803.1820 [gr-qc] .
  * Schmidt  _et al._ [2015] P. Schmidt, F. Ohme, and M. Hannam, Towards models of gravitational waveforms from generic binaries II: Modelling precession effects with a single effective precession parameter, Phys. Rev. D 91, 024043 (2015), arXiv:1408.1810 [gr-qc] .
  * Apostolatos  _et al._ [1994] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Spin-induced orbital precession and its modulation of the gravitational waveforms from merging binaries, Phys. Rev. D 49, 6274 (1994).
  * Miller  _et al._ [2020] S. Miller, T. A. Callister, and W. Farr, The Low Effective Spin of Binary Black Holes and Implications for Individual Gravitational-Wave Events, Astrophys. J. 895, 128 (2020), arXiv:2001.06051 [astro-ph.HE] .
  * Roulet  _et al._ [2021] J. Roulet, H. S. Chia, S. Olsen, L. Dai, T. Venumadhav, B. Zackay, and M. Zaldarriaga, Distribution of effective spins and masses of binary black holes from the LIGO and Virgo O1–O3a observing runs, Phys. Rev. D 104, 083010 (2021), arXiv:2105.10580 [astro-ph.HE] .
  * Callister  _et al._ [2022] T. A. Callister, S. J. Miller, K. Chatziioannou, and W. M. Farr, No Evidence that the Majority of Black Holes in Binaries Have Zero Spin, Astrophys. J. Lett. 937, L13 (2022), arXiv:2205.08574 [astro-ph.HE] .
  * Ng  _et al._ [2018] K. K. Y. Ng, S. Vitale, A. Zimmerman, K. Chatziioannou, D. Gerosa, and C.-J. Haster, Gravitational-wave astrophysics with effective-spin measurements: asymmetries and selection biases, Phys. Rev. D 98, 083007 (2018), arXiv:1805.03046 [gr-qc] .
  * Chatziioannou and Farr [2020] K. Chatziioannou and W. M. Farr, Inferring the maximum and minimum mass of merging neutron stars with gravitational waves, Phys. Rev. D 102, 064063 (2020), arXiv:2005.00482 [astro-ph.HE] .
  * Vitale  _et al._ [2022] S. Vitale, S. Biscoveanu, and C. Talbot, Spin it as you like: The (lack of a) measurement of the spin tilt distribution with LIGO-Virgo-KAGRA binary black holes, Astron. Astrophys. 668, L2 (2022), arXiv:2209.06978 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific Collaboration & Virgo Collaboration), Binary black hole population properties inferred from the first and second observing runs of advanced LIGO and advanced virgo, Astrophys. J. 882, L24 (2019).
  * Golomb and Talbot [2022] J. Golomb and C. Talbot, Searching for structure in the binary black hole spin distribution, (2022), arXiv:2210.12287 [astro-ph.HE] .
  * Tong  _et al._ [2022] H. Tong, S. Galaudage, and E. Thrane, Population properties of spinning black holes using the gravitational-wave transient catalog 3, Phys. Rev. D 106, 103019 (2022), arXiv:2209.02206 [astro-ph.HE] .
  * Edelman  _et al._ [2022] B. Edelman, B. Farr, and Z. Doctor, Cover Your Basis: Comprehensive Data-Driven Characterization of the Binary Black Hole Population, (2022), arXiv:2210.12834 [astro-ph.HE] .
  * Poisson [2004] E. Poisson, Absorption of mass and angular momentum by a black hole: Time-domain formalisms for gravitational perturbations, and the small-hole / slow-motion approximation, Phys. Rev. D 70, 084044 (2004), arXiv:gr-qc/0407050 .
  * Gerosa  _et al._ [2015] D. Gerosa, M. Kesden, U. Sperhake, E. Berti, and R. O’Shaughnessy, Multi-timescale analysis of phase transitions in precessing black-hole binaries, Phys. Rev. D 92, 064016 (2015), arXiv:1506.03492 [gr-qc] .
  * Gerosa  _et al._ [2021] D. Gerosa, M. Mould, D. Gangardt, P. Schmidt, G. Pratten, and L. M. Thomas, A generalized precession parameter $\chi_{\mathrm{p}}$ to interpret gravitational-wave data, Phys. Rev. D 103, 064067 (2021), arXiv:2011.11948 [gr-qc] .
  * Thomas  _et al._ [2021] L. M. Thomas, P. Schmidt, and G. Pratten, New effective precession spin for modeling multimodal gravitational waveforms in the strong-field regime, Phys. Rev. D 103, 083022 (2021), arXiv:2012.02209 [gr-qc] .
  * Schmidt  _et al._ [2012] P. Schmidt, M. Hannam, and S. Husa, Towards models of gravitational waveforms from generic binaries: A simple approximate mapping between precessing and non-precessing inspiral signals, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Callister and Farr [2023] T. A. Callister and W. M. Farr, A Parameter-Free Tour of the Binary Black Hole Population, (2023), arXiv:2302.07289 [astro-ph.HE] .
  * Collaboration and Collaboration [2020] L. S. Collaboration and V. Collaboration, Noise curves used for simulations in the update of the observing scenarios paper (2020).
  * Speagle [2020] J. S. Speagle, DYNESTY: a dynamic nested sampling package for estimating Bayesian posteriors and evidences, MNRAS 493, 3132 (2020), arXiv:1904.02180 [astro-ph.IM] .
  * Ashton  _et al._ [2019] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, _et al._ , Bilby: A user-friendly bayesian inference library for gravitational-wave astronomy, The Astrophysical Journal Supplement Series 241, 27 (2019).
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw _et al._ , Bayesian inference for compact binary coalescences with bilby: validation and application to the first LIGO–Virgo gravitational-wave transient catalogue, Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Fishbach  _et al._ [2020] M. Fishbach, W. M. Farr, and D. E. Holz, The Most Massive Binary Black Hole Detections and the Identification of Population Outliers, Astrophys. J. Lett. 891, L31 (2020), arXiv:1911.05882 [astro-ph.HE] .
  * Fairhurst  _et al._ [2023] S. Fairhurst, C. Hoy, R. Green, C. Mills, and S. A. Usman, Simple parameter estimation using observable features of gravitational-wave signals, Phys. Rev. D 108, 082006 (2023), arXiv:2304.03731 [gr-qc] .
  * Farah  _et al._ [2023] A. M. Farah, B. Edelman, M. Zevin, M. Fishbach, J. M. Ezquiaga, B. Farr, and D. E. Holz, Things That Might Go Bump in the Night: Assessing Structure in the Binary Black Hole Mass Spectrum, Astrophys. J. 955, 107 (2023), arXiv:2301.00834 [astro-ph.HE] .
  * Foreman-Mackey  _et al._ [2013] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, Publications of the Astronomical Society of the Pacific 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Phan  _et al._ [2019] D. Phan, N. Pradhan, and M. Jankowiak, Composable effects for flexible and accelerated probabilistic programming in numpyro, arXiv:1912.11554 (2019).
  * Bingham  _et al._ [2019] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, _et al._ , Pyro: Deep universal probabilistic programming, J. Mach. Learn. Res. 20, 28:1 (2019).
  * Loredo [2004] T. J. Loredo, Accounting for source uncertainties in analyses of astronomical survey data, AIP Conf. Proc. 735, 195 (2004), arXiv:astro-ph/0409387 .
  * Mandel  _et al._ [2019] I. Mandel, W. M. Farr, and J. R. Gair, Extracting distribution parameters from multiple uncertain observations with selection biases, Mon. Not. Roy. Astron. Soc. 486, 1086 (2019), arXiv:1809.02063 [physics.data-an] .
  * Vitale  _et al._ [2020] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, Inferring the properties of a population of compact binaries in presence of selection effects, arXiv:2007.05579 [astro-ph.IM] (2020).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On information and sufficiency, Annals of Mathematical Statistics 22, 79 (1951).
  * Samantha R Cook and Rubin [2006] A. G. Samantha R Cook and D. B. Rubin, Validation of software for bayesian models using posterior quantiles, Journal of Computational and Graphical Statistics 15, 675 (2006).
  * Talts  _et al._ [2018] S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman, Validating Bayesian Inference Algorithms with Simulation-Based Calibration, arXiv e-prints , arXiv:1804.06788 (2018), arXiv:1804.06788 [stat.ME] .
  * Biscoveanu  _et al._ [2022] S. Biscoveanu, C. Talbot, and S. Vitale, The effect of spin mismodelling on gravitational-wave measurements of the binary neutron star mass distribution, Mon. Not. Roy. Astron. Soc. 511, 4350 (2022), arXiv:2111.13619 [astro-ph.HE] .
  * Fishbach and Holz [2020] M. Fishbach and D. E. Holz, Minding the gap: GW190521 as a straddling binary, Astrophys. J. Lett. 904, L26 (2020), arXiv:2009.05472 [astro-ph.HE] .
  * Essick and Fishbach [2023] R. Essick and M. Fishbach, DAGnabbit! Ensuring Consistency between Noise and Detection in Hierarchical Bayesian Inference, (2023), arXiv:2310.02017 [gr-qc] .
  * Sinharay and Stern [2003] S. Sinharay and H. S. Stern, Posterior predictive model checking in hierarchical models, Journal of Statistical Planning and Inference 111, 209 (2003).
  * Bayarri and Castellanos [2007] M. J. Bayarri and M. E. Castellanos, Bayesian checking of the second levels of hierarchical models, Statistical Science 22, 322 (2007).
  * Callister [2021] T. Callister, Reweighting single event posteriors with hyperparameter marginalization (2021).
  * Talbot and Thrane [2017] C. Talbot and E. Thrane, Determining the population properties of spinning black holes, Phys. Rev. D 96, 023012 (2017), arXiv:1704.08370 [astro-ph.HE] .
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW150914: First results from the search for binary black hole coalescence with Advanced LIGO, Phys. Rev. D 93, 122003 (2016), arXiv:1602.03839 [gr-qc] .
  * Hannam  _et al._ [2014] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Smith  _et al._ [2020] R. J. E. Smith, C. Talbot, F. Hernandez Vivanco, and E. Thrane, Inferring the population properties of binary black holes from unresolved gravitational waves, Mon. Not. Roy. Astron. Soc. 496, 3281 (2020), arXiv:2004.09700 [astro-ph.HE] .
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357 (2020).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, _et al._ , SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods 17, 261 (2020).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2d graphics environment, Computing In Science & Engineering 9, 90 (2007).
  * Waskom [2021] M. L. Waskom, seaborn: statistical data visualization, Journal of Open Source Software 6, 3021 (2021).
  * Astropy Collaboration [2013] Astropy Collaboration, Astropy: A community Python package for astronomy, Astronomy & Astrophysics 558, A33 (2013), arXiv:1307.6212 [astro-ph.IM] .
  * Astropy Collaboration [2018] Astropy Collaboration, The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package, Astronomical J. 156, 123 (2018), arXiv:1801.02634 [astro-ph.IM] .
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, _et al._ , JAX: composable transformations of Python+NumPy programs (2018).
  * Talbot and Thrane [2018] C. Talbot and E. Thrane, Measuring the binary black hole mass spectrum with an astrophysically motivated parameterization, Astrophys. J. 856, 173 (2018), arXiv:1801.02699 [astro-ph.HE] .
  * Ade  _et al._ [2014] P. A. R. Ade _et al._ (Planck), Planck 2013 results. XVI. Cosmological parameters, Astron. Astrophys. 571, A16 (2014), arXiv:1303.5076 [astro-ph.CO] .
  * Thorne [1987] K. S. Thorne, Gravitational radiation., in  _Three Hundred Years of Gravitation_ (1987) pp. 330–458.
  * Finn and Chernoff [1993] L. S. Finn and D. F. Chernoff, Observing binary inspiral in gravitational radiation: One interferometer, Phys. Rev. D 47, 2198 (1993), arXiv:gr-qc/9301003 .
  * Allen  _et al._ [2012] B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. E. Creighton, FINDCHIRP: An Algorithm for detection of gravitational waves from inspiraling compact binaries, Phys. Rev. D 85, 122006 (2012), arXiv:gr-qc/0509116 .
  * Essick [2023] R. Essick, Semianalytic sensitivity estimates for catalogs of gravitational-wave transients, Phys. Rev. D 108, 043011 (2023), arXiv:2307.02765 [gr-qc] .
  * Peters [1964] P. C. Peters, Gravitational radiation and the motion of two point masses, Phys. Rev. 136, B1224 (1964).
  * Blanchet  _et al._ [1995] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will, and A. G. Wiseman, Gravitational-radiation damping of compact binary systems to second post-newtonian order, Phys. Rev. Lett. 74, 3515 (1995).
  * Thrane and Talbot [2019] E. Thrane and C. Talbot, An introduction to bayesian inference in gravitational-wave astronomy: Parameter estimation, model selection, and hierarchical models, Publications of the Astronomical Society of Australia 36, e010 (2019).
  * Farr [2019] W. M. Farr, Accuracy Requirements for Empirically Measured Selection Functions, Research Notes of the American Astronomical Society 3, 66 (2019), arXiv:1904.10879 [astro-ph.IM] .
  * Talbot and Golomb [2023] C. Talbot and J. Golomb, Growing pains: understanding the impact of likelihood uncertainty on hierarchical Bayesian inference for gravitational-wave astronomy, Mon. Not. Roy. Astron. Soc. 526, 3495 (2023), arXiv:2304.06138 [astro-ph.IM] .
  * Callister  _et al._ [2021] T. A. Callister, C.-J. Haster, K. K. Y. Ng, S. Vitale, and W. M. Farr, Who Ordered That? Unequal-mass Binary Black Hole Mergers Have Larger Effective Spins, Astrophys. J. Lett. 922, L5 (2021), arXiv:2106.00521 [astro-ph.HE] .