Hierarchical analysis of gravitational-wave measurements of binary black hole spin-orbit misalignments

TL;DR

This paper develops a hierarchical analysis method to interpret gravitational-wave data, enabling the identification of sub-populations of binary black holes with different spin-orbit misalignments, which sheds light on their formation channels.

Contribution

It introduces a hierarchical analysis framework applied to simulated gravitational-wave data to distinguish and quantify sub-populations of binary black holes based on spin orientations.

Findings

Tens of observations can differentiate subpopulations based on spin orientations.

100 observations allow estimating the fraction of isotropic spin directions with ~40% uncertainty.

Only 5 observations suffice to distinguish extreme models of spin alignment.

Abstract

Binary black holes may form both through isolated binary evolution and through dynamical interactions in dense stellar environments. The formation channel leaves an imprint on the alignment between the black hole spins and the orbital angular momentum. Gravitational waves from these systems directly encode information about the spin--orbit misalignment angles, allowing them to be (weakly) constrained. Identifying sub-populations of spinning binary black holes will inform us about compact binary formation and evolution. We simulate a mixed population of binary black holes with spin--orbit misalignments modelled under a range of assumptions. We then develop a hierarchical analysis and apply it to mock gravitational-wave observations of these populations. Assuming a population with dimensionless spin magnitudes of $\chi = 0.7$, we show that tens of observations will make it possible to distinguish the presence of subpopulations of coalescing binary black holes based on their spin orientations. With $100$ observations it will be possible to infer the relative fraction of coalescing binary black holes with isotropic spin directions (corresponding to dynamical formation in our models) with a fractional uncertainty of $\sim 40\%$. Meanwhile, only $\sim 5$ observations are sufficient to distinguish between extreme models---all binary black holes either having exactly aligned spins or isotropic spin directions.

Equations

$$m_i=30{\left(\frac{m_i^{\mathrm{ZAMS}}}{50M_{\odot }}\right)}^{\alpha }{M}_{\odot },$$

$$m_i=0.6m_i^{\mathrm{ZAMS}}.$$

$$\cos {\theta }_i={\hat{\mathbf{S}}}_i\cdot \hat{\mathbf{L}},$$

$${\mathbf{S}}_i={\chi }_i{m}_i^2{\hat{\mathbf{S}}}_i,$$

$$d(t)=n(t)+h(\mathbf{\Theta },t).$$

$$p(\mathbf{\Theta }|d)=\frac{p(d|\mathbf{\Theta })p(\mathbf{\Theta })}{p(d)},$$

$$p(d|\mathbf{\Theta })\propto \exp \left[-\frac{1}{2}\left(d-h(\mathbf{\Theta })|d-h(\mathbf{\Theta })\right)\right],$$

$$\left(g|h\right)=4\Re {\int }_{f_{\mathrm{low}}}^{\infty }\frac{\tilde{g}(f){\tilde{h}}^{\ast }(f)}{S_n(f)}\mathrm{d}f,$$

$$p(\mathbf{\Theta }|\mathbf{\lambda })=\sum _{\ell =1}^4{\lambda }_{\ell }p(\mathbf{\Theta }|{\Lambda }_{\ell }).$$

$$\sum _{\ell =1}^4{\lambda }_{\ell }=1,$$

$$p(\mathbf{\Theta }|d_{\alpha })=\frac{p(d_{\alpha }|\mathbf{\Theta })p(\mathbf{\Theta })}{p(d_{\alpha })},$$

$$p\left({\left\{d_{\alpha }\right\}}_{\alpha =1}^N|\mathbf{\lambda }\right)$$

$$=\prod _{\alpha =1}^N\int \mathrm{d}{\mathbf{\Theta }}_{\alpha }p(d_{\alpha }|{\mathbf{\Theta }}_{\alpha })p({\mathbf{\Theta }}_{\alpha }|\mathbf{\lambda })$$

$$=\prod _{\alpha =1}^{N}p(d_{\alpha })\int \mathrm{d}{\mathbf{\Theta }}_{\alpha }\frac{p({\mathbf{\Theta }}_{\alpha }|d_{\alpha })}{p({\mathbf{\Theta }}_{\alpha })}p({\mathbf{\Theta }}_{\alpha }|\mathbf{\lambda }),$$

$$\int \mathrm{d}{\mathbf{\Theta }}_{\alpha }p({\mathbf{\Theta }}_{\alpha }|d_{\alpha })f({\mathbf{\Theta }}_{\alpha })=\frac{1}{{\nu }_{\alpha }}\sum _{k=1}^{{\nu }_{\alpha }}f({\mathbf{\Theta }}_{\alpha }^k),$$

$$p\left({\left\{d_{\alpha }\right\}}_{\alpha =1}^N|\mathbf{\lambda }\right)=\prod _{\alpha =1}^N\frac{p(d_{\alpha })}{{\nu }_{\alpha }}\sum _{k=1}^{{\nu }_{\alpha }}\frac{p({\mathbf{\Theta }}_{\alpha }^k|\mathbf{\lambda })}{p({\mathbf{\Theta }}_{\alpha }^k)}.$$

$$p\left(\mathbf{\lambda }|{\left\{d_{\alpha }\right\}}_{\alpha =1}^N\right)\propto p\left({\left\{d_{\alpha }\right\}}_{\alpha =1}^N|\mathbf{\lambda }\right)p(\mathbf{\lambda }),$$

$${\chi }_{\mathrm{eff}}=\frac{{\chi }_1\cos {\theta }_1+q{\chi }_2\cos {\theta }_2}{(1+q)}.$$

$${\chi }_{\mathrm{eff}}^{\mathrm{data}}$$

$$\cos {{\theta }_1}^{\mathrm{data}}$$

$${\chi }_{\mathrm{eff}}^{\mathrm{sample}}$$

$$\cos {{\theta }_1}^{\mathrm{sample}}$$

$${\sigma }_{\cos {\theta }_1}$$

$${\sigma }_{{\chi }_{\mathrm{eff}}}$$

$$\frac{Z_{\mathrm{aligned}}}{Z_{\mathrm{isotropic}}}=\frac{p\left({\left\{d_{\alpha }\right\}}_{\alpha =1}^N|{\lambda }_1=1\right)}{p\left({\left\{d_{\alpha }\right\}}_{\alpha =1}^N|{\lambda }_2=1\right)}.$$

