A Trend in the Effective Spin Distribution of LIGO Binary Black Holes with Mass

The spin probability distribution function of stellar mass black holes at birth is unknown. The efficiency of angular momentum (AM) transfer from the core of a dying star to outer shell layers through magnetic fields sets the expected spin of the newly born compact object. The efficiency of the mechanism is debated in the literature, for example, the Geneva stellar evolution model (Eggenberger et al. 2007; Ekström et al. 2011) assumes moderate efficiency of AM transport through meridional currents and therefore permits BHs to be born with nonnegligible spin, while efficient transport by the Tayler–Spruit magnetic dynamo (Spruit 1999, 2001), as implemented in stellar evolution calculations (Fuller et al. 2019; Fuller & Ma 2019) predicts all isolated BHs to be born very slowly rotating.

The models in which BHs are assumed to be nonspinning at birth have a difficult time explaining the observed high spin of the BHs in high-mass X-ray binaries (HMXBs) (Batta et al. 2017; Qin et al. 2019) since these black holes are wind fed and therefore gas accretion (Fragos & McClintock 2014) or tidal locking (Zaldarriaga et al. 2017b; Schrøder et al. 2018) cannot explain their high spins. Moreover, one of the BHs in GW 151226 has spin greater than 0.2 (Abbott et al. 2016a) which might challenge the zero spin scenario. Given that the BHs in HMXBs are less massive compared to the LIGO BHs, it is possible there may be a mass trend for the spin of the black holes. In addition to the effects of AM transport discussed above, such a trend could be due to the supernova explosion mechanisms that form the BHs or secondary astrophysical mechanisms such as tidal locking (Zaldarriaga et al. 2017b; Schrøder et al. 2018) or gas accretion that can change the spin of a BH (Fragos & McClintock 2014). This latter mechanism is operative in the case of low-mass X-ray binaries.

The effective spin of a binary black hole (BBH) system is defined as

$$\begin{eqnarray}&&{\chi }_{\mathrm{eff}}\equiv \displaystyle \frac{{m}_{1}{a}_{1}\cos ({\theta }_{1})+{m}_{2}{a}_{2}\cos ({\theta }_{2})}{{m}_{1}+{m}_{2}},\end{eqnarray} \tag{ 1 }$$

where m₁ and m₂ are the masses of the primary and secondary black holes, and a₁ and a₂ are their associated dimensionless spin magnitudes defined as:

$$\begin{eqnarray}&&a=\displaystyle \frac{{{cJ}}_{\mathrm{BH}}}{{{GM}}_{\mathrm{BH}}^{2}}.\end{eqnarray} \tag{ 2 }$$

Here c is the speed of light, G is the gravitational constant, and M and J are the mass and AM of the BH. θ is the angle between the direction of each BH's spin and the orbital AM of the BBH. The effective spin parameter is the best-measured spin-related parameter from gravitational wave observations (Farr et al. 2017, and references therein), so here we focus on this one-dimensional summary of the full, six-dimensional space of BBH spins.

The spin distribution of the LIGO black holes therefore carries crucial information that illuminates the formation process of these systems (Farr et al. 2017; Stevenson et al. 2017; Vitale et al. 2017b). Broadly, two different mechanisms have been proposed for the formation of the BBHs: (i) assembled in the field through stellar evolution and a potential common envelope phase, (ii) assembled dynamically, either in globular or nuclear star clusters or hierarchical triple or higher order stellar systems (Rodriguez & Antonini 2018; Safarzadeh et al. 2020). Each of these channels predicts a different spin–mass distribution: field binaries are expected have their BH spins preferentially aligned with the orbital AM of the binary (Belczynski et al. 2002; Dominik et al. 2012; Zaldarriaga et al. 2017b; Gerosa et al. 2018; Schrøder et al. 2018; Qin et al. 2019; Bavera et al. 2020), while dynamically assembled binaries (Zwart et al. 2004; Samsing et al. 2014, 2018; Chatterjee et al. 2016; Rodriguez et al. 2016, 2018; Antonini et al. 2017) are expected to have their spin isotropically distributed with respect to the AM of the binary and therefore result in a large fraction of the systems having ${\chi }_{\mathrm{eff}}$  < 0.

