Evidence for a correlation between binary black hole mass ratio and black-hole spins

The astrophysical origins of the binary black hole systems seen with gravitational waves are still not well understood. However, features in the distribution of black-hole masses, spins, redshifts, and eccentricities provide clues into how these systems form. Much has been learned by investigating these distributions one parameter at a time. However, we can extract additional information by studying the covariance between pairs of parameters. Previous work has shown preliminary support for an anti-correlation between mass ratio $q \equiv m_2/m_1$ and effective inspiral spin $\chi_\mathrm{eff}$ in the binary black hole population. In this study, we test for the existence of this anti-correlation using updated data from the third gravitational wave transient catalogue (GWTC-3) and improve our copula-based framework to employ a more robust model for black-hole spins. We find evidence for an anti-correlation in $(q, \chi_\mathrm{eff})$ with 99.7% credibility. This may imply high common-envelope efficiencies, stages of super-Eddington accretion, or a tendency for binary black hole systems to undergo mass-ratio reversal during isolated evolution. Covariance in $(q,\chi_\mathrm{eff})$ may also be used to investigate the physics of tidal spin-up as well as the properties of binary-black-hole-forming active galactic nuclei.


INTRODUCTION
There are a number of theorised formation mechanisms for merging binary black hole (BBH) systems (see reviews by Mapelli 2018;Mandel & Farmer 2022;Spera et al. 2022).These scenarios typically fall into one of two broad categories.The first is dynamical formation, in which black holes are brought together through interactions with other black holes to form merging binaries in dense stellar environments.The second is isolated (or "field") formation, in which stellar black-hole progenitors already in a binary evolve into black holes together and precede to merge afterwards.However, it remains unclear the relative proportions these channels contribute to the BBH population observed in gravitational waves (see Abbott et al. 2021a) by the LIGO-Virgo-KAGRA collaboration (LVK; Aasi et al. 2015;Acernese et al. 2015;Akutsu et al. 2021).
In our previous study (Adamcewicz & Thrane 2022), we suggested that investigations like the one carried out by Callister et al. (2021) (and consequently Abbott et al. 2023) are in-principle susceptible to a subtle statistical effect: if the marginal distribution of one parameter implicitly varies with the correlation parameter, then it is difficult to know if the apparent correlation is genuinely present in the data, or if it only seems so because of an improved fit of the marginal distribution.In other words, one may ask: does LVK data support an anticorrelation between (q, χ eff ), or does the data merely prefer the marginal distribution for χ eff implied by anticorrelation?Adamcewicz & Thrane (2022) addresses this concern using a copula density function (Sklar 1996): a statistical tool that imposes a variable correlation κ between two parameters but keeps their corresponding marginal distributions fixed.Adamcewicz & Thrane (2022) corroborates the results of Callister et al. (2021), inferring an anti-correlation in (q, χ eff ) with 98.7% credibility with GWTC-2 events and demonstrating that the evidence for anti-correlation is not due merely to a better model for χ eff . 23 However, Adamcewicz & Thrane (2022) does not provide an extensive investigation of model misspecification to determine if the anti-correlation can be explained away as arising due to some assumption in our models. 4Callister et al. (2021) and Adamcewicz & Thrane (2022) mention the potential for the inferred anti-correlation to be a result of two separate, independent sub-populations, neither of which exhibit the observed anti-correlation, but which produce a global anticorrelation via the amalgamation paradox (see Simpson 1951;Pearson et al. 1899;Yule 1903).
Perhaps a more pressing concern, however, is that the measurement could be affected by an unmodeled subpopulation of non-spinning χ eff = 0 systems.Given the inconclusive evidence for a non-spinning sub-population of BBH systems found in the literature (Tong et al. 2022;Callister et al. 2022;Mould et al. 2022;Galaudage et al. 2021;Roulet et al. 2021), it is important to consider how such a feature may influence the inferred shape of the distribution of mass ratio and χ eff .
In this work, we aim to more thoroughly investigate the purported (q, χ eff ) anti-correlation from Callister et al. (2021), Adamcewicz &Thrane (2022), andAbbott et al. (2023).We begin in Section 2 by constructing a more comprehensive population model than the one used in Adamcewicz & Thrane (2022), which allows for a non-spinning sub-population.We then fit this model to BBH data from the third gravitational wave transient catalogue (GWTC-3;Abbott et al. 2021a) and present the inferred (q, χ eff ) correlation, as well as other population properties, in Section 4. We conclude in Section 5, exploring the implications of these results.