Full text

††pubyear: 2017††pagerange: Hierarchical analysis of gravitational-wave measurements of binary black hole spin–orbit misalignments–Hierarchical analysis of gravitational-wave measurements of binary black hole spin–orbit misalignments

1 Introduction

Compact binaries containing two stellar-mass black holes can form as the end point of isolated binary evolution, or via dynamical interactions in dense stellar environments (see, e.g., Mandel & O’Shaughnessy, 2010; Abbott et al., 2016g, for a review). These binary black holes are a promising source of gravitational waves for ground-based detectors such as Advanced LIGO (aLIGO; Aasi et al., 2015), Advanced Virgo (AdV; Acernese et al., 2015) and KAGRA (Aso et al., 2013). Searches of data from the first observing run (Aasi et al., 2016) of Advanced LIGO (aLIGO) yielded three likely BBH coalescences: GW150914 (Abbott et al., 2016d), GW151226 (Abbott et al., 2016f) and LVT151012 (Abbott et al., 2016c, b). \colorblackA further BBH coalescence, GW170104, has been reported from the on-going second observing run (Abbott et al., 2017). GW observations give a unique insight into the properties of BBHs. We will examine one of the ways in which black hole spin measurements can be used to constrain formation mechanisms.

GW observations inform our understanding of BBH evolution in two ways: from the merger rate, and from the properties of the individual systems. The merger rate of BBHs is inferred from the number of detections; it is uncertain as a consequence of the small number of BBH observations so far. \colorblackCurrently, merger rates are estimated to be \colorblack$12$–$213~{}\mathrm{Gpc^{-3}\,yr^{-1}}$ (Abbott et al., 2016h, 2017). These rates are broadly consistent with predictions from both population synthesis models of isolated binary evolution (e.g., Voss & Tauris, 2003; Lipunov et al., 2009; Dominik et al., 2012; Postnov & Yungelson, 2014) and dynamical formation models (e.g., Sigurdsson & Hernquist, 1993; Sadowski et al., 2008; Rodriguez et al., 2015). Possible progenitors systems of BBHs, including Cyg X-3 (Belczynski et al., 2013), IC 10 X-1 (Bulik et al., 2011) and NGC 300 X-1 (Crowther et al., 2010) provide some additional limits on BBH merger rates, but extrapolation is hindered by current observational uncertainties.

The parameters of individual systems can be estimated by comparing the measured GW signal with template waveforms (Abbott et al., 2016e). The masses and spins of the BHs can be measured through their influence on the inspiral, merger and ringdown of the system (Cutler & Flanagan, 1994; Poisson & Will, 1995; Finn & Chernoff, 1993). The distribution of parameters observed by aLIGO will encode information about the population of BBHs, and may also help to shed light on their formation channels (O’Shaughnessy et al., 2008; Kelley et al., 2010; Gerosa et al., 2013, 2014; Mandel et al., 2015; Stevenson et al., 2015; Mandel et al., 2017; Vitale et al., 2017)

Stellar-mass BHs are expected to be born spinning, with observations suggesting their dimensionless spin parameters $\chi$ take the full range of allowed values between [math] and $1$ (McClintock et al., 2011; Fragos & McClintock, 2015; Miller & Miller, 2015). Stars formed in binaries are expected to have their rotational axis aligned with the orbital angular momentum (e.g., Boss, 1988; Albrecht et al., 2007), although there is observational evidence this is not always the case (e.g., Albrecht et al., 2009; Albrecht et al., 2014). Even if binaries are born with misaligned spins, there are many processes in binary evolution which can act to align the spin of stars, such as realignment during a stable mass accretion phase (Bardeen & Petterson, 1975; Kalogera, 2000; King et al., 2005), accretion onto a BH passing through a common-envelope (CE) event (Ivanova et al., 2013), and realignment through tidal interactions in close binaries (e.g., Zahn, 1977; Hut, 1981).

On the other hand, asymmetric mass loss during supernova explosions can tilt the orbital plane in binaries (Brandt & Podsiadlowski, 1995; Kalogera, 2000), leading to BH spins being misaligned with respect to the orbital angular momentum vector. Population synthesis studies of BH X-ray binaries predict that these misalignments are generally small (Kalogera, 2000), with Fragos et al. (2010) finding that the primary BH is typically misaligned by $\lesssim 10\degree$. However, electromagnetic observations of high mass X-ray binaries containing BHs have hinted that the BHs may be more significantly misaligned (Miller & Miller, 2015). One such system is the microquasar V4641 Sgr (Orosz et al., 2001; Martin et al., 2008) where the primary BH is has been interpreted to be misaligned by $>55\degree$.

Alternatively, BBHs can form dynamically in dense stellar environments such as globular clusters. In these environments, it is expected that the distribution of BBH spin–orbit misalignment angles is isotropic (e.g., Rodriguez et al., 2016). The distribution of BBH spin–orbit misalignments therefore contains information about their formation mechanisms.

Constraints on spin alignment from GW observations so far are weak (Abbott et al., 2016e, a, b, 2017). Some configurations, such as anti-aligned spins for GW151226 (Abbott et al., 2016f), are disfavoured; however, there is considerable uncertainty in the spin magnitude and orientation. Determining the spins precisely is difficult because their effects on the waveform can be intrinsically small (especially if the the source is viewed face on), and because of degeneracies between the spin and mass parameters (Poisson & Will, 1995; Baird et al., 2013; Farr et al., 2016). Although the spins of individual systems are difficult to measure, here we show it is possible to use inferences from multiple systems to build a statistical model for the population (cf. Vitale et al., 2017).

This paper describes how to combine posterior probability density functions on spin–orbit misalignment angles from multiple GW events to explore the underlying population.

\color