The effective spin parameter for the 10 LIGO/Virgo GWTC-1 BBHs is consistent with being clustered all around zero (Abbott et al. 2019a; Roulet & Zaldarriaga 2019; Belczynski et al. 2020) which could be due the fact that LIGO black holes are mostly nonspinning or their spins lie in the orbital plane of the binary. Here we examine this population for mass-dependent effects on effective spin.

Several additional BBH merger events have been claimed in the same LIGO/Virgo data used to generate GWTC-1 (Venumadhav et al. 2019, 2020), including some with (large) positive (Zackay et al. 2019) and negative effective spin (Venumadhav et al. 2020). Piran & Piran (2020) argue that the larger catalog is more consistent with field than dynamical formation (using models where the entire population comes from a single channel). Here we focus only on the BBH systems in GWTC-1 for two reasons: (1) full posterior distributions for the parameters of the additional events have not been made available, and a Gaussian approximation to the ${\chi }_{\mathrm{eff}}$ posterior may not be adequate for population analysis in such a large catalog (Ng et al. 2018) and (2) there is no publicly released procedure to characterize the sensitivity of the pipelines used in Venumadhav et al. (2019, 2020) well enough to account for selection biases in the population.

We summarize how different formation channels of the BBHs populate different regions in the ${\chi }_{\mathrm{eff}}$ –mass plane in Section 2. In Section 3 we analyze the joint mass-effective spin distribution of the 10 LIGO/Virgo BBHs to search for possible correlations of the mean and dispersion of the effective spin with mass, where mass can be either the primary mass, the chirp mass, or the total mass of the binary. In Section 4 we summarize our results and suggest alternative joint distribution studies that could carry similar information as ${\chi }_{\mathrm{eff}}$ –mass distribution.

Different models for the evolution of field binaries predict different ${\chi }_{\mathrm{eff}}$ –mass distribution: there are models that predict all isolated black holes should be born slowly rotating (e.g., Fuller et al. 2019; Fuller & Ma 2019) and, therefore, secondary astrophysical mechanisms such as tidal interactions are invoked to explain fast rotating BHs in either HMXBs or GW 151226. Effectively such models predict a distribution of BHs in mass–spin at birth similar to the blue band in Figure 1 for the BBHs.

Zoom In Zoom Out Reset image size

In the case of moderate efficiency of AM as implemented in the MESA stellar evolution model (Eggenberger et al. 2007), low metallicity stars are expected to not lose mass through winds, as their opacity for EUV/UV photons is small (Kudritzki & Puls 2000; Vink et al. 2001). Therefore, the collapse of such stars is expected to result in both massive, and highly spinning BHs (although in such cases feedback is likely to limit their mass; Batta & Ramirez-Ruiz 2019). The predicted locus of such objects are depicted by the green circle in Figure 1. Since low metallicity environments are thought to be the underlying requirement for the formation of such BBHs, given the metallicity evolution of the universe, these systems are expected to be likely born at high redshifts (although it is possible to form such systems in pockets of the low metallicity regions in the local universe), and therefore, a long delay time (large separations at birth) are thought to make them merge at z < 0.2 such that LIGO can see them. This channel has been proposed to explain GW 170729 (Bavera et al. 2020).

Other models with inefficient AM transport can produce a distribution in ${\chi }_{\mathrm{eff}}\,$ space with a mean value that decreases and a dispersion that increases with increasing mass (Belczynski et al. 2020), as we find in the GWTC-1 catalog. However, such models generally predict higher values of ${\chi }_{\mathrm{eff}}\,$ at low masses than are observed in GWTC-1, and are therefore disfavored (Belczynski et al. 2020). Nevertheless, it may be possible to produce the trends we observe in mean and dispersion with mass through field formation models invoking a combination of efficient AM transport, tidal spin-up, and metallicity-dependent mass loss from high-mass stellar winds that differs from the one explored in Bavera et al. (2020).

Dynamical assembly can also make massive spinning BBH mergers; however, the expected distribution in the case of dynamical assembly is symmetric and therefore future BBH detections can tell us whether there is a locus of BBHs at high mass and high spin, or whether the massive BBHs are symmetrically distributed in ${\chi }_{\mathrm{eff}}$ . This can shed light on the underlying formation mechanism of such systems and whether their spins have been altered by subsequent stellar encounters (e.g., Lopez et al. 2019).