Effective coordinates versus physical coordinates
2 By coincidence, the analyses in Callister et al. (2021) and Adamcewicz & Thrane (2022) result in the same 98.7% credibility for anti-correlation, despite utilising different statistical frameworks with different assumptions. 3Adamcewicz & Thrane (2022) also shows that the anti-correlation between (q, χ eff ) is not an inadvertent consequence of implicit assumptions about the distribution of physical spin parameters {χ 1,2 , cos t 1,2 } (see, for example Abbott et al. 2023). 4For an overview of misspecification in gravitational-wave astronomy, see Romero-Shaw et al. (2022).
Before we define our model for binary black hole spin, it is useful to distinguish between two different coordinate systems.Binary black hole spins are defined by six parameters, three for each black hole: the spin magnitudes of each black hole χ 1,2 , (the cosines of) the spin tilts of each black hole cos t 1,2 , and the the azimuthal angles of each component ϕ 1,2 .While the distribution of spin magnitudes and tilts is believed to vary depending on the formation channel, the azimuthal angles are very likely to be uniformly distributed no matter the formation scenario.5Hence, models for BBH spin are fully characterised by four dimensions: {χ 1,2 , cos t 1,2 }.
Both Callister et al. (2021) and Adamcewicz & Thrane (2022) model black hole spin using only the effective inspiral spin χ eff .The advantage of using this parameter is that χ eff is an approximate constant of motion (Damour 2001) and is relatively well measured from gravitationalwave signals.In contrast, the spin tilts cos t 1,2 evolve as the binary precesses, and are relatively poorly measured.However, modelling only the distribution of χ eff does not fully determine the distribution of binary black hole spins.The distributions of the remaining three spin dimensions are defined implicitly based on the priors used for initial parameter estimation (see Thrane & Talbot 2019), which may not be physically well-motivated.
To address this, we explicitly model χ eff as well as three additional spin parameters so that the distribution of black-hole spins is fully specified.Following Roulet et al. (2021), we define as well as the in-plane spin magnitudes of each black hole where i = {1, 2}.We refer to the set of {χ eff , χ diff , ρ 1,2 } as the effective coordinate system, which can be contrasted with the physical coordinates {χ 1,2 , cos t 1,2 }.Either set can be used to fully specify the distribution of BBH spin.

Effective spin coordinate model
In this subsection we describe a phenomenological model for the marginal distributions of the four effective spin coordinates.Our population model is inspired by the so-called Extended model recently used in Galaudage et al. (2021) and Tong et al. (2022), which has the flexibility for two sub-populations: one with spinning black holes, and one with non-spinning black holes.For the moment, we set aside the (trivially defined) non-spinning sub-population (which we revisit in Section 2.4), and describe our model for spinning black holes.Since the Extended model employs physical coordinates, our first step is to define an effectivecoordinate model that is phenomenologically similar.The distributions of effective spin parameters implied by the Extended model are shown in Fig. 1 (blue).We parameterise the distributions of χ eff , χ diff and ρ 1,2 , with the intent of capturing the salient features from the Extended model.Our chosen effective spin models are over-plotted in Fig. 1 (orange).A summary of the models is provided in Table 1.
We model χ eff and χ diff using a skewed normal distribution truncated on the interval [−1, 1]: where Λ denotes the set of all hyper-parameters controlling the shape of the population distributions.The (truncated) skewed normal distribution is Here, N (x|ξ, ω) is a normal distribution with mean ξ and width ω.Meanwhile, a controls the skew, and Θ is a Heaviside step function.The mean of N is not the mean of N s and so one should think of ξ, ω, and a as shape parameters.In Appendix A, we derive the mean µ, variance σ 2 , and skewness γ of the N s distribution in terms of these three shape parameters.This allows us to rewrite our population models for χ eff and χ diff in their final form: We model ρ 1,2 using beta distributions: where µ ρ is the Beta distribution mean and σ 2 ρ is the variance.Comparing the blue and orange intervals in Fig. 1, we observe that our new effective-coordinate parameterisation produces similar distributions to the ones implied by the Extended model.