blackWe develop a hierarchical analysis in order to combine multiple GW observations of BBH spin–orbit misalignments to give constraints on the fractions of BBHs forming through different channels. We consider different populations of potential spin–orbit misalignments, each representing different assumptions about binary formation, and use the GW observations to infer the fraction of binaries from each population. In the field of exoplanets, similar hierarchical analyses have been used to make inference on the frequency of Earth-like exoplanets from measurements of the period and radius of individual exoplanet candidates (e.g., Foreman-Mackey et al., 2014; Farr et al., 2014a). Other examples of the use of hierarchical analyses in astrophysics include modelling a population of trans-Neptunian objects (Loredo, 2004), measurements of spin–orbit misalignments in exoplanets (Naoz et al., 2012), measurements of the eccentricity distribution of exoplanets (Hogg et al., 2010) and the measurement of the mass distribution of galaxy clusters (Lieu et al., 2017).

In Section 2 we introduce our simplified population synthesis models for BBHs, paying special attention to the BH spins. We briefly describe in Section 3 the parameter estimation (PE) pipeline that will be employed to infer the properties (such as misalignment angles) of real GW events, and discuss previous spin-misalignment studies in the literature. We introduce a framework for combining posterior probability density functions on spin–orbit misalignment angles from multiple GW events to explore the underlying population in Section 4. We demonstrate the method using a set of mock GW events in Section 5, and show that tens of observations will be sufficient to distinguish subpopulations of coalescing binary black holes, \colorblackassuming spin magnitudes of $\sim 0.7$. Lower typical BH spin magnitudes would reduce the accuracy with which the spin-orbit misalignment angle can be measured, therefore requiring more observations to extract information about subpopulations. We also show that more extreme models, such as the hypothesis that all BBHs have their spins exactly aligned with the orbital angular momentum, can be ruled out at a $5~{}\sigma$ confidence level with only $\mathcal{O}(5)$ observations of rapidly spinning BBHs. Finally, we conclude and suggest areas which require further study in Section 6.

2 BH spin misalignment models

Owing to the many uncertainties pertaining to stellar spins and their evolution in a binary, and the fact that keeping track of stellar spin vectors can be computationally intensive, many population synthesis models choose not to include spin evolution. However, the distribution of spins of the final merging BHs is one of the observables that can be measured with the advanced GW detectors. In this section, we therefore implement a simplified population synthesis model to evolve an ensemble of binaries that will be detectable with aLIGO and Advanced Virgo (AdV) and predict their distributions of spin–orbit misalignments. We describe the assumed mass distribution, spin distribution, and spin evolution below.

2.1 Mass and spin magnitude distribution

We assume the same simplified mass distribution for all of our models, so that any differences in the final spin distributions are purely due to our assumptions about the spin–orbit misalignments described in the next section. There are many uncertainties in the evolution of isolated massive binaries, including (but not limited to) uncertainties in the initial distributions of the orbital elements (de Mink & Belczynski, 2015), the strength of stellar winds in massive stars (Belczynski et al., 2010), the effect of rotation of massive stars on stellar evolution (de Mink et al., 2013; Ramírez-Agudelo et al., 2015), the natal kicks (if any) given to BHs (Repetto et al., 2012; Mandel, 2016; Mirabel, 2016; Repetto et al., 2017) and the efficiency of the CE (Ivanova et al., 2013; Kruckow et al., 2016). Population synthesis methods are large Monte-Carlo simulations using semi-analytic prescriptions in order to explore the effect these uncertainties have on the predicted distributions of compact binaries. Instead, we adopt a number of simplifications that allow us to produce an astrophysically plausible distribution which should not, however, be considered representative of the actual mass distribution of BBHs.

We simulate massive binaries with semimajor axis $a$ drawn from a distribution uniform in $\ln a$ (Abt, 1983). The components of the binary are a massive primary BH with mass $m_{1}$ and a secondary star at the end of its main sequence lifetime.

The primary BH was formed from a massive star with zero-age main sequence (ZAMS) mass $m_{1}^{\mathrm{ZAMS}}$ drawn from the initial-mass function (IMF) with a power law index of $-2.35$ (Salpeter, 1955; Kroupa, 2001). The mass ratio of the binary at ZAMS $q^{\mathrm{ZAMS}}$ is drawn from a flat distribution $[0,1]$. The mass of the secondary star is given by $m_{2}^{\mathrm{ZAMS}}=q^{\mathrm{ZAMS}}m_{1}^{\mathrm{ZAMS}}$.

We calculate the final remnant mass $m_{i}$ as a function of the ZAMS mass $m_{i}^{\mathrm{ZAMS}}$ for each star using a fit to Figure 12 in Woosley et al. (2002). For stars with $30<m_{i}^{\mathrm{ZAMS}}/M_{\odot}<50$, in which range BHs are formed after some delay by fall-back of ejecta, we use

[TABLE]

with $\alpha=3.9$. For more massive stars with $m_{i}^{\mathrm{ZAMS}}>50M_{\odot}$, which are massive enough to directly collapse during the iron-core collapse to form BHs, we use

[TABLE]

We only consider BBHs with component masses above $10M_{\odot}$ below, consistent with aLIGO detections to date \colorblack(Abbott et al., 2016b, 2017), and hence omit stars with ZAMS masses below $30M_{\odot}$ from our population.

We assume that the binary has negligible eccentricity $e=0$, appropriate for post-CE systems. In all models we have assumed that both main-sequence stars in the binary are born with their rotation axis aligned with the orbital angular momentum axis. In general, the first supernova will misalign the spins due to any natal kick imparted on the remnant. There are expected to be mass-transfer phases between the first and second supernovae which may realign both the spins of the primary BH and the secondary BH progenitor; we vary the assumed degree of realignment in our models.

We assume that BHs receive natal kicks comparable to those received by neutron stars (Hobbs et al., 2005), namely drawn from a Maxwellian with a root-mean-square velocity of $\sim~{}250~{}\mathrm{km\,s}^{-1}$. This assumption will lead to the maximum amount of spin misalignment, and may be consistent with neutrino-driven kicks; if the natal kicks are due to asymmetric ejection of baryonic matter, then any fall-back (Fryer et al., 2012) onto BHs during formation will reduce the kick magnitude and thus the spin misalignment.