Under dynamical assembly, the predicted ${\chi }_{\mathrm{eff}}$ –mass distribution is more well defined in its structure: first generation black holes could be born with zero spin; however, the final merger product of such black holes will have high spins (a ≈ 0.7) (Fishbach et al. 2017; Gerosa & Berti 2017). The merger of these second generation BHs with either other second generation BHs, or first generation BHs, forms a BBH with at least one of the BHs highly spinning. The spin orientation of the BHs in the dynamical assembly would be random with respect to the binary's orbital plane and therefore the final merger product is expected to show a symmetric ${\chi }_{\mathrm{eff}}$ distribution around zero that widens at larger masses. The widening at larger masses is due to the merger products of higher generation black holes in a dense cluster-like environment (Rodriguez et al. 2019; Doctor et al. 2020).

Figure 2 shows the expected distribution in ${\chi }_{\mathrm{eff}}$ –mass from two separate categories: field binaries shown as the blue banana region, and the dynamically assembled binaries shown with the red pear-like region. These are rough sketches of the expected distributions and not necessarily to scale. The lower orange band indicates the debated lower mass gap (between 2 and 5 ${M}_{\odot }$). The green band and the light green region above it indicates the presence of a drop in mass function of the BBHs given the 10 LIGO/Virgo BBHs (Fishbach & Holz 2017; Talbot & Thrane 2018; Roulet & Zaldarriaga 2019; Safarzadeh & Farr 2019) due to pulsational pair instability supernovae (Woosley 2017) where BHs with mass between ≈50 and 150 ${M}_{\odot }$ are expected to not form. The field binaries can provide the negative trend with mass, and the dynamically assembled binaries provide the increase in dispersion with mass, the combination of which can potentially explain our findings.

Zoom In Zoom Out Reset image size

In this section we analyze the joint ${\chi }_{\mathrm{eff}}$ –mass distribution of the 10 LIGO/Virgo BBHs to search for possible correlations of the mean and dispersion of the effective spin with mass, where mass can be either the primary mass, the chirp mass, or the total mass of the binary.

We assume that the population distribution for χ_eff, conditioned on mass, follows

$$\begin{eqnarray}&&p\left({\chi }_{\mathrm{eff}}| m\right)={ \mathcal N }\left(\mu \left(m\right),\sigma \left(m\right)\right)T[-1,1],\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal N }\left(\mu ,\sigma \right)$ is a normal distribution with mean μ and dispersion σ. The notation T[−1, 1] means to truncate the normal distribution to the range [−1, 1]. The parameters μ and σ are set by the mass of the BBH via

$$\begin{eqnarray}&&\mu (m)={\mu }_{0}+\alpha \times \left(\displaystyle \frac{m}{30\,{M}_{\odot }}-1\right),\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&\sigma (m)={\sigma }_{0}\exp \left(\beta \times \left(\displaystyle \frac{m}{30\,{M}_{\odot }}-1\right)\right).\end{eqnarray} \tag{ 5 }$$

Here m can be any measure of the mass of the BBHs: the primary mass of the BBH (m₁), the total mass of the system (M_tot = m₁ + m₂), or the chirp mass of the system,

$$\begin{eqnarray}&&{M}_{\mathrm{chirp}}={\eta }^{3/5}{M}_{\mathrm{tot}},\end{eqnarray} \tag{ 6 }$$

where η is the symmetric mass ratio given by

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{m}_{1}{m}_{2}}{{({m}_{1}+{m}_{2})}^{2}}.\end{eqnarray} \tag{ 7 }$$

We assume a fixed mass and redshift distribution consistent with the population analysis in Abbott et al. (2019a) throughout this work, with

$$\begin{eqnarray}&&p\left({m}_{1}\right)\propto \displaystyle \frac{1}{{m}_{1}},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&p\left({m}_{2}| {m}_{1}\right)\propto \mathrm{const}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&p\left(z\right)\propto {\left(1+z\right)}^{1.7}\displaystyle \frac{{dV}}{{dz}},\end{eqnarray} \tag{ 10 }$$

where V(z) is the comoving volume (Hogg 1999). (The redshift distribution gives an observed merger rate that is consistent with the comoving merger rate tracking the low-redshift evolution of the star formation rate (Madau & Dickinson 2014; Fishbach et al. 2018).) We have verified that our results do not depend on the choice of mass and redshift distribution within the ranges permitted by the population analysis of Abbott et al. (2019a).