Correlating q and χ eff
The next step is to combine our effective spin coordinate model with a model for black-hole masses.We take In this demonstration, the hyper-parameters for our model are chosen to resemble the Extended model distributions in blue.The correlation parameter κ is set to zero in this example.Varying κ changes the degree of correlation in the (q, χ eff ) plane, but leaves every other panel unchanged.The non-spinning sub-population is also excluded from this demonstration.
primary mass m 1 and mass ratio q to be distributed according to the Power-Law + Peak mass model from Talbot & Thrane (2018).For completeness, we define this model in Appendix B. With all of the marginal distributions defined, we employ a Frank copula density function (Frank 1979) in order to introduce a correlation between q and χ eff : e κ − e κu − e κv + e κ(u+v) 2 . (11) Here, u and v are copula coordinates determined by the marginal distributions of q and χ eff respectively: For a Frank copula, κ < 0 implies correlation, κ > 0 implies anti-correlation, and κ = 0 implies no correlation.
For a more detailed introduction to copulas in population studies, see Adamcewicz & Thrane (2022).

Final touches
For spinning black holes, the joint distribution of mass and spin is where we denote the full set of parameters θ = {m 1 , q, χ eff , χ diff , ρ 1,2 }.Here, π(m 1 , q|Λ) is the Power-Law + Peak mass model described in Appendix B, and π(χ eff , χ diff , ρ 1,2 |Λ) is the product of the onedimensional spin distributions defined in Section 2.2.Note that while the distribution of effective spin coordinates factorizes, the distribution of physical spin coordinates does not.
Our model includes a second sub-population of nonspinning black holes, described as where δ is a Dirac delta function, and π(m 1 , q|Λ) is the same Power-Law + Peak mass model used for the spinning sub-population.Putting everything together, the full model is Here, λ 0 is a mixing fraction corresponding to the fraction of binaries with non-spinning black holes.The priors for each hyper-parameter in Λ, along with κ and λ 0 , are listed in Table 1.

METHOD
We measure the population hyper-parameters using hierarchical Bayesian inference; see Appendix C and Thrane & Talbot (2019).We employ the population inference package GWPopulation (Talbot et al. 2019) -which is built on top of Bilby (Ashton et al. 2019;Romero-Shaw et al. 2020) -and employ the nested sampler DYNSETY (Speagle 2020).We use BBH events from GWTC-3 (Abbott et al. 2021a;LVK 2021LVK , 2023) ) that were considered reliable for inclusion in the LVK O3b population analysis (Abbott et al. 2023).However, following Tong et al. (2022), we exclude the events GW200129 065458 and GW191109 010717 (hereby shortened to GW200129 and GW191109 respectively) from our analysis over data-quality concerns.The former event GW200129 was coincident with a detector glitch.While parameter estimation on cleaned data shows strong evidence for spin-precession (Abbott et al. 2021a;Hannam et al. 2022), the results may be sensitive to the glitch subtraction procedure (Payne et al. 2022), which is still under discussion (Davis et al. 2022;Macas & Lundgren 2023).The latter event GW191109 is still under investigation (Tong 2023).We repeat our analysis both with and without GW191109, bringing the total number of BBH events surveyed to 68 and 67 in each case.For more details on how parameter estimation is utilised, how our population likelihood is constructed, and how we account for selection biases in our analyses see Appendix C.
To correct for this, unphysical samples should be discarded and the population distribution re-normalised accordingly.We carry out this step in post-processingcomputing a weight for each hyper-posterior sample according to the appropriate normalisation factor.Before this correction, values of χ 1,2 > 1 account for ≲ 1% of the population for any given draw of Λ and κ from the hyper-posterior.As such, the effect of this normalisation is small.
The correlation hypothesis is preferred over the nocorrelation hypothesis with a natural log Bayes factor of ln BF = 2.2 excluding GW191109 or ln BF = 3.2 including GW191109.Moreover, the correlation hypothesis yields a maximum natural log likelihood that is ∆ ln L max = 6.4 (9.5) higher than the no-correlation hypothesis including (excluding) GW191109, suggesting that the data are better fit by the correlation model.The Bayes factor is smaller than the ratio of maximum likelihoods because the correlation hypothesis incurs an Occam's factor as a result of the additional correlation parameter κ.Thus, we report increased evidence for anti-correlation between (q, χ eff ) from Callister et al. (2021) and Adamcewicz & Thrane (2022), both of which analysed only the 44 BBH events from GWTC-2 (Abbott et al. 2021c), and from Abbott et al. (2023), which analysed all events from GWTC-3 (Abbott et al. 2021a).
The anti-correlation in (q, χ eff ) is illustrated by the reconstructed population distribution (population predictive distribution), plotted in Fig. 3. Superimposed over this distribution, is a set of weighted delta functions π,  2023), consisting of a set of delta functions that attempt to maximise the population likelihood.Darker tones indicate a higher weight, thus higher importance, for the given feature.Our model matches qualitative features present in the π result.Note that samples drawn from the non-spinning sub-population have been randomly off-set by χ eff ∼ ±0.005 in order to help with the visualisation.
which attempt to maximise the population likelihood, taken from Payne & Thrane (2023).Essentially, π is a data-driven population model, in which the constituent delta functions latch onto the most important features in the data, indicating where most of the population weight should lie (see Payne & Thrane 2023, for a more detailed explanation).Our phenomenological model reproduces the most highly-weighted features.
In Appendix D, we explore the impact that changing our choice of copula density function has on our measurement of a correlation κ.Regardless of the copula we choose, we find that the anti-correlation in (q, χ eff ) is still strongly favoured.
Finally, in Appendix E, we share results peripheral to the goal of this work, including posterior distributions on all hyper-parameters Λ, and full reconstructed population distributions.