BH spins magnitudes can take any value $0\leq\chi_{i}<1$, but we set $\chi_{i}=0.7$ for all our BHs. High spin magnitudes are consistent with measurements from X-ray observations (cf. Miller & Miller, 2015), and lie toward the upper end of the range allowed by current GW observations \colorblack(Abbott et al., 2016b, 2017). Such spins are large enough to ensure spin effects on the gravitational waveform are significant, providing an opportunity for us to demonstrate our hierarchical approach, but small enough that we do not have to worry about the validity of the model gravitational waveform. Uncertainties in the relationship between pre-supernova stellar spins and BH spins mean it is not currently possible to produce a realistic distribution of spins from first principles, although a direct translation is often assumed, e.g., by Kushnir et al. (2016). If the distribution of BH spin magnitudes in nature favours smaller values, then more observations will be required to draw the conclusions we find here. The methodology we use here can be extended to models including BH spin magnitudes, which could potentially give us further information regarding formation mechanisms.

After a supernova, we establish whether the binary remains bound and, for those that do, find the new orbital elements (Blaauw, 1961; Brandt & Podsiadlowski, 1995; Kalogera, 1996). Of the remaining bound systems, we are only interested in those binaries which merge due to the emission of gravitational radiation within a Hubble time, as these are the binaries that are potentially observable with GW detectors.

2.2 Models for spin–orbit misalignment distributions

We model the overall population of BBHs as a mixture of $4$ subpopulations, each of which makes differing assumptions leading to distinct spin–orbit misalignment distributions.

We define the spin–orbit misalignment angle as the angle between the spin vector $\hat{\boldsymbol{S}}_{i}$ of binary component $i\in 1,2$ and the (Newtonian) orbital angular momentum vector $\hat{\boldsymbol{L}}$,

[TABLE]

where

[TABLE]

and $m_{i}$ is the component mass ($m_{1}\geq m_{2}$).111Throughout this paper we use geometric units $G=c=1$ unless otherwise stated. We will consider how a set of spin-misalignment measurements could be used to infer BBH formation mechanisms.

Subpopulation 1: Exactly aligned

We assume that irrespective of all prior processes, both BHs have their spins aligned with the orbital angular momentum vector after the second supernova, such that $\cos{\theta_{1}}=\cos{\theta_{2}}=1$. This may be the case if BHs receive no kicks. GW searches often assume BBHs have aligned spins as this simplification makes the search less computationally demanding (Abbott et al., 2016c).

Subpopulation 2: Isotropic/dynamical formation

We assume that BBHs are formed dynamically, such that the distribution of spin angles is isotropic. Initially isotropic distributions of spins remain isotropic (Schnittman, 2004). We still generate the binary mass distribution with our standard approach, so that the only difference in BBHs between this model and the others is the spin distribution.

Subpopulation 3: Alignment before second SN

Motivated by Kalogera (2000), we assume that the spins of both components are aligned with the orbital angular momentum after a CE event and prior to the second supernova. The tilt of the orbital plane caused by the second supernova is then taken to be the spin misalignment angle of both components, i.e. $\cos{\theta_{1}}=\cos{\theta_{2}}$. As we discuss in Section 2.3, these spins freely precess from the time of the second supernova up until merger. This precession somewhat scatters these angles, but leaves them with generally similar values, as seen in Figure 1.

Subpopulation 4: Alignment of secondary

We follow the standard mass-ratio model with effective tides presented in Gerosa et al. (2013), which assumes that after the first supernova, the secondary is realigned via tides or the CE prior to the second supernova. However, the primary BH is not realigned. Because the binary’s orbit shrinks during CE ejection, the kick velocity of the secondary is small relative to its pre-supernova velocity, causing the secondary to be only mildly misaligned (in general $\theta_{1}>\theta_{2}$).

We generate several thousand samples from each of these models. We plot samples from the four subpopulations in the $\{\cos{\theta_{1}},\cos{\theta_{2}}\}$ plane of the misalignment angles in Figure 1. From each of these models we randomly select $20$ mock detected systems (for a total of $N=80$) with component masses between $10M_{\odot}$ and $40M_{\odot}$ (cf. Abbott et al., 2016b); we describe the analysis of these mock GW signals in Section 3. We show the true values of the spin misalignment angles for the mock detected systems in Figure 2.

We use a \colorblackroughly astrophysical distribution of systems with sky positions and inclinations randomly chosen, and distances $D_{\mathrm{L}}$ distributed uniformly in volume, $p(D_{\mathrm{L}})\propto D_{\mathrm{L}}^{2}$, such that the distribution of signal-to-noise ratio (SNR) $\rho$ is $p(\rho)\propto\rho^{-4}$ (Schutz, 2011). \colorblackThis approximates the Universe as being static, with constant merger rates per unit time, at the distance scales probed by current detectors. We use a detection threshold (minimum network SNR) of $\rho_{\mathrm{min}}=12$ (Aasi et al., 2016; Berry et al., 2015).

2.3 Precession and spin–orbit resonances

After the second supernova, the evolution of the BBH is purely driven by relativistic effects and the orbit decays due to the emission of gravitational radiation (Peters & Mathews, 1963; Peters, 1964). As the BHs orbit, their spins precess around the total angular momentum (Apostolatos et al., 1994; Blanchet, 2014). In order to predict the spin misalignment angles when the frequency of GWs emitted by the binary are high enough (or equivalently when the orbital separation of the binary is sufficiently small) to be in the aLIGO band ($f_{\mathrm{GW}}>10~{}\mathrm{Hz}$), we take into account the post-Newtonian (PN) evolution of the spins by evolving the ten coupled differential equations given by Equations (14)–(17) in Gerosa et al. (2013). We begin our integrations at an orbital separation $a=1000M$, and integrate up until $f_{\mathrm{GW}}=10~{}\mathrm{Hz}$.222A more efficient method of evolving binaries from wide orbital separations to the frequencies where they enter the aLIGO band was introduced in Kesden et al. (2015) and Gerosa et al. (2015). This exploits the hierarchy of timescales in the problem and integrates precession averaged equations of motion on the radiation reaction timescale, rather than integrating the orbit-averaged equations we use here.

Some of these binaries are attracted to spin–orbit resonances (Schnittman, 2004). \colorblackIn particular, the binaries from subpopulation 4 are attracted to the $\Delta\Phi=\pm 180\degree$ resonance, where $\Delta\Phi$ is the angle between the projection of the two spins on the orbital plane. The \colorblackcurrent generation of ground-based GW observations are generally insensitive to this angle \colorblackfor binary black holes (Schmidt et al., 2015; Abbott et al., 2016a), \colorblackand the waveform model we use does not include it, so we focus on distinguishing subpopulations through the better-measured $\theta_{1}$ and $\theta_{2}$ angles.