The marginal likelihood for the catalog data d_GW given population parameters ${\rm{\Theta }}\equiv \left({\mu }_{0},{\sigma }_{0},\alpha ,\beta \right)$ is (Mandel et al. 2018)

$$\begin{eqnarray}\begin{array}{rcl}p\left({d}_{\mathrm{GW}}| {\rm{\Theta }}\right) & = & \displaystyle \prod _{i=1}^{N}\displaystyle \frac{1}{\alpha \left({\rm{\Theta }}\right)}\displaystyle \int d{\chi }_{\mathrm{eff}}^{(i)}\,{{dm}}_{1}^{(i)}\,{{dm}}_{2}^{(i)}\,{dz}\\ & & \times \,p\left({d}_{\mathrm{GW}}^{(i)}| {\chi }_{\mathrm{eff}}^{(i)},{m}_{1},{m}_{2},z\right)\\ & & \times \,p\left({\chi }_{\mathrm{eff}}^{(i)}| {\rm{\Theta }},{m}_{1},{m}_{2}\right)p\left({m}_{1},{m}_{2},z\right),\end{array}\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}\alpha \left({\rm{\Theta }}\right) & = & \displaystyle \int d{\chi }_{\mathrm{eff}}\,{{dm}}_{1}\,{{dm}}_{2}\,{dz}\,{P}_{\det }\left({\chi }_{\mathrm{eff}},{m}_{1},{m}_{2},z\right)\\ & & \times \,p\left({\chi }_{\mathrm{eff}}| {\rm{\Theta }},{m}_{1},{m}_{2}\right)p\left({m}_{1},{m}_{2},z\right)\end{array}\end{eqnarray} \tag{ 12 }$$

is the population-averaged detection probability.

Here we model the detection process semianalytically, using a method similar to the one described in Abbott et al. (2016b), but with a three-detector network (two Advanced LIGO, one Advanced Virgo, all assumed to operate with "early high-sensitivity" noise (Abbott et al. 2018)) and a correspondingly enhanced (noisy) SNR threshold of $\rho \gt 8\sqrt{2}\simeq 11.3$. For all these calculations we use the IMRPHENOMPV2 waveform family (Hannam et al. 2014). The detection probabilities produced for ${\chi }_{\mathrm{eff}}$ , m₁, and m₂ from our analytic model are shown in Figure 3.

Zoom In Zoom Out Reset image size

The integral in the product can be approximated as a weighted sum over samples drawn from the likelihood function $p\left({d}_{\mathrm{GW}}^{(i)}| {\chi }_{\mathrm{eff}}^{(i)},{m}^{(i)}\right)$:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \int d{\chi }_{\mathrm{eff}}\,{{dm}}_{1}\,{{dm}}_{2}\,{dz}\,p\left({d}_{\mathrm{GW}}^{(i)}| {\chi }_{\mathrm{eff}},{m}_{1},{m}_{2},z\right)\\ \quad \times \,p\left({\chi }_{\mathrm{eff}}| {\rm{\Theta }},{m}_{1},{m}_{2}\right)p\left({m}_{1},{m}_{2},z\right)\\ \quad \propto \,\displaystyle \frac{1}{{N}_{\mathrm{samp}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{samp}}}p\left({\chi }_{\mathrm{eff}}^{(j)}| {\rm{\Theta }},{m}_{1}^{(j)},{m}_{2}^{(j)}\right)p\left({m}_{1},{m}_{2},z\right).\end{array}\end{eqnarray} \tag{ 13 }$$

Similarly, α can be approximated as a weighted sum over samples drawn from a canonical distribution, ${p}_{\mathrm{draw}}\left({\chi }_{\mathrm{eff}},{m}_{1},{m}_{2},z\right)$, and "detected" by our synthetic pipeline (Farr 2019):

$$\begin{eqnarray}&&\alpha \left({\rm{\Theta }}\right)\propto \sum _{k=1}^{{N}_{\mathrm{detected}}}\displaystyle \frac{p\left({\chi }_{\mathrm{eff}}\,| {\rm{\Theta }},{m}_{1},{m}_{2}\right)p\left({m}_{1},{m}_{2},z\right)}{{p}_{\mathrm{draw}}\left({\chi }_{\mathrm{eff}},{m}_{1},{m}_{2},z\right)}.\end{eqnarray} \tag{ 14 }$$