DISCUSSION
Although we have shown stronger evidence for an anticorrelation in (q, χ eff ), a modeler's job is never done: while we have addressed the most obvious sources of model-dependency, additional work should be carried out to investigate other possible sources of misspecification.
The possibility of the observed anti-correlation being a result of the amalgamation paradox (see the discussion in Section 1) remains a point of interest.Baibhav et al. (2023), for instance, show that when modelling field-like and dynamical-like sub-populations in physical spin coordinates, the two sub-populations prefer different distributions of mass ratio, which could drive such an effect.Although Baibhav et al. (2023) emphasise that we are unable to distinguish between multiple sub-populations in the marginal distribution of χ eff (a point further exemplified in non-parametric population studies; Callister & Farr 2023;Edelman et al. 2023), it remains possible that more clearly separated sub-populations could be seen in the q-χ eff plane.
This point aside, there are multiple astrophysical explanations for the apparent (q, χ eff ) anti-correlation that also merit future investigation.Simulations of massratio reversed field binaries by Broekgaarden et al. (2022) yield a distribution of (q, χ eff ), qualitatively similar to our reconstructed population (Fig. 3).We speculate that the observed anti-correlation could be an imprint of mass-ratio reversal in field binaries.
On the other hand, Bavera et al. (2021) and Zevin & Bavera (2022) predict a similar anti-correlation in isolated BBH systems that undergo common-envelope evolution with high common-envelope efficiency or stable mass transfer involving super-Eddington accretion.These scenarios, however, would not be expected to dominate the observed BBH population.Ma & Fuller (2023) predict that tidal spin-up rate (see also Qin et al. 2018;Hu et al. 2022) should scale with q 2 .Given that isolated systems have preferentially aligned black hole spins (positive χ eff ), such a scaling would presumably result in a positive correlation in (q, χ eff ).This seems to be in tension with our results.However, Ma & Fuller (2023) also cite a number of caveats that could complicate the relationship between mass ratio and tidal spin-up, which are yet to be explored.
An anti-correlated (q, χ eff ) might arise from some dynamical channels.For example, McKernan et al. (2022) show how such a feature may arise due to BBH formation in active galactic nuclei (AGN), provided a number of assumptions about the properties of the AGN in question hold.Namely, that black holes in the disk of the AGN are both more massive and have orbitally aligned spins, while those outside are lighter and isotropically distributed in spin orientation.Black holes outside the disk are then assumed to undergo turbulent migration into the AGN disk where they form merging binaries.Therefore, if it can be verified that this pathway is re-sponsible for the observed anti-correlation in (q, χ eff ), it may be possible infer properties of AGN environments using BBH mergers (Vajpeyi et al. 2022).
Moving forward, it will be beneficial to construct population models that have the capacity to distinguish between the various scenarios outlined above.Doing so may provide evidence for, or falsify, the aforementioned hypotheses on the origin of the purported (q, χ eff ) anticorrelation -therefore informing our understanding of BBH formation mechanisms more generally.
where N (m 1 |µ m , σ m ) is a Gaussian distribution with mean µ m and width σ m , λ m is the fraction of masses in this Gaussian component, and is a power-law distribution with spectral-index α.The primary mass m 1 is restricted to lie between hyper-parameters [m min , m max ] using a Heaviside step function Θ at the upper mass limit and a smoothing function at low mass: where δ m is the smoothing length.The mass ratio is then power-law distributed with a dependence on m 1 : Usually, upon constructing a population model, one can infer the population hyper-parameters by inserting it into, then minimising, the likelihood function defined in equation 32 of Thrane & Talbot (2019).This likelihood function effectively re-weights individual event parameter samples -which have been obtained using an uninformative fiducial prior -according to the newly proposed population model.These re-weighted samples are then combined to form a likelihood.However, the aforementioned fiducial samples generally do not have sufficient prior volume around χ 1,2 = 0, meaning population-level inferences would be biased away from cases of zero-spin.
To remedy this, we use the methodology from Galaudage et al. (2021) and Tong et al. (2022), in which the likelihood uses two sets of fiducial samples.One set is produced using the standard uninformative prior on χ 1,2 , being uniform on [0, 1].The other is produced by setting χ 1,2 = 0. We then create a pair of corresponding, event-level likelihoods.For the non-spinning case, given an event i with associated data d i , we have where Z i 0 is the evidence obtained during fiducial sampling for event i, assuming χ 1,2 = 0, n is the number of samples re-weighted, π 0 (θ i k |Λ) is the model defined in equation 15 for the non-spinning sub-population, evaluated at a sample θ i k , and π Ø (θ i k |χ 1,2 = 0) is the fiducial prior used in non-spinning sample generation, evaluated for the same sample.For the spinning case, we then have the familiar likelihood where Z i s is the fiducial evidence for event i, assuming a uniform spin prior, π s (θ i k |Λ, κ) is the spinning sub-population model defined in equation 14, evaluated at a sample θ i k , and π Ø (θ i k ) is the fiducial prior with uniform spin, evaluated at the same sample.
We can then combine these event-level likelihoods with a mixing fraction λ 0 , and produce a total likelihood by taking the product over all N events Here, the factor of p det (Λ) −N accounts for mass and redshift-based effects on BBH detection efficiency (see Messenger & Veitch 2013;Thrane & Talbot 2019;Mandel et al. 2019;Abbott et al. 2019Abbott et al. , 2021bAbbott et al. , 2023)).We opt to exclude selection effects based on BBH spins as there are a number of complications that arise due to sampling issues about spin magnitudes of χ = 0.These spin-based selection biases are generally thought to be minor (Abbott et al. 2023) and should not have a significant impact on our inferences.Plots consist of 2 × 10 4 draws from the population predictive distribution for mass ratio and χ eff , assuming the maximum posterior probability hyper-parameters.These plots use results from analyses excluding GW191109.Note that samples drawn from the non-spinning sub-population have been randomly off-set by χ eff ∼ ±0.005 in order to help with the visualisation.