3 GW parameter estimation

3.1 Signal analysis and inference

The strain measured by a GW detector is a combination of detector noise and (possibly) a GW signal $h(\boldsymbol{\Theta},t)$,

[TABLE]

Here $\boldsymbol{\Theta}$ is the vector of parameters describing the GW signal; for a general spinning circular BBH, there are $15$ parameters.333These parameters are (e.g., Veitch et al., 2015a): two component masses $\{m_{i}\}$; six spin parameters describing $\{\boldsymbol{S}_{i}\}$; two sky coordinates; distance $D_{\mathrm{L}}$; inclination and polarization angles; a reference time, and the orbital phase at this time. Given a data stream, we want to infer the most probable set of parameters for that data. To estimate the properties of the signal, waveform templates are matched to the data (Cutler & Flanagan, 1994; Veitch et al., 2015a; Abbott et al., 2016e).

The posterior probability for the parameters is given by Bayes’ theorem,

[TABLE]

where $p(d|\boldsymbol{\Theta})$ is the likelihood of observing the data given a choice of parameters, $p(\boldsymbol{\Theta})$ is the prior on those parameters, and the evidence $p(d)$ is a normalisation constant for the purposes of PE. The prior encodes our belief about the parameters before we considered the data: we assume that sources are uniformly distributed across the sky and in volume; that spin magnitudes are uniformly distributed between [math] and $1$; that spin orientations and the binary orientation are uniformly distributed across the surface of the sphere, and that component masses are uniformly distributed up to a maximum of $150M_{\odot}$ (cf. Abbott et al., 2016e). The likelihood is calculated from the residuals between the data and the signal template, assuming that the noise is Gaussian (Cutler & Flanagan, 1994):

[TABLE]

where the inner product $\left(g|h\right)$ is given by (Finn, 1992)

[TABLE]

and $S_{n}(f)$ is the (one-sided) noise power spectral density (Moore et al., 2015), which we take to be the design sensitivities for aLIGO and AdV respectively, with $f_{\mathrm{low}}=10~{}\mathrm{Hz}$ as is appropriate for the advanced detectors.

We sample the posterior distribution using the publicly available, Bayesian PE code LALInference (Veitch et al., 2015a).444Available as part of the LIGO Scientific Collaboration Algorithm Library (LAL) https://wiki.ligo.org/DASWG/LALSuite. For each event we obtain $\nu\sim 5000$ independent posterior samples. We show an example of the marginalised posterior distribution in $\{\cos{\theta_{1}},\cos{\theta_{2}}\}$ space for one of our 80 events in Figure 3. \colorblackUnless otherwise stated, we quote all parameters at a reference frequency $f_{\mathrm{ref}}=10~{}\mathrm{Hz}$.

We sample in the system frame (Farr et al., 2014b), where the binary is parametrised by the masses and spin magnitudes of the two component BHs $\{m_{i}\}$ and $\{\chi_{i}\}$, the spin misalignment angles $\{\theta_{i}\}$, the angle $\Delta\Phi$ between the projections of the two spins on the orbital plane, and the angle $\beta$ between the total and orbital angular momentum vectors. We find, in agreement with similar studies such as Littenberg et al. (2015) and Miller et al. (2015), that there is a strong preference for detecting GWs from nearly face-on binaries, since GW emission is strongest perpendicular to the orbital plane.

Following common practice in PE studies, we use a special realisation of Gaussian noise which is exactly zero in each frequency bin (Veitch et al., 2015b; Trifirò et al., 2016). Real GW detector noise will be non-Gaussian and non-stationary, and events will be recovered with non-zero noise.555Non-stationary, non-Gaussian noise has been shown not to affect average PE performance for binary neutron stars (Berry et al., 2015); however, these noise features could be more significant in analysing the shorter duration BBH signals. Non-zero noise-realisations will mean that in general the maximum likelihood parameters do not match the injection parameters; in the Gaussian limit, however, using a zero noise realisation is equivalent to averaging over a large number of random noise realisations, such that these offsets approximately cancel out (cf. Vallisneri, 2011). This assumption makes it straightforward to compare the posterior distributions, as differences only arise from the input parameters and not any the specific noise realisation.

3.2 Previous studies

Vitale et al. (2017) study GW measurements of BH spin misalignments in compact binaries containing at least one BH. They consider both BBHs (using IMRPhenomPv2 waveforms as we do here) and neutron star–black hole (NSBH) binaries (using inspiral-only SpinTaylorT4 waveforms). They fit a mixture model allowing for both a preferentially aligned/anti-aligned component and an isotropically misaligned component, excluding aligned/anti-aligned systems. They find that $\sim 100$ detections yield a $\sim 10\%$ precision for the measured aligned fraction. One of the main limitations of the analysis performed by Vitale et al. (2017) is that they only consider models which are mutually exclusive, \colorblackalthough this should not affect their results since the excluded region for their nonaligned model is negligible. Here, all of our formation models overlap in the parameter space of spin–orbit misalignment angles. \colorblackTherefore, we cannot directly apply the formalism of Vitale et al. (2017). The framework we develop here is able to correctly determine the relative contributions of multiple models, even when those models overlap in parameter space significantly, as expected in practice.

There have also been significant advances in the past few years in the understanding of PN spin–orbit resonances. These resonances occur when BH spins become aligned or anti-aligned with one another and precess in a common plane around the total angular momentum (Schnittman, 2004). This causes binaries to be attracted to different points in parameter space identified by $\Delta\Phi$, the angle between the projections of the two BH spins onto the orbital plane. Kesden et al. (2010) have shown that these resonances are effective at capturing binaries with mass ratios $0.4<q<1$ and spins $\chi_{i}>0.5$. For equal-mass binaries, \colorblackspin morphologies remain locked with binaries trapped in or out of resonance (Gerosa et al., 2017); however, it is unlikely for astrophysical formation scenarios to produce exactly equal mass binaries, although Marchant et al. (2016) predict nearly equal masses for the chemically homogeneous evolution channel.