We reweight the posterior samples of ${\chi }_{\mathrm{eff}}$ provided by the LVC for each of the BBH systems in GWTC-1 (Abbott et al. 2019b) by the inverse of the LALInference prior (Veitch et al. 2015) to draw samples from the likelihood function. The default prior from LALInference assumes a uniform mass distribution for m₁ and m₂ in the detector frame, uniform distribution for the magnitude of the spin parameter for each black hole, and uniform prior on the $\cos (\theta )$ from −1 to 1, and a prior on the luminosity distance proportional to d_L². Figure 4 shows the difference between the posterior samples from the GWTC-1 catalog and the likelihood function for ${\chi }_{\mathrm{eff}}\,$ marginalized over our assumed mass and redshift distribution (Equations (8)–(10)).

Zoom In Zoom Out Reset image size

We apply a flat prior on the population parameters Θ, and sample from their posterior distribution given the catalog data (i.e., we draw samples of Θ proportional to the function in Equation (11)) using the emcee stochastic sampler (Foreman-Mackey et al. 2013) with various choices of mass parameter controlling the ${\chi }_{\mathrm{eff}}\,$ distribution. Our results are summarized in Figures 5–7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The left panel of Figure 5 shows a contour plot of the posterior ${\chi }_{\mathrm{eff}}\,$ distribution in the ${\chi }_{\mathrm{eff}}$ –m_chirp plane after analyzing the GWTC-1 BBHs and marginalizing over Θ:

$$\begin{eqnarray}&&p\left({\chi }_{\mathrm{eff}}| {d}_{\mathrm{GW}}\right)=\int d{\rm{\Theta }}\,p\left({\chi }_{\mathrm{eff}}| {\rm{\Theta }}\right)p\left({\rm{\Theta }}| {d}_{\mathrm{GW}}\right).\end{eqnarray} \tag{ 15 }$$

The crosses show the 1σ error in chirp mass and ${\chi }_{\mathrm{eff}}$ of the 10 LIGO BBHs. The blue cross represents GW 170729, which has the highest observed ${\chi }_{\mathrm{eff}}$ among the 10 BBH systems. The red contours show the posterior PDF, $p\left({\chi }_{\mathrm{eff}}| {m}_{\mathrm{chirp}},{d}_{\mathrm{GW}}\right)$ without considering GW 170729, while the blue contours show the result when all of the 10 LIGO events are analyzed. We have singled out GW 170729 since it has the highest significant false alarm rate (FAR) of 0.18 yr⁻¹ in the GstLAL pipeline (Abbott et al. 2019a), which nearly matches the threshold of 0.1 yr⁻¹ , while the other BBH systems have significantly lower FARs (<10⁻³).

The right panel of Figure 5 shows the posterior PDF of the four parameters in our model used to describe the ${\chi }_{\mathrm{eff}}$ –mass relation. A negative trend with mass (α < 0) and a positive trend of ${\chi }_{\mathrm{eff}}$ dispersion with mass (β > 0) are favored at 95% and 60%, respectively, when ignoring GW 170729, and 74% and 78% when analyzing all 10 LIGO events.

Figure 6 shows the same result but when the mass variable is taken to be the primary mass of the system. When the mass scale is the primary mass of the BBHs, a negative trend with mass (α < 0) and a positive trend of ${\chi }_{\mathrm{eff}}$ dispersion with mass (β > 0) is favored at 80% and 90%, respectively, when ignoring GW 170729, and 70% and 95% when analyzing all 10 LIGO events.

Figure 7 shows the same result but when the mass is taken to be the total mass of the system. When the mass scale is the total mass of the BBHs, a negative trend with mass (α < 0) and a positive trend of ${\chi }_{\mathrm{eff}}$ dispersion with mass (β > 0) is favored at 30% and 60%, respectively, when ignoring GW 170729, and 70% and 70% when analyzing all 10 LIGO events.