D. IMPACT OF ALTERNATIVE COPULAS
When choosing a copula density function, one must consider the types of correlations that the given function has the flexibility for.Currently in the literature, there appear to be only three copula density functions that allow for a smooth transition from correlation to anti-correlation by varying κ.These are the Frank copula density function (Frank 1979), Gaussian copula density function (Klaassen & Wellner 1997), and Farlie-Gumbel-Morgenstern (FGM) copula density function (Farlie 1960;Gumbel 1960;Morgenstern 1956).Here, we investigate the effect of replacing the Frank copula density function, assumed in Section 2, with the two latter copulas.That is, we replace the function π c (u, v|κ) in equation 14, and re-run our analyses to get alternative inferences on the correlation κ.We do so excluding GW191109 from our dataset.
Firstly, we explore the effects of a Gaussian copula on our inferences.A Gaussian copula density function can be given as and u and v are as defined in equations 12 and 13.Here, κ is restricted to [−1, 1], where κ = 1 implies maximum correlation and κ = −1 implies maximum anti-correlation.6With a Gaussian copula, we find anti-correlation (κ < 0) with 97.7% credibility.Compared to the Frank copula model, the Gaussian copula model is disfavoured by a natural log Bayes factor of only ln BF = 0.1.However, the maximum natural log likelihood of the Frank copula model is ∆ ln L max = 1.2 higher.This implies that the Gaussian copula model is favoured over the no-correlation hypothesis by a natural log Bayes factor of ln BF = 2.1 (∆ ln L max = 5.2).The (q, χ eff ) population distribution for this model is reconstructed using the maximum posterior probability hyper-parameters in the left panel of Fig. 4. We can see that the Gaussian copula forces the distribution of χ eff to narrow substantially as q approaches one, implying only a small range of χ eff for equal-mass binaries.This feature may explain the data's slight preference for the Frank copula model.Next, we experiment with an FGM copula density function, defined as where κ is again defined on [−1, 1], with κ = 1 and κ = −1 giving maximum correlation and anti-correlation respectively.With an FGM copula, we find anti-correlation (κ < 0) with 79.1% credibility.Again comparing to the Frank copula model, the FGM model is disfavoured by a natural log Bayes factor of ln BF = 1.8 (difference in maximum natural log likelihood of ∆ ln L max = 2.2).Hence, the FGM model is less distinguishable from the uncorrelated model, favoured by a natural log Bayes factor of only ln BF = 0.3, but with a notable improvement in maximum natural log likelihood ∆ ln L max = 4.3.The (q, χ eff ) population distribution for this model is again plotted assuming the maximum posterior probability hyper-parameters in the right panel of Fig. 4. We can see that correlations manifest themselves more subtly using an FGM copula -making it unsurprising that this model is less distinguishable from an uncorrelated model.