Gerosa et al. (2013) show how the family of resonances that BBHs are attracted to can act as a diagnostic of the formation scenario for those binaries. Trifirò et al. (2016) demonstrate that GW measurements of spin misalignments can be used to distinguish between the two resonant families of $\Delta\Phi=0\degree$ and $\Delta\Phi=\pm 180\degree$. They use a full PE study to show that they can distinguish two families of PN resonances. However, they only consider a small corner of parameter space which contains binaries which will become locked in these PN resonances.

Our study extends on those discussed in several ways:

Rather than focusing on specific systems preferred in previous studies, we use injections from an astrophysically motivated population. Our injections have total masses $M=m_{1}+m_{2}$ in the range $10$–$40M_{\odot}$ and an astrophysical distribution of SNRs. 2. 2.

The misalignment angles of our BHs are given by simple but astrophysically motivated models introduced in Section 2. 3. 3.

For performing PE on individual GW events, we use the inspiral-merger-ringdown gravitational waveform IMRPhenomPv2 model, rather than the inspiral-only waveforms used in some of the earlier studies. 4. 4.

Most importantly, we perform a hierarchical Bayesian analysis on the posterior probability density functions of a mock catalog of detected events in order to make inferences about the underlying population.

4 Hierarchical analysis for population inference

PE on individual GW events yields samples from the posterior distributions for parameters under astrophysical prior constraints. We now wish to combine these individual measurements of BH spin misalignment angles in order to learn about the underlying population, which may act as a diagnostic for binary formation channels and binary evolution scenarios. Importantly, we are able to do this without reanalysing the data for the individual events.

Given a set of reasonable population synthesis model predictions for BBH spin misalignment angles, we would like to learn what mixture of those subpopulations best explains the observed data. Here we assume that the subpopulation distributions representing different formation channels are known perfectly, and use the same subpopulations that we drew our injections from to set these distributions. Thus, each subpopulation model $\Lambda_{\ell}$ ($\ell\in 1\dots 4$) corresponds to a known distribution of source parameters $p(\boldsymbol{\Theta}|\Lambda_{\ell})$. In practice, the uncertainty in the subpopulation models will be one of the challenges in carrying out accurate hierarchical inference.666The clustering approach of Mandel et al. (2015, 2017), which eschews assumptions about the subpopulation distributions, could provide an alternative pathway for robust but less informative inference on the data alone.

The overall mixture model is described by hyperparameters $\lambda_{\ell}$, corresponding to the fraction of each of the four subpopupulations, such that

[TABLE]

We assume that each event comes from one of these subpopulations:

[TABLE]

i.e., $\boldsymbol{\lambda}$ is a unit simplex.

For any individual event $\alpha$ ($\alpha=1,\ldots,N$) we have the posterior on $\boldsymbol{\Theta}$ given by

[TABLE]

where $p(\boldsymbol{\Theta})$ is the prior used by LALInference, $p(d_{\alpha})$ is the evidence (which is only a normalising factor in our analysis), and we represent $p(\boldsymbol{\Theta}|d_{\alpha})$ by a set of discrete samples $\{\boldsymbol{\Theta}_{i}^{k}\}$ where $k=1,\ldots,\nu_{\alpha}$.

We can write the likelihood for obtaining all of the events as the product over the individual likelihoods (Mandel, 2010; Hogg et al., 2010),

[TABLE]

where we have marginalised over the physical parameters of the individual events, and used Bayes’ theorem to obtain the final line. Since we have samples drawn from the posterior $p(\boldsymbol{\Theta}_{\alpha}|d_{\alpha})$, we can approximate posterior-weighted integrals (posterior averages) as a sum over samples (MacKay, 2003, chapter 29),

[TABLE]

where $f(\boldsymbol{\Theta})$ is some general function. Thus, we can rewrite Equation (14) as

[TABLE]

In effect, for each event we reweigh the evidence calculated using our general PE prior to what it would have been using a prior for the model of interest, and then combine these probabilities together to form a likelihood.

The posterior for $\boldsymbol{\lambda}$ is then

[TABLE]

for a choice of prior $p(\boldsymbol{\lambda})$. We assume a flat Dirichlet prior as shown in Figure 4. We sample from this posterior on $\boldsymbol{\lambda}$ using emcee (Foreman-Mackey et al., 2013), an affine-invariant ensemble sampler (Goodman & Weare, 2010).777Available from http://dan.iel.fm/emcee/.

5 Results

To gain a qualitative understanding of hierarchical modelling on the spin–orbit misalignment angles, we first consider inference under the assumption of perfect measurement accuracy for individual observations, and then introduce realistic measurement uncertainties. We then analyse the scaling of the inference accuracy with the number of observations.

5.1 Perfect measurement accuracy

Here, we assume that aLIGO–AdV GW observations could perfectly measure the spin–orbit misalignment angles of merging BBHs. In this case, the posterior is simply a delta function centered at the true value. Since our underlying astrophysical models have significant overlap in the $\{\cos(\theta_{1})$, $\cos(\theta_{2})\}$ plane, as shown in Figure 1, there is still ambiguity about which model a given event comes from.

We sample Equation (17), where our data consist of $80$ events with perfectly measured spin–orbit misalignments (as seen in Figure 2). This number of detections could be available by the end of the third observing run under optimistic assumptions about detector sensitivity improvements (Abbott et al., 2016b). The results of this analysis are shown in Figure 5.

We find that after $80$ BBH observations with perfect measurement accuracy, we would be able to confidently establish the presence of all four subpopulations. From this analysis, we can already understand some of the features of the posterior on the hyperparameters. For example, we see that there is a strong degeneracy between $\lambda_{1}$ and $\lambda_{3}$, since both of these models predict a large (nearly) aligned ($\theta_{1}=\theta_{2}=0$) population. There is a similar degeneracy between $\lambda_{2}$ and $\lambda_{4}$. We can also see that the fraction of exactly aligned systems ($\lambda_{1}$) and the fraction of systems with isotropically distributed spin–orbit misalignments ($\lambda_{2}$) are not strongly correlated. Both fractions are measured with to be between $\sim 0.15$ and $\sim 0.45$ at the $90\%$ credible level with $80$ BBH observations, corresponding to a fractional uncertainty of $\sim 50\%$.

5.2 Realistic measurement accuracy