Although the choice of mass scale somewhat affects the confidence by which a trend with mass for either the mean ${\chi }_{\mathrm{eff}}$ or its dispersion could be detected, the data suggests that if there is a trend, it is a negative trend with mass for the mean ${\chi }_{\mathrm{eff}}$ and a positive trend for its dispersion. Dynamical assembly alone can account for the larger dispersion in ${\chi }_{\mathrm{eff}}$ with mass; however, the negative trend of the mean ${\chi }_{\mathrm{eff}}$ with mass cannot be accounted for based on dynamical assembly. Field formation of the BBHs on the other hand can potentially explain the negative trend with mass (through a combination of tidally spun-up systems at low masses and massive BBHs with about zero effective spin at birth); however, the increase of dispersion with mass would be hard to accommodate based on field evolution alone. Thus we suggest that the observed trends in ${\chi }_{\mathrm{eff}}\,$ with mass could be indicating the operation of both formation channels in the GWTC-1 observations.

A given formation channel for a BBH system would predict a certain distribution in ${\chi }_{\mathrm{eff}}$ –mass plane for the final merger event that LIGO/Virgo would observe. Field binaries tend to predict a banana shaped region that encompasses massive systems with negligible ${\chi }_{\mathrm{eff}}$ magnitude or low-mass systems with positive ${\chi }_{\mathrm{eff}}$ magnitudes if the AM transport is weak. The distribution could be a combination of three main mechanisms: (i) formation with negligible spin at all masses due to efficient AM transport, (ii) tidal spin of the binaries in which the second born compact object forms from a tidally spun-up progenitor in a close orbit with another BH, (iii) formation with moderate efficiency of AM transport while including feedback from disk formation after the core of the progenitor star has collapsed and formed a BH.

On the other hand, dynamical assembly of BBHs in dense stellar clusters leads to a symmetric distribution of BBHs in ${\chi }_{\mathrm{eff}}$ at all masses with larger dispersion at higher masses. The increase of dispersion is due to a random walk in ${\chi }_{\mathrm{eff}}$ –mass that higher generation BHs follow.

In this study we find a tentative negative correlation between ${\chi }_{\mathrm{eff}}$ and chirp mass for the 10 LIGO/Virgo BBHs with ∼75% confidence. The negative correlation could in principle be explained by field formation alone. However, standard field formation consistent with the observed small effective spins at low mass would predict that the dispersion should decrease with mass, the opposite of what the data suggests. We find that the dispersion in ${\chi }_{\mathrm{eff}}$ grows with mass with 80% confidence. These trends are consistent with a combined channel of dynamically assembled BBHs that provide the positive trend of dispersion with mass, and a field formation channel that provides the negative mean trend with mass could explain our findings.

Given the public alerts released by Advanced LIGO and Virgo⁵ in the first half of the third observing run, the O3a catalog should contain ∼30 BBH mergers; therefore, the statistical uncertainties on our parameters Θ should decrease by about a factor of $\sqrt{{N}_{{\rm{O}}3{\rm{a}}}/{N}_{{\rm{O}}2}}\sim 2$ on incorporating the O3a catalog; additional detections in O3b (Abbott et al. 2018, 2019b) should produce a further factor of $\sim \sqrt{2}$ reduction in uncertainty. Thus we can expect the O3 catalog of BBH systems to confidently confirm or overturn the trends observed here.

M.T.S. is grateful to the Center for Computational Astrophysics for hospitality during the course of this work. E.R.R. and M.T.S. thank the Heising-Simons Foundation, the Danish National Research Foundation (DNRF132) and NSF (AST-1911206 and AST-1852393) for support. The authors thank the LIGO Scientific Collaboration for access to public data and gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, and the Max Planck-Society (MPS) for support of the construction of Advanced LIGO. Additional support for Advanced LIGO was provided by the Australian Research Council. We are thankful to Daniel Wysocki, Richard O'shaughnessy, and Jim Fuller for useful comments on the earlier version of this work.

Software: Numpy (Walt et al. 2011), Scipy (Virtanen et al. 2020), IPython (Pérez & Granger 2007), Matplotlib (Hunter 2007), astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), PyStan (Carpenter et al. 2017; Stan Development Team 2018), Seaborn (Waskom et al. 2018), Arviz (Kumar et al. 2019), emcee (Foreman-Mackey et al. 2013), corner (Foreman-Mackey 2016).