E. FULL POSTERIOR AND POPULATION DISTRIBUTIONS
Here, we share results that are not relevant to the key question of this study (are mass ratio and χ eff correlated?), but may still be of interest for some readers.
In Fig. 5, we plot the inferred posterior distributions of the hyper-parameters Λ, which govern the shapes of the marginal population distributions.This includes the fraction of binaries with non-spinning black holes λ 0 , which we find to be consistent with Tong et al. (2022).That is, while λ 0 > 0 is preferred, the data is consistent with λ 0 = 0.
We plot the full reconstructed population distributions of mass ratio and effective spin parameters in Fig. 6, as well as the implied distributions of physical spin parameters in Fig. 7. .Implied population predictive distribution of mass ratio q, and physical spin parameters χ1,2 and cos t1,2, reconstructed by marginalising over the hyper-parameter posterior distributions.From darkest to lightest, the shaded regions in the two-dimensional plots encapsulate 50, 90, and 99% of the inferred BBH population.This plot uses results excluding GW191109.

Figure 1 .
Figure 1.Comparison between the Extended population model from Tong et al. (2022) (blue), and this work (orange).In this demonstration, the hyper-parameters for our model are chosen to resemble the Extended model distributions in blue.The correlation parameter κ is set to zero in this example.Varying κ changes the degree of correlation in the (q, χ eff ) plane, but leaves every other panel unchanged.The non-spinning sub-population is also excluded from this demonstration.

Figure 3 .
Figure 3. Plot of 3 × 10 4 draws from the population predictive distribution for mass ratio and χ eff , marginalised over the posterior distributions of the population hyperparameters and excluding GW191109.In red, is the purely data-driven π from Payne & Thrane (2023), consisting of a set of delta functions that attempt to maximise the population likelihood.Darker tones indicate a higher weight, thus higher importance, for the given feature.Our model matches qualitative features present in the π result.Note that samples drawn from the non-spinning sub-population have been randomly off-set by χ eff ∼ ±0.005 in order to help with the visualisation.

Figure 4 .
Figure 4. Comparison of (q, χ eff ) population distributions for the Gaussian copula model (left) and FGM copula model (right).Plots consist of 2 × 10 4 draws from the population predictive distribution for mass ratio and χ eff , assuming the maximum posterior probability hyper-parameters.These plots use results from analyses excluding GW191109.Note that samples drawn from the non-spinning sub-population have been randomly off-set by χ eff ∼ ±0.005 in order to help with the visualisation.

Figure 5 .Figure 6 .
Figure5.Posterior distributions for hyper-parameters governing the marginal distributions relevant to this study.Values listed above the posteriors are the medians, accompanied by 90% credible intervals.These intervals are plotted as dotted lines in the one-dimensional posteriors.From darkest to lightest, the shaded regions in the two-dimensional plots represent 50, 90, and 99% credibility regions.This plot uses results excluding GW191109.