We know that in practice GW detectors will not perfectly measure the spin–orbit misalignments of merging BBHs (see Figure 3 for a typical marginalised posterior). We now use the full set of $80$ LALInference posteriors, each containing $\sim 5000$ posterior samples as our input data, when sampling Equations (16) and (17).

We show the results of this analysis in Figure 6. Many of the features seen in the posteriors on the hyperparameters are the same as those seen in Figure 5, such as the strong anti-correlation between $\lambda_{1}$ and $\lambda_{3}$. We see that the posterior is not perfectly centred on the true $\lambda$ values, though the true values do have posterior support. While the hierarchical modelling unambiguously points to the presence of multiple subpopulations, with no single subpopulation able to explain the full set of observations, the data no longer require all four subpopulations to be present.

We have checked that the structure of this posterior is typical given the limited number of observations and the large measurement uncertainties. In the next section we show that our posteriors converge to the true values in the limit of a large number of detections.

5.3 Dependence on number of observations

The LALInference PE pipeline used to compute the posterior distributions for our 80 injections in Section 3 is computationally expensive. However, we would like to generate a larger catalogue of mock observations. First, this allows us to check that our analysis is self consistent by running many tests, such as confirming that the true result lies within the $P\%$ credible interval in $P\%$ of trials. Second, it allows us to predict how the accuracy of the inferred fractions of the subpopulations evolves as a function of the number of GW observations.

We develop approximations to these posteriors, similar to Mandel et al. (2017), based on the $80$ posterior distributions generated in Section 3. The best measured spin parameter is a combination of the two component spins called the effective inspiral spin $\chi_{\mathrm{eff}}\in[-1,1]$ (Ajith et al., 2011; Abbott et al., 2016e; Vitale et al., 2017):

[TABLE]

Having information about a single spin parameter makes it challenging to extract information about the spin distribution, but not impossible; for example, GW151226’s positive $\chi_{\mathrm{eff}}$ means that at least one spin must have non-zero magnitude and $\theta_{i}<90\degree$ (Abbott et al., 2016f).

To compute the approximate posteriors, we represent each observation with true parameter values $\chi_{\mathrm{eff}}^{\mathrm{true}}$ and $\cos{\theta_{1}}^{\mathrm{true}}$ by data which are maximum-likelihood estimates in a random noise realization:

[TABLE]

where $N(\mu,\sigma^{2})$ indicates a normal distribution. Posterior samples are then drawn using the same likelihood functions, centred on the maximum-likelihood data value,

[TABLE]

Here

[TABLE]

with $A=-0.2$, $B=0.3$ and $C=0.2$ based on a fit to our $80$ posteriors; the measurement uncertainty scales inversely with the signal-to-noise ratio $\rho$ (Cutler & Flanagan, 1994). In all cases, we only consider $\cos{\theta_{1}}$, $\cos{\theta_{2}}$ and $\chi_{\mathrm{eff}}$ in the permitted range of $[-1,1]$.

\color

blackThese posterior models are designed to mimic the typical properties of the posterior on the $\theta_{1}$ and $\theta_{2}$. We have chosen a particular strategy that allows us to make faithful models; starting with models for posteriors on $\chi_{\mathrm{eff}}$ and $\theta_{1}$ before converting them into posteriors on $\theta_{1}$ and $\theta_{2}$ is desirable because $\chi_{\mathrm{eff}}$ and $\theta_{1}$ are only weakly correlated, making it possible to write down independent expressions for the individual likelihoods, while $\theta_{1}$ and $\theta_{2}$ are strongly correlated, so would require a more complex joint likelihood model. The posteriors already implicitly account for marginalization over correlated parameters such as spin magnitudes and the mass ratio, which will not be known in practice. These mock posteriors will not correctly reproduce actual posteriors for specific events, since they do not include the noise realization, except probabilistically; they only reproduce the average properties of the posterior distributions.

We draw $\nu=5000$ posterior samples of $\chi_{\mathrm{eff}}$ and $\cos{\theta_{1}}$ independently, and calculate the values of $\cos{\theta_{2}}$ using Equation (18), fixing the mass ratio $q$ and spin magnitudes $\chi_{i}$ to their true values. This builds the correct degeneracies between $\cos{\theta_{1}}$ and $\cos{\theta_{2}}$ into the mock posteriors. We show an example of a posterior distribution generated this way in Figure 7.

Using this method, we generate spin–orbit misalignment posteriors for $400$ BBHs drawn in equal fractions ($\lambda_{i}=0.25$) from the four subpopulation models introduced in Section 2. Using the method introduced in Section 4, we calculate the posteriors on the $\lambda$ parameters after [math] (prior), $10$, $20$, $40$, $80$, $160$ and $400$ observations, similar to Mandel et al. (2017). In Figure 8 we show the $68\%$ and $95\%$ credible intervals for the fraction $\lambda_{2}$ of observed BBHs coming from an isotropic distribution (subpopulation 2) corresponding to dynamically formed binaries. Given our models and incorporating realistic measurement uncertainties, we find that this fraction can be measured with a $\sim 40\%$ fractional uncertainty after $100$ observations. Since subpopulation 1 and 3 are somewhat degenerate in our model, we find that the combined fraction $\lambda_{1}+\lambda_{3}$ is a well measured parameter (as shown in Figure 9), whilst the individual components are measured less well. For $N\gtrsim 100$ observations the uncertainties in the $\lambda_{\ell}$ scale as the inverse square root of the number of observations; e.g., $\sigma_{\lambda_{1}+\lambda_{3}}\approx 0.8\,N^{-1/2}$.

Although $\gtrsim 100$ observations are required to accurately measure the contribution of each of the four subpopulations, it is possible to test for more extreme models with fewer observations. For example, $\sim 20$ observations are sufficient to demonstrate the presence of an isotropic subpopulation at the 95% credible level.

Even fewer observations are needed to confidently rule out the hypothesis that all observations come from the exactly aligned or isotropic subpopulations. We draw observations from the isotropic subpopulation and calculate the ratio of the evidence (Bayes factor) $Z_{\mathrm{aligned}}$ for the model under which all BBH spins are exactly aligned ($\lambda_{1}=1$) to the evidence $Z_{\mathrm{isotropic}}$ for the model under which all BBH spins are isotropically distributed ($\lambda_{2}=1$):

[TABLE]

\color

blackInference using small numbers of observations is sensitive to the exact choice of these observations (both true parameters and measurement errors), which we randomly draw from the relevant distributions. We \colorblacktherefore repeat this test $100$ times to account for the random nature of the mock catalogue. In all cases, we find that with only $5$ observations of BBHs with component spin magnitudes $\chi=0.7$, the exactly aligned model $\lambda_{1}=1$ can be ruled out at more than $5~{}\sigma$ confidence. Similarly, when drawing from the exactly aligned model, we find that the hypothesis that all events come from an isotropic population $\lambda_{2}=1$ can be ruled out at more than $5~{}\sigma$ confidence in all tests with $5$ observations.

6 Discussion and conclusions

With the first direct observations of GWs from merging BBHs, the era of GW astronomy has begun. GW observations provide a new and unique insight into the properties of BBHs and their progenitors. For individual systems, we can infer the masses and spins of the component BHs; combining these measurements we can learn about the population, and place constraints on the formation mechanisms for these systems, whether as the end point of isolated binary evolution or as the results of dynamical interactions.

In this work, we investigated how measurements of BBH spin–orbit misalignments could inform our understanding of the BBH population. We chose the properties of our sources to match those we hope to observe with aLIGO and AdV (at design sensitivity), using four different astrophysically motivated subpopulations for the distribution of spin–orbit misalignment angles, each reflecting a different formation scenario. We performed a Bayesian analysis of GW signals (using full inspiral–merger–ringdown waveforms) for a population of BBHs. We assumed a mixture model for the overall population of BBH spin–orbit misalignments and combined the full PE results from our GW analysis in a hierarchical framework to infer the fraction of the population coming from each subpopulation. A similar analysis could be performed following the detection of real signals.

\color

blackAdopting a population with spins of $\chi=0.7$, we demonstrate that the fraction of BBHs with spins preferentially aligned with the orbital angular momentum ($\lambda_{1}+\lambda_{3}$) is well measured and can be measured with an uncertainty of $\sim 10\%$ with $100$ observations, scaling as the inverse square root of the number of observations. We also show that after $100$ observations, we can measure the fraction $\lambda_{2}$ of the subpopulation with isotropic spins (assumed to correspond to dynamical formation) with a fractional uncertainty of $\sim 40\%$. Extreme hypotheses can be tested and ruled out with even fewer observations. For example, with just $5$ observations we can rule out the hypothesis that all BBHs have their spins exactly aligned with high confidence ($>5~{}\sigma$) if the true population has isotropically distributed spins, and vice versa. This number of observations may be reached by the end of the second aLIGO observing run.

One limitation of the current approach is the assumption that the subpopulation distributions are known perfectly. This will not be the case in practice, but the simplified models considered here are still relevant as parametrizable proxies for astrophysical scenarios. Hierarchical modelling with strong population assumptions could lead to systematic biases in the interpretation of the observations if those assumptions are not representative of the true populations; this can be mitigated by coupling such analysis with weakly modelled approaches, such as observation-based clustering (Mandel et al., 2015, 2017).

In this work we have not taken into account observational selection effects, \colorblackparticularly the differences in the detectabilty of different subpopulations because their different spin parameter distributions impact the SNR (Reisswig et al., 2009). These must be incorporated in the analysis \colorblackto correctly infer the intrinsic subpopulation fractions (Farr et al., 2015; Mandel et al., 2016), but the impact is not expected to be significant in this case. Care must also be taken to avoid biases when performing an hierarchical analysis with real observations, since the observations will not be drawn from the same distribution as the priors used for the analysis of individual events. Our framework accounts for the differences in the priors on the parameters of interest (spin–orbit misalignment angles) between the original PE and model predictions, but not for any discrepancy in the priors of the parameters we marginalize over (e.g., masses); this is likely a second-order effect.

Neither theoretical models nor observations can currently place tight constraints on the spin magnitudes of BHs. We therefore chose to give all BHs a spin magnitude of $0.7$ in this study. This choice is clearly ad hoc; we expect a distribution of BH spin magnitudes in nature. Since GW events with BHs with low spin will not constrain the spin–orbit misalignment angles well, a distribution of spin magnitudes containing lower BH spins will act to increase the requirements for the numbers of observations quoted here. \colorblackEven with significantly smaller BH spin magnitudes, a few to a few tens of observations can distinguish between sufficiently different subpopulation models for spin–orbit misalignment (perfectly aligned versus isotropic) as shown in Farr et al. (2017).

Here we have assumed that BHs receive large natal kicks, comparable to neutron stars, leading to relatively large spin–orbit misalignments even for isolated binary evolution. We further assume that the effect of the kick is simply to tilt the orbital plane, but not the BH spin. There is, however, evidence from the Galactic double pulsar PSR J0737$-$3039 that the second born pulsar received a spin tilt at birth (Farr et al., 2011). \colorblackFurthermore, alignment through mass transfer prior to BH formation may be imperfect. Our models of BBHs that are preferentially but imperfectly aligned because of high natal kicks can be viewed as proxies for misalignment through a combination of these effects.

Optimal hierarchical modelling should fold in all available information, including component masses \colorblack(cf. Rodriguez et al., 2016) and spin magnitudes \colorblack(cf. Gerosa & Berti, 2017; Fishbach et al., 2017), into a single analysis. Complementary electromagnetic observations of high-mass X-ray binaries, Galactic radio pulsars, short gamma ray bursts, supernovae and luminous red novae will contribute to a concordance model of massive binary formation and evolution.

Acknowledgements

This work was supported in part by the Science and Technology Facilities Council. IM acknowledges support from the Leverhulme Trust. We are grateful for computational resources provided by Cardiff University, and funded by an STFC grant supporting UK Involvement in the Operation of Advanced LIGO. SS and IM acknowledge partial support by the National Science Foundation under Grant No. NSF PHY11-25915. We are grateful to colleagues from the Institute of Gravitational-wave Astronomy at the University of Birmingham, as well as Ben Farr, M. Coleman Miller, Chris Pankow and Salvatore Vitale for fruitful discussions. We thank Davide Gerosa and the anonymous referee for comments on the paper. Corner plots were made using triangle.py available from https://github.com/dfm/triangle.py. This is LIGO Document P1700007.