Hierarchical analysis of gravitational-wave measurements of binary black hole spin-orbit misalignments
MMNRAS , 1–12 (2017) Preprint 13 September 2017 Compiled using MNRAS L A TEX style file v3.0
Hierarchical analysis of gravitational-wave measurementsof binary black hole spin–orbit misalignments
Simon Stevenson, (cid:63) Christopher P. L. Berry, and Ilya Mandel School of Physics and Astronomy and Institute of Gravitational Wave Astronomy, University of Birmingham, Edgbaston, BirminghamB15 2TT, United Kingdom
13 September 2017
ABSTRACT
Binary black holes may form both through isolated binary evolution and throughdynamical interactions in dense stellar environments. The formation channel leaves animprint on the alignment between the black hole spins and the orbital angular momen-tum. Gravitational waves from these systems directly encode information about thespin–orbit misalignment angles, allowing them to be (weakly) constrained. Identifyingsub-populations of spinning binary black holes will inform us about compact binaryformation and evolution. We simulate a mixed population of binary black holes withspin–orbit misalignments modelled under a range of assumptions. We then developa hierarchical analysis and apply it to mock gravitational-wave observations of thesepopulations. Assuming a population with dimensionless spin magnitudes of χ = 0 . ∼ ∼ Key words: black hole physics – gravitational waves – methods: data analysis –stars: evolution – stars: black hole
Compact binaries containing two stellar-mass black holes(BHs) 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 (BBHs) are apromising source of gravitational waves (GWs) for ground-based detectors such as Advanced LIGO (aLIGO; Aasi et al.2015), Advanced Virgo (AdV; Acernese et al. 2015) andKAGRA (Aso et al. 2013). Searches of data from the firstobserving run (Aasi et al. 2016) of aLIGO yielded threelikely BBH coalescences: GW150914 (Abbott et al. 2016d),GW151226 (Abbott et al. 2016f) and LVT151012 (Abbottet al. 2016c,b). A further BBH coalescence, GW170104, hasbeen reported from the on-going second observing run (Ab-bott et al. 2017). GW observations give a unique insight intothe properties of BBHs. We will examine one of the ways in (cid:63)
E-mail: [email protected] which black hole spin measurements can be used to constrainformation mechanisms.GW observations inform our understanding of BBHevolution in two ways: from the merger rate, and from theproperties of the individual systems. The merger rate ofBBHs is inferred from the number of detections; it is un-certain as a consequence of the small number of BBH ob-servations so far. Currently, merger rates are estimated tobe 12–213 Gpc − yr − (Abbott et al. 2016h; Abbott et al.2017). These rates are broadly consistent with predictionsfrom both population synthesis models of isolated binaryevolution (e.g., Voss & Tauris 2003; Lipunov et al. 2009; Do-minik et al. 2012; Postnov & Yungelson 2014) and dynam-ical formation models (e.g., Sigurdsson & Hernquist 1993;Sadowski et al. 2008; Rodriguez et al. 2015). Possible pro-genitors systems of BBHs, including Cyg X-3 (Belczynskiet al. 2013), IC 10 X-1 (Bulik et al. 2011) and NGC 300 X-1 (Crowther et al. 2010) provide some additional limits onBBH merger rates, but extrapolation is hindered by currentobservational uncertainties.The parameters of individual systems can be estimated c (cid:13) a r X i v : . [ a s t r o - ph . H E ] S e p S. Stevenson et al. by comparing the measured GW signal with template wave-forms (Abbott et al. 2016e). The masses and spins of theBHs can be measured through their influence on the inspi-ral, merger and ringdown of the system (Cutler & Flanagan1994; Poisson & Will 1995; Finn & Chernoff 1993). The dis-tribution of parameters observed by aLIGO will encode in-formation about the population of BBHs, and may also helpto shed light on their formation channels (O’Shaughnessyet 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, withobservations suggesting their dimensionless spin parameters χ take the full range of allowed values between 0 and 1 (Mc-Clintock et al. 2011; Fragos & McClintock 2015; Miller &Miller 2015). Stars formed in binaries are expected to havetheir rotational axis aligned with the orbital angular mo-mentum (e.g., Boss 1988; Albrecht et al. 2007), althoughthere is observational evidence this is not always the case(e.g., Albrecht et al. 2009, 2014). Even if binaries are bornwith misaligned spins, there are many processes in binaryevolution which can act to align the spin of stars, such asrealignment during a stable mass accretion phase (Bardeen& Petterson 1975; Kalogera 2000; King et al. 2005), accre-tion onto a BH passing through a common-envelope (CE)event (Ivanova et al. 2013), and realignment through tidalinteractions in close binaries (e.g., Zahn 1977; Hut 1981).On the other hand, asymmetric mass loss during super-nova explosions can tilt the orbital plane in binaries (Brandt& Podsiadlowski 1995; Kalogera 2000), leading to BH spinsbeing misaligned with respect to the orbital angular mo-mentum vector. Population synthesis studies of BH X-raybinaries predict that these misalignments are generally small(Kalogera 2000), with Fragos et al. (2010) finding that theprimary BH is typically misaligned by (cid:46) ◦ . However, elec-tromagnetic observations of high mass X-ray binaries con-taining BHs have hinted that the BHs may be more signifi-cantly misaligned (Miller & Miller 2015). One such system isthe microquasar V4641 Sgr (Orosz et al. 2001; Martin et al.2008) where the primary BH is has been interpreted to bemisaligned by > ◦ .Alternatively, BBHs can form dynamically in dense stel-lar environments such as globular clusters. In these environ-ments, it is expected that the distribution of BBH spin–orbitmisalignment angles is isotropic (e.g., Rodriguez et al. 2016).The distribution of BBH spin–orbit misalignments thereforecontains information about their formation mechanisms.Constraints on spin alignment from GW observationsso far are weak (Abbott et al. 2016e; Abbott et al. 2016a;Abbott et al. 2016b, 2017). Some configurations, such asanti-aligned spins for GW151226 (Abbott et al. 2016f), aredisfavoured; however, there is considerable uncertainty inthe spin magnitude and orientation. Determining the spinsprecisely is difficult because their effects on the waveform canbe intrinsically small (especially if the the source is viewedface on), and because of degeneracies between the spin andmass parameters (Poisson & Will 1995; Baird et al. 2013;Farr et al. 2016). Although the spins of individual systemsare difficult to measure, here we show it is possible to useinferences from multiple systems to build a statistical modelfor the population (cf. Vitale et al. 2017).This paper describes how to combine posterior probabil- ity density functions on spin–orbit misalignment angles frommultiple GW events to explore the underlying population.We develop a hierarchical analysis in order to combinemultiple GW observations of BBH spin–orbit misalignmentsto give constraints on the fractions of BBHs forming throughdifferent channels. We consider different populations of po-tential spin–orbit misalignments, each representing differentassumptions about binary formation, and use the GW ob-servations to infer the fraction of binaries from each popula-tion. In the field of exoplanets, similar hierarchical analyseshave been used to make inference on the frequency of Earth-like exoplanets from measurements of the period and radiusof individual exoplanet candidates (e.g., Foreman-Mackeyet al. 2014; Farr et al. 2014a). Other examples of the useof hierarchical analyses in astrophysics include modelling apopulation of trans-Neptunian objects (Loredo 2004), mea-surements of spin–orbit misalignments in exoplanets (Naozet al. 2012), measurements of the eccentricity distribution ofexoplanets (Hogg et al. 2010) and the measurement of themass distribution of galaxy clusters (Lieu et al. 2017).In Section 2 we introduce our simplified population syn-thesis models for BBHs, paying special attention to the BHspins. We briefly describe in Section 3 the parameter estima-tion (PE) pipeline that will be employed to infer the prop-erties (such as misalignment angles) of real GW events, anddiscuss previous spin-misalignment studies in the literature.We introduce a framework for combining posterior proba-bility density functions on spin–orbit misalignment anglesfrom multiple GW events to explore the underlying popula-tion in Section 4. We demonstrate the method using a setof mock GW events in Section 5, and show that tens of ob-servations will be sufficient to distinguish subpopulations ofcoalescing binary black holes, assuming spin magnitudes of ∼ .
7. Lower typical BH spin magnitudes would reduce theaccuracy with which the spin-orbit misalignment angle canbe measured, therefore requiring more observations to ex-tract information about subpopulations. We also show thatmore extreme models, such as the hypothesis that all BBHshave their spins exactly aligned with the orbital angular mo-mentum, can be ruled out at a 5 σ confidence level withonly O (5) observations of rapidly spinning BBHs. Finally,we conclude and suggest areas which require further studyin Section 6. Owing to the many uncertainties pertaining to stellar spinsand their evolution in a binary, and the fact that keepingtrack of stellar spin vectors can be computationally inten-sive, many population synthesis models choose not to includespin evolution. However, the distribution of spins of the finalmerging BHs is one of the observables that can be measuredwith the advanced GW detectors. In this section, we there-fore implement a simplified population synthesis model toevolve an ensemble of binaries that will be detectable withaLIGO and AdV and predict their distributions of spin–orbitmisalignments. We describe the assumed mass distribution,spin distribution, and spin evolution below.
MNRAS000
MNRAS000 , 1–12 (2017)
BH spin misalignments We assume the same simplified mass distribution for all ofour models, so that any differences in the final spin distri-butions are purely due to our assumptions about the spin–orbit misalignments described in the next section. There aremany uncertainties in the evolution of isolated massive bina-ries, including (but not limited to) uncertainties in the initialdistributions of the orbital elements (de Mink & Belczynski2015), the strength of stellar winds in massive stars (Bel-czynski et al. 2010), the effect of rotation of massive starson stellar evolution (de Mink et al. 2013; Ram´ırez-Agudeloet al. 2015), the natal kicks (if any) given to BHs (Repettoet al. 2012; Mandel 2016; Mirabel 2016; Repetto et al. 2017)and the efficiency of the CE (Ivanova et al. 2013; Kruckowet al. 2016). Population synthesis methods are large Monte-Carlo simulations using semi-analytic prescriptions in orderto explore the effect these uncertainties have on the pre-dicted distributions of compact binaries. Instead, we adopt anumber of simplifications that allow us to produce an astro-physically plausible distribution which should not, however,be considered representative of the actual mass distributionof BBHs.We simulate massive binaries with semimajor axis a drawn from a distribution uniform in ln a (Abt 1983). Thecomponents of the binary are a massive primary BH withmass m and a secondary star at the end of its main sequencelifetime.The primary BH was formed from a massive star withzero-age main sequence (ZAMS) mass m ZAMS1 drawn fromthe initial-mass function (IMF) with a power law index of − .
35 (Salpeter 1955; Kroupa 2001). The mass ratio of thebinary at ZAMS q ZAMS is drawn from a flat distribution[0 , m ZAMS2 = q ZAMS m ZAMS1 .We calculate the final remnant mass m i as a function ofthe ZAMS mass m ZAMS i for each star using a fit to Figure 12in Woosley et al. (2002). For stars with 30 < m ZAMS i /M (cid:12) <
50, in which range BHs are formed after some delay by fall-back of ejecta, we use m i = 30 (cid:18) m ZAMS i M (cid:12) (cid:19) α M (cid:12) , (1)with α = 3 .
9. For more massive stars with m ZAMS i > M (cid:12) ,which are massive enough to directly collapse during theiron-core collapse to form BHs, we use m i = 0 . m ZAMS i . (2)We only consider BBHs with component masses above10 M (cid:12) below, consistent with aLIGO detections to date (Ab-bott et al. 2016b, 2017), and hence omit stars with ZAMSmasses below 30 M (cid:12) from our population.We assume that the binary has negligible eccentricity e = 0, appropriate for post-CE systems. In all models wehave assumed that both main-sequence stars in the binaryare born with their rotation axis aligned with the orbitalangular momentum axis. In general, the first supernova willmisalign the spins due to any natal kick imparted on theremnant. There are expected to be mass-transfer phases be-tween the first and second supernovae which may realignboth the spins of the primary BH and the secondary BHprogenitor; we vary the assumed degree of realignment inour models. We assume that BHs receive natal kicks comparable tothose received by neutron stars (Hobbs et al. 2005), namelydrawn from a Maxwellian with a root-mean-square velocityof ∼
250 km s − . This assumption will lead to the maximumamount of spin misalignment, and may be consistent withneutrino-driven kicks; if the natal kicks are due to asym-metric ejection of baryonic matter, then any fall-back (Fryeret al. 2012) onto BHs during formation will reduce the kickmagnitude and thus the spin misalignment.BH spins magnitudes can take any value 0 (cid:54) χ i <
1, butwe set χ i = 0 . We model the overall population of BBHs as a mixture of 4subpopulations, each of which makes differing assumptionsleading to distinct spin–orbit misalignment distributions.We define the spin–orbit misalignment angle as the an-gle between the spin vector ˆ S i of binary component i ∈ , L ,cos θ i = ˆ S i · ˆ L , (3)where S i = χ i m i ˆ S i , (4)and m i is the component mass ( m (cid:62) m ). We will considerhow a set of spin-misalignment measurements could be usedto infer BBH formation mechanisms.
Subpopulation 1: Exactly aligned
We assume that ir-respective of all prior processes, both BHs have their spins Throughout this paper we use geometric units G = c = 1 unlessotherwise stated.MNRAS , 1–12 (2017) S. Stevenson et al. aligned with the orbital angular momentum vector after thesecond supernova, such that cos θ = cos θ = 1. This maybe the case if BHs receive no kicks. GW searches often as-sume BBHs have aligned spins as this simplification makesthe search less computationally demanding (Abbott et al.2016c). Subpopulation 2: Isotropic/dynamical formation
We assume that BBHs are formed dynamically, such thatthe distribution of spin angles is isotropic. Initially isotropicdistributions of spins remain isotropic (Schnittman 2004).We still generate the binary mass distribution with ourstandard approach, so that the only difference in BBHsbetween this model and the others is the spin distribution.
Subpopulation 3: Alignment before second SN
Motivated by Kalogera (2000), we assume that the spinsof both components are aligned with the orbital angularmomentum after a CE event and prior to the secondsupernova. The tilt of the orbital plane caused by thesecond supernova is then taken to be the spin misalignmentangle of both components, i.e. cos θ = cos θ . As we discussin Section 2.3, these spins freely precess from the timeof the second supernova up until merger. This precessionsomewhat scatters these angles, but leaves them withgenerally similar values, as seen in Figure 1. Subpopulation 4: Alignment of secondary
We fol-low the standard mass-ratio model with effective tidespresented in Gerosa et al. (2013), which assumes that afterthe first supernova, the secondary is realigned via tides orthe CE prior to the second supernova. However, the primaryBH is not realigned. Because the binary’s orbit shrinksduring CE ejection, the kick velocity of the secondary issmall relative to its pre-supernova velocity, causing thesecondary to be only mildly misaligned (in general θ > θ ).We generate several thousand samples from each ofthese models. We plot samples from the four subpopula-tions in the { cos θ , cos θ } plane of the misalignment anglesin Figure 1. From each of these models we randomly select20 mock detected systems (for a total of N = 80) with com-ponent masses between 10 M (cid:12) and 40 M (cid:12) (cf. Abbott et al.2016b); we describe the analysis of these mock GW signals inSection 3. We show the true values of the spin misalignmentangles for the mock detected systems in Figure 2.We use a roughly astrophysical distribution of systemswith sky positions and inclinations randomly chosen, anddistances D L distributed uniformly in volume, p ( D L ) ∝ D ,such that the distribution of signal-to-noise ratio (SNR) ρ is p ( ρ ) ∝ ρ − (Schutz 2011). This approximates the Universeas being static, with constant merger rates per unit time, atthe distance scales probed by current detectors. We use adetection threshold (minimum network SNR) of ρ min = 12(Aasi et al. 2016; Berry et al. 2015). After the second supernova, the evolution of the BBH ispurely driven by relativistic effects and the orbit decays dueto the emission of gravitational radiation (Peters & Math-ews 1963; Peters 1964). As the BHs orbit, their spins pre-cess around the total angular momentum (Apostolatos et al.1994; Blanchet 2014). In order to predict the spin misalign-ment angles when the frequency of GWs emitted by the binary are high enough (or equivalently when the orbitalseparation of the binary is sufficiently small) to be in theaLIGO band ( f GW >
10 Hz), we take into account thepost-Newtonian (PN) evolution of the spins by evolving theten coupled differential equations given by Equations (14)–(17) in Gerosa et al. (2013). We begin our integrations atan orbital separation a = 1000 M , and integrate up until f GW = 10 Hz. Some of these binaries are attracted to spin–orbit reso-nances (Schnittman 2004). In particular, the binaries fromsubpopulation 4 are attracted to the ∆Φ = ± ◦ reso-nance, where ∆Φ is the angle between the projection of thetwo spins on the orbital plane. The current generation ofground-based GW observations are generally insensitive tothis angle for binary black holes (Schmidt et al. 2015; Ab-bott et al. 2016a), and the waveform model we use doesnot include it, so we focus on distinguishing subpopulationsthrough the better-measured θ and θ angles. The strain measured by a GW detector is a combination ofdetector noise and (possibly) a GW signal h ( Θ , t ), d ( t ) = n ( t ) + h ( Θ , t ) . (5)Here Θ is the vector of parameters describing the GW sig-nal; for a general spinning circular BBH, there are 15 pa-rameters. Given a data stream, we want to infer the mostprobable set of parameters for that data. To estimate theproperties of the signal, waveform templates are matchedto the data (Cutler & Flanagan 1994; Veitch et al. 2015a;Abbott et al. 2016e).The posterior probability for the parameters is given byBayes’ theorem, p ( Θ | d ) = p ( d | Θ ) p ( Θ ) p ( d ) , (6)where p ( d | Θ ) is the likelihood of observing the data given achoice of parameters, p ( Θ ) is the prior on those parameters,and the evidence p ( d ) is a normalisation constant for the pur-poses of PE. The prior encodes our belief about the param-eters before we considered the data: we assume that sourcesare uniformly distributed across the sky and in volume; thatspin magnitudes are uniformly distributed between 0 and 1;that spin orientations and the binary orientation are uni-formly distributed across the surface of the sphere, and thatcomponent masses are uniformly distributed up to a maxi-mum of 150 M (cid:12) (cf. Abbott et al. 2016e). The likelihood is A more efficient method of evolving binaries from wide orbitalseparations to the frequencies where they enter the aLIGO bandwas introduced in Kesden et al. (2015) and Gerosa et al. (2015).This exploits the hierarchy of timescales in the problem and in-tegrates precession averaged equations of motion on the radia-tion reaction timescale, rather than integrating the orbit-averagedequations we use here. These parameters are (e.g., Veitch et al. 2015a): two compo-nent masses { m i } ; six spin parameters describing { S i } ; two skycoordinates; distance D L ; inclination and polarization angles; areference time, and the orbital phase at this time.MNRAS000
10 Hz), we take into account thepost-Newtonian (PN) evolution of the spins by evolving theten coupled differential equations given by Equations (14)–(17) in Gerosa et al. (2013). We begin our integrations atan orbital separation a = 1000 M , and integrate up until f GW = 10 Hz. Some of these binaries are attracted to spin–orbit reso-nances (Schnittman 2004). In particular, the binaries fromsubpopulation 4 are attracted to the ∆Φ = ± ◦ reso-nance, where ∆Φ is the angle between the projection of thetwo spins on the orbital plane. The current generation ofground-based GW observations are generally insensitive tothis angle for binary black holes (Schmidt et al. 2015; Ab-bott et al. 2016a), and the waveform model we use doesnot include it, so we focus on distinguishing subpopulationsthrough the better-measured θ and θ angles. The strain measured by a GW detector is a combination ofdetector noise and (possibly) a GW signal h ( Θ , t ), d ( t ) = n ( t ) + h ( Θ , t ) . (5)Here Θ is the vector of parameters describing the GW sig-nal; for a general spinning circular BBH, there are 15 pa-rameters. Given a data stream, we want to infer the mostprobable set of parameters for that data. To estimate theproperties of the signal, waveform templates are matchedto the data (Cutler & Flanagan 1994; Veitch et al. 2015a;Abbott et al. 2016e).The posterior probability for the parameters is given byBayes’ theorem, p ( Θ | d ) = p ( d | Θ ) p ( Θ ) p ( d ) , (6)where p ( d | Θ ) is the likelihood of observing the data given achoice of parameters, p ( Θ ) is the prior on those parameters,and the evidence p ( d ) is a normalisation constant for the pur-poses of PE. The prior encodes our belief about the param-eters before we considered the data: we assume that sourcesare uniformly distributed across the sky and in volume; thatspin magnitudes are uniformly distributed between 0 and 1;that spin orientations and the binary orientation are uni-formly distributed across the surface of the sphere, and thatcomponent masses are uniformly distributed up to a maxi-mum of 150 M (cid:12) (cf. Abbott et al. 2016e). The likelihood is A more efficient method of evolving binaries from wide orbitalseparations to the frequencies where they enter the aLIGO bandwas introduced in Kesden et al. (2015) and Gerosa et al. (2015).This exploits the hierarchy of timescales in the problem and in-tegrates precession averaged equations of motion on the radia-tion reaction timescale, rather than integrating the orbit-averagedequations we use here. These parameters are (e.g., Veitch et al. 2015a): two compo-nent masses { m i } ; six spin parameters describing { S i } ; two skycoordinates; distance D L ; inclination and polarization angles; areference time, and the orbital phase at this time.MNRAS000 , 1–12 (2017) BH spin misalignments cos( θ ) c o s ( θ ) cos( θ ) c o s ( θ ) cos( θ ) c o s ( θ ) Figure 1.
Three of the four astrophysically motivated subpopulations making up our mixture model for BBH spin misalignment angles θ and θ described in Section 2. Subpopulation 1 (not shown) has both spins perfectly aligned (cos θ = cos θ = 1), so all pointswould lie in the top right corner. In subpopulation 2 (left), both spins are drawn from an isotropic distribution, and so the samples aredistributed uniformly in the plane. In subpopulation 3 BH spins are aligned with the orbital angular momentum just prior to the secondsupernova. In subpopulation 4, the secondary BH has its spin aligned with the orbital angular momentum prior to the second supernova,whilst the primary is misaligned. See Section 2 for more details. Spins are quoted at a GW frequency of f ref = 10 Hz. Figure 2.
True values of BH spin–orbit misalignment anglescos θ and cos θ for a mixture of 20 draws from each of our foursubpopulations. Exactly aligned systems from subpopulation 1 sitin the upper right corner of this diagram and thus are not shown.Systems drawn from subpopulation 2 are shown as blue crosses,those from subpopulation 3 as red squares and those from sub-population 4 as green triangles. The injection plotted in Figures 3and 7 is circled in magenta. calculated from the residuals between the data and the sig-nal template, assuming that the noise is Gaussian (Cutler &Flanagan 1994): p ( d | Θ ) ∝ exp (cid:20) −
12 ( d − h ( Θ ) | d − h ( Θ )) (cid:21) , (7)where the inner product ( g | h ) is given by (Finn 1992)( g | h ) = 4 (cid:60) (cid:90) ∞ f low ˜ g ( f )˜ h ∗ ( f ) S n ( f ) d f, (8)and S n ( f ) is the (one-sided) noise power spectral density(Moore et al. 2015), which we take to be the design sensi-tivities for aLIGO and AdV respectively, with f low = 10 Hzas is appropriate for the advanced detectors.We sample the posterior distribution using the publiclyavailable, Bayesian PE code LALInference (Veitch et al.2015a). For each event we obtain ν ∼ Available as part of the LIGO Scientific Collaboration Algo-rithm Library (LAL) https://wiki.ligo.org/DASWG/LALSuite.
Figure 3.
Marginalised posterior samples for one of the 80 eventsshown in Figure 2, generated by analysing mock GW data us-ing
LALInference . The true spin–orbit misalignments (thick blackplus) for this event were cos θ = − .
08 and cos θ = 0 .
35, with anetwork SNR of 15 .
35. The dashed black diagonal line shows theline of constant χ eff . posterior distribution in { cos θ , cos θ } space for one of our80 events in Figure 3. Unless otherwise stated, we quote allparameters at a reference frequency f ref = 10 Hz.We sample in the system frame (Farr et al. 2014b),where the binary is parametrised by the masses and spinmagnitudes of the two component BHs { m i } and { χ i } , thespin misalignment angles { θ i } , the angle ∆Φ between theprojections of the two spins on the orbital plane, and theangle β between the total and orbital angular momentumvectors. We find, in agreement with similar studies such asLittenberg et al. (2015) and Miller et al. (2015), that thereis a strong preference for detecting GWs from nearly face-on binaries, since GW emission is strongest perpendicularto the orbital plane.Following common practice in PE studies, we use aspecial realisation of Gaussian noise which is exactly zeroin each frequency bin (Veitch et al. 2015b; Trifir`o et al.2016). Real GW detector noise will be non-Gaussian andnon-stationary, and events will be recovered with non-zeronoise. Non-zero noise-realisations will mean that in general Non-stationary, non-Gaussian noise has been shown not to af-MNRAS , 1–12 (2017)
S. Stevenson et al. the maximum likelihood parameters do not match the in-jection parameters; in the Gaussian limit, however, using azero noise realisation is equivalent to averaging over a largenumber of random noise realisations, such that these off-sets approximately cancel out (cf. Vallisneri 2011). This as-sumption makes it straightforward to compare the posteriordistributions, as differences only arise from the input param-eters and not any the specific noise realisation.
Vitale et al. (2017) study GW measurements of BH spinmisalignments in compact binaries containing at least oneBH. They consider both BBHs (using
IMRPhenomPv2 wave-forms as we do here) and neutron star–black hole (NSBH) bi-naries (using inspiral-only
SpinTaylorT4 waveforms). Theyfit a mixture model allowing for both a preferentiallyaligned/anti-aligned component and an isotropically mis-aligned component, excluding aligned/anti-aligned systems.They find that ∼
100 detections yield a ∼
10% precision forthe measured aligned fraction. One of the main limitationsof the analysis performed by Vitale et al. (2017) is that theyonly consider models which are mutually exclusive, althoughthis should not affect their results since the excluded regionfor their nonaligned model is negligible. Here, all of our for-mation models overlap in the parameter space of spin–orbitmisalignment angles. Therefore, we cannot directly applythe formalism of Vitale et al. (2017). The framework we de-velop here is able to correctly determine the relative contri-butions of multiple models, even when those models overlapin parameter space significantly, as expected in practice.There have also been significant advances in the pastfew years in the understanding of PN spin–orbit resonances.These resonances occur when BH spins become aligned oranti-aligned with one another and precess in a common planearound the total angular momentum (Schnittman 2004).This causes binaries to be attracted to different points inparameter space identified by ∆Φ, the angle between theprojections of the two BH spins onto the orbital plane. Kes-den et al. (2010) have shown that these resonances are effec-tive at capturing binaries with mass ratios 0 . < q < χ i > .
5. For equal-mass binaries, spin morphologiesremain locked with binaries trapped in or out of resonance(Gerosa et al. 2017); however, it is unlikely for astrophysicalformation scenarios to produce exactly equal mass binaries,although Marchant et al. (2016) predict nearly equal massesfor the chemically homogeneous evolution channel.Gerosa et al. (2013) show how the family of resonancesthat BBHs are attracted to can act as a diagnostic of theformation scenario for those binaries. Trifir`o et al. (2016)demonstrate that GW measurements of spin misalignmentscan be used to distinguish between the two resonant fam-ilies of ∆Φ = 0 ◦ and ∆Φ = ± ◦ . They use a full PEstudy to show that they can distinguish two families of PNresonances. However, they only consider a small corner ofparameter space which contains binaries which will becomelocked in these PN resonances. fect average PE performance for binary neutron stars (Berry et al.2015); however, these noise features could be more significant inanalysing the shorter duration BBH signals. Our study extends on those discussed in several ways:(i) Rather than focusing on specific systems preferredin previous studies, we use injections from an astrophysi-cally motivated population. Our injections have total masses M = m + m in the range 10–40 M (cid:12) and an astrophysicaldistribution of SNRs.(ii) The misalignment angles of our BHs are given by sim-ple but astrophysically motivated models introduced in Sec-tion 2.(iii) For performing PE on individual GW events, we usethe inspiral-merger-ringdown gravitational waveform IMR-PhenomPv2 model, rather than the inspiral-only waveformsused in some of the earlier studies.(iv) Most importantly, we perform a hierarchicalBayesian analysis on the posterior probability density func-tions of a mock catalog of detected events in order to makeinferences about the underlying population.
PE on individual GW events yields samples from the pos-terior distributions for parameters under astrophysical priorconstraints. We now wish to combine these individual mea-surements of BH spin misalignment angles in order to learnabout the underlying population, which may act as a diag-nostic for binary formation channels and binary evolutionscenarios. Importantly, we are able to do this without re-analysing the data for the individual events.Given a set of reasonable population synthesis modelpredictions for BBH spin misalignment angles, we would liketo learn what mixture of those subpopulations best explainsthe observed data. Here we assume that the subpopulationdistributions representing different formation channels areknown perfectly, and use the same subpopulations that wedrew our injections from to set these distributions. Thus,each subpopulation model Λ (cid:96) ( (cid:96) ∈ . . .
4) corresponds to aknown distribution of source parameters p ( Θ | Λ (cid:96) ). In prac-tice, the uncertainty in the subpopulation models will beone of the challenges in carrying out accurate hierarchicalinference. The overall mixture model is described by hyperparam-eters λ (cid:96) , corresponding to the fraction of each of the foursubpopupulations, such that p ( Θ | λ ) = (cid:88) (cid:96) =1 λ (cid:96) p ( Θ | Λ (cid:96) ) . (9)We assume that each event comes from one of these subpop-ulations: (cid:88) (cid:96) =1 λ (cid:96) = 1 , (10)i.e., λ is a unit simplex. The clustering approach of Mandel et al. (2015, 2017), which es-chews assumptions about the subpopulation distributions, couldprovide an alternative pathway for robust but less informativeinference on the data alone. MNRAS000
4) corresponds to aknown distribution of source parameters p ( Θ | Λ (cid:96) ). In prac-tice, the uncertainty in the subpopulation models will beone of the challenges in carrying out accurate hierarchicalinference. The overall mixture model is described by hyperparam-eters λ (cid:96) , corresponding to the fraction of each of the foursubpopupulations, such that p ( Θ | λ ) = (cid:88) (cid:96) =1 λ (cid:96) p ( Θ | Λ (cid:96) ) . (9)We assume that each event comes from one of these subpop-ulations: (cid:88) (cid:96) =1 λ (cid:96) = 1 , (10)i.e., λ is a unit simplex. The clustering approach of Mandel et al. (2015, 2017), which es-chews assumptions about the subpopulation distributions, couldprovide an alternative pathway for robust but less informativeinference on the data alone. MNRAS000 , 1–12 (2017)
BH spin misalignments For any individual event α ( α = 1 , . . . , N ) we have theposterior on Θ given by p ( Θ | d α ) = p ( d α | Θ ) p ( Θ ) p ( d α ) , (11)where p ( Θ ) is the prior used by LALInference , p ( d α ) is theevidence (which is only a normalising factor in our analysis),and we represent p ( Θ | d α ) by a set of discrete samples { Θ ki } where k = 1 , . . . , ν α .We can write the likelihood for obtaining all of theevents as the product over the individual likelihoods (Man-del 2010; Hogg et al. 2010), p (cid:16) { d α } Nα =1 (cid:12)(cid:12)(cid:12) λ (cid:17) = N (cid:89) α =1 p ( d α | λ ) (12)= N (cid:89) α =1 (cid:90) d Θ α p ( d α | Θ α ) p ( Θ α | λ ) (13)= N (cid:89) α =1 p ( d α ) (cid:90) d Θ α p ( Θ α | d α ) p ( Θ α ) p ( Θ α | λ ) , (14)where we have marginalised over the physical parameters ofthe individual events, and used Bayes’ theorem to obtain thefinal line. Since we have samples drawn from the posterior p ( Θ α | d α ), we can approximate posterior-weighted integrals(posterior averages) as a sum over samples (MacKay 2003,chapter 29), (cid:90) d Θ α p ( Θ α | d α ) f ( Θ α ) = 1 ν α ν α (cid:88) k =1 f ( Θ kα ) , (15)where f ( Θ ) is some general function. Thus, we can rewriteEquation (14) as p (cid:16) { d α } Nα =1 (cid:12)(cid:12)(cid:12) λ (cid:17) = N (cid:89) α =1 p ( d α ) ν α ν α (cid:88) k =1 p ( Θ kα | λ ) p ( Θ kα ) . (16)In effect, for each event we reweigh the evidence calculatedusing our general PE prior to what it would have been usinga prior for the model of interest, and then combine theseprobabilities together to form a likelihood.The posterior for λ is then p (cid:16) λ (cid:12)(cid:12)(cid:12) { d α } Nα =1 (cid:17) ∝ p (cid:16) { d α } Nα =1 (cid:12)(cid:12)(cid:12) λ (cid:17) p ( λ ) , (17)for a choice of prior p ( λ ). We assume a flat Dirichlet prior asshown in Figure 4. We sample from this posterior on λ us-ing emcee (Foreman-Mackey et al. 2013), an affine-invariantensemble sampler (Goodman & Weare 2010). To gain a qualitative understanding of hierarchical mod-elling on the spin–orbit misalignment angles, we first con-sider inference under the assumption of perfect measure-ment accuracy for individual observations, and then intro-duce realistic measurement uncertainties. We then analysethe scaling of the inference accuracy with the number ofobservations. Available from http://dan.iel.fm/emcee/.
Figure 4.
Marginalised 1D and 2D probability density func-tions for the Dirichlet prior used for the analysis of the λ pa-rameters, which describe the fractional contribution of each ofthe four subpopulations introduced in Section 2. The constraint λ + λ + λ + λ ≡ Here, we assume that aLIGO–AdV GW observations couldperfectly measure the spin–orbit misalignment angles ofmerging BBHs. In this case, the posterior is simply a deltafunction centered at the true value. Since our underlying as-trophysical models have significant overlap in the { cos( θ ),cos( θ ) } plane, as shown in Figure 1, there is still ambiguityabout which model a given event comes from.We sample Equation (17), where our data consist of 80events with perfectly measured spin–orbit misalignments (asseen in Figure 2). This number of detections could be avail-able by the end of the third observing run under optimisticassumptions about detector sensitivity improvements (Ab-bott et al. 2016b). The results of this analysis are shown inFigure 5.We find that after 80 BBH observations with perfectmeasurement accuracy, we would be able to confidently es-tablish the presence of all four subpopulations. From thisanalysis, we can already understand some of the featuresof the posterior on the hyperparameters. For example, wesee that there is a strong degeneracy between λ and λ ,since both of these models predict a large (nearly) aligned( θ = θ = 0) population. There is a similar degeneracybetween λ and λ . We can also see that the fraction of ex-actly aligned systems ( λ ) and the fraction of systems withisotropically distributed spin–orbit misalignments ( λ ) arenot strongly correlated. Both fractions are measured withto be between ∼ .
15 and ∼ .
45 at the 90% credible levelwith 80 BBH observations, corresponding to a fractional un-certainty of ∼ MNRAS , 1–12 (2017)
S. Stevenson et al.
Figure 5.
Marginalised 1D and 2D probability density functionsfor the λ parameters describing the fractional contribution of eachof the four subpopulations introduced in Section 2. The thin blacklines indicate the true injection fraction from each model, whichis 0 .
25 for all models. The data used were the 80 mock GW eventsshown in Figure 2, assumed to have perfect measurements of thespin–orbit misalignment angles cos θ and cos θ . Colours are thesame as Figure 4 We know that in practice GW detectors will not perfectlymeasure the spin–orbit misalignments of merging BBHs (seeFigure 3 for a typical marginalised posterior). We now usethe full set of 80
LALInference posteriors, each containing ∼ λ and λ . We see that the posterioris not perfectly centred on the true λ values, though thetrue values do have posterior support. While the hierarchicalmodelling unambiguously points to the presence of multiplesubpopulations, with no single subpopulation able to explainthe full set of observations, the data no longer require all foursubpopulations to be present.We have checked that the structure of this posterior istypical given the limited number of observations and thelarge measurement uncertainties. In the next section weshow that our posteriors converge to the true values in thelimit of a large number of detections. The
LALInference
PE pipeline used to compute the poste-rior distributions for our 80 injections in Section 3 is com-putationally expensive. However, we would like to generatea larger catalogue of mock observations. First, this allows usto check that our analysis is self consistent by running manytests, such as confirming that the true result lies within the
Figure 6.
Marginalised 1D and 2D probability density functionsfor the λ parameters describing the fractional contribution of eachof the four subpopulations introduced in Section 2. The thin blacklines indicate the true injection fraction from each model, whichis 0 .
25 for all models. The data used were the full
LALInference posteriors of the 80 mock GW events shown in Figure 2. Coloursare the same as Figure 4 P % credible interval in P % of trials. Second, it allows usto predict how the accuracy of the inferred fractions of thesubpopulations evolves as a function of the number of GWobservations.We develop approximations to these posteriors, similarto Mandel et al. (2017), based on the 80 posterior distribu-tions generated in Section 3. The best measured spin pa-rameter is a combination of the two component spins calledthe effective inspiral spin χ eff ∈ [ − ,
1] (Ajith et al. 2011;Abbott et al. 2016e; Vitale et al. 2017): χ eff = χ cos θ + qχ cos θ (1 + q ) . (18)Having information about a single spin parameter makes itchallenging to extract information about the spin distribu-tion, but not impossible; for example, GW151226’s positive χ eff means that at least one spin must have non-zero mag-nitude and θ i < ◦ (Abbott et al. 2016f).To compute the approximate posteriors, we representeach observation with true parameter values χ trueeff andcos θ by data which are maximum-likelihood estimatesin a random noise realization: χ dataeff ∼ N (cid:0) χ trueeff , σ χ eff ( χ trueeff ) (cid:1) , (19)cos θ ∼ N (cid:0) cos θ , σ θ (cos θ ) (cid:1) , (20)where N ( µ, σ ) indicates a normal distribution. Posteriorsamples are then drawn using the same likelihood functions,centred on the maximum-likelihood data value, χ sampleeff ∼ N (cid:16) χ dataeff , σ χ eff ( χ sampleeff ) (cid:17) , (21)cos θ ∼ N (cid:16) cos θ , σ θ (cos θ ) (cid:17) . (22) MNRAS , 1–12 (2017)
BH spin misalignments Figure 7.
Marginalised posterior samples for the same event asshown in Figure 3, with the same notation. This posterior distri-bution was approximated using the model described in section 5.3.
Here σ cos θ = 12 ρ ( A cos θ + B ) , (23) σ χ eff = 12 ρ C, (24)with A = − . B = 0 . C = 0 . ρ (Cutler & Flanagan 1994).In all cases, we only consider cos θ , cos θ and χ eff in thepermitted range of [ − , θ and θ . We havechosen a particular strategy that allows us to make faithfulmodels; starting with models for posteriors on χ eff and θ before converting them into posteriors on θ and θ is desir-able because χ eff and θ are only weakly correlated, makingit possible to write down independent expressions for the in-dividual likelihoods, while θ and θ are strongly correlated,so would require a more complex joint likelihood model.The posteriors already implicitly account for marginaliza-tion over correlated parameters such as spin magnitudes andthe mass ratio, which will not be known in practice. Thesemock posteriors will not correctly reproduce actual posteri-ors for specific events, since they do not include the noiserealization, except probabilistically; they only reproduce theaverage properties of the posterior distributions.We draw ν = 5000 posterior samples of χ eff and cos θ independently, and calculate the values of cos θ using Equa-tion (18), fixing the mass ratio q and spin magnitudes χ i totheir true values. This builds the correct degeneracies be-tween cos θ and cos θ into the mock posteriors. We showan example of a posterior distribution generated this way inFigure 7.Using this method, we generate spin–orbit misalignmentposteriors for 400 BBHs drawn in equal fractions ( λ i = 0 . λ parameters after 0 (prior), 10, 20, 40, 80,160 and 400 observations, similar to Mandel et al. (2017). InFigure 8 we show the 68% and 95% credible intervals for thefraction λ of observed BBHs coming from an isotropic dis-tribution (subpopulation 2) corresponding to dynamicallyformed binaries. Given our models and incorporating real- Figure 8.
Marginalised posterior on λ , the fraction of BBHsfrom the subpopulation with isotropic spin distribution (repre-senting dynamical formation), as a function of the number ofGW observations. The posterior converges to the injected valueof λ = 0 .
25 (dashed black horizontal line) after ∼
100 obser-vations. The coloured bands show the 68% (darkest) and 95%(lightest) credible intervals.
Figure 9.
Marginalised posterior on λ + λ , the combined frac-tion of BBHs formed through subpopulations 1 and 3, as a func-tion of the number of GW observations. These subpopulationscorrespond to spins preferentially aligned with the orbital angu-lar momentum. The posterior converges to the injected value of λ + λ = 0 . ∼
20 obser-vations. The coloured bands show the 68% (darkest) and 95%(lightest) credible intervals. istic measurement uncertainties, we find that this fractioncan be measured with a ∼
40% fractional uncertainty after100 observations. Since subpopulation 1 and 3 are somewhatdegenerate in our model, we find that the combined fraction λ + λ is a well measured parameter (as shown in Figure 9),whilst the individual components are measured less well. For N (cid:38)
100 observations the uncertainties in the λ (cid:96) scale asthe inverse square root of the number of observations; e.g., σ λ + λ ≈ . N − / .Although (cid:38)
100 observations are required to accuratelymeasure the contribution of each of the four subpopulations,it is possible to test for more extreme models with fewerobservations. For example, ∼
20 observations are sufficientto demonstrate the presence of an isotropic subpopulationat the 95% credible level.Even fewer observations are needed to confidently rule
MNRAS , 1–12 (2017) S. Stevenson et al. out the hypothesis that all observations come from the ex-actly aligned or isotropic subpopulations. We draw obser-vations from the isotropic subpopulation and calculate theratio of the evidence (Bayes factor) Z aligned for the modelunder which all BBH spins are exactly aligned ( λ = 1) tothe evidence Z isotropic for the model under which all BBHspins are isotropically distributed ( λ = 1): Z aligned Z isotropic = p (cid:16) { d α } Nα =1 (cid:12)(cid:12)(cid:12) λ = 1 (cid:17) p (cid:16) { d α } Nα =1 (cid:12)(cid:12)(cid:12) λ = 1 (cid:17) . (25)Inference using small numbers of observations is sen-sitive to the exact choice of these observations (both trueparameters and measurement errors), which we randomlydraw from the relevant distributions. We therefore repeatthis test 100 times to account for the random nature of themock catalogue. In all cases, we find that with only 5 obser-vations of BBHs with component spin magnitudes χ = 0 . λ = 1 can be ruled out at morethan 5 σ confidence. Similarly, when drawing from the ex-actly aligned model, we find that the hypothesis that allevents come from an isotropic population λ = 1 can beruled out at more than 5 σ confidence in all tests with 5observations. With the first direct observations of GWs from mergingBBHs, the era of GW astronomy has begun. GW observa-tions provide a new and unique insight into the properties ofBBHs and their progenitors. For individual systems, we caninfer the masses and spins of the component BHs; combiningthese measurements we can learn about the population, andplace constraints on the formation mechanisms for these sys-tems, whether as the end point of isolated binary evolutionor as the results of dynamical interactions.In this work, we investigated how measurements of BBHspin–orbit misalignments could inform our understanding ofthe BBH population. We chose the properties of our sourcesto match those we hope to observe with aLIGO and AdV(at design sensitivity), using four different astrophysicallymotivated subpopulations for the distribution of spin–orbitmisalignment angles, each reflecting a different formationscenario. We performed a Bayesian analysis of GW signals(using full inspiral–merger–ringdown waveforms) for a popu-lation of BBHs. We assumed a mixture model for the overallpopulation of BBH spin–orbit misalignments and combinedthe full PE results from our GW analysis in a hierarchicalframework to infer the fraction of the population comingfrom each subpopulation. A similar analysis could be per-formed following the detection of real signals.Adopting a population with spins of χ = 0 .
7, we demon-strate that the fraction of BBHs with spins preferentiallyaligned with the orbital angular momentum ( λ + λ ) iswell measured and can be measured with an uncertainty of ∼
10% with 100 observations, scaling as the inverse squareroot of the number of observations. We also show that af-ter 100 observations, we can measure the fraction λ of thesubpopulation with isotropic spins (assumed to correspondto dynamical formation) with a fractional uncertainty of ∼ > σ ) ifthe true population has isotropically distributed spins, andvice versa. This number of observations may be reached bythe end of the second aLIGO observing run.One limitation of the current approach is the assump-tion that the subpopulation distributions are known per-fectly. This will not be the case in practice, but the simplifiedmodels considered here are still relevant as parametrizableproxies for astrophysical scenarios. Hierarchical modellingwith strong population assumptions could lead to system-atic biases in the interpretation of the observations if thoseassumptions are not representative of the true populations;this can be mitigated by coupling such analysis with weaklymodelled approaches, such as observation-based clustering(Mandel et al. 2015, 2017).In this work we have not taken into account observa-tional selection effects, particularly the differences in thedetectabilty of different subpopulations because their dif-ferent spin parameter distributions impact the SNR (Reis-swig et al. 2009). These must be incorporated in the anal-ysis to correctly infer the intrinsic subpopulation fractions(Farr et al. 2015; Mandel et al. 2016), but the impact is notexpected to be significant in this case. Care must also betaken to avoid biases when performing an hierarchical anal-ysis with real observations, since the observations will not bedrawn from the same distribution as the priors used for theanalysis of individual events. Our framework accounts forthe differences in the priors on the parameters of interest(spin–orbit misalignment angles) between the original PEand model predictions, but not for any discrepancy in thepriors of the parameters we marginalize over (e.g., masses);this is likely a second-order effect.Neither theoretical models nor observations can cur-rently place tight constraints on the spin magnitudes of BHs.We therefore chose to give all BHs a spin magnitude of 0 . − MNRAS , 1–12 (2017)
BH spin misalignments driguez et al. 2016) and spin magnitudes (cf. Gerosa & Berti2017; Fishbach et al. 2017), into a single analysis. Comple-mentary electromagnetic observations of high-mass X-raybinaries, Galactic radio pulsars, short gamma ray bursts,supernovae and luminous red novae will contribute to a con-cordance model of massive binary formation and evolution. ACKNOWLEDGEMENTS
This work was supported in part by the Science and Tech-nology Facilities Council. IM acknowledges support from theLeverhulme Trust. We are grateful for computational re-sources provided by Cardiff University, and funded by anSTFC grant supporting UK Involvement in the Operationof Advanced LIGO. SS and IM acknowledge partial supportby the National Science Foundation under Grant No. NSFPHY11-25915. We are grateful to colleagues from the Insti-tute of Gravitational-wave Astronomy at the University ofBirmingham, as well as Ben Farr, M. Coleman Miller, ChrisPankow and Salvatore Vitale for fruitful discussions. Wethank Davide Gerosa and the anonymous referee for com-ments on the paper. Corner plots were made using trian-gle.py available from https://github.com/dfm/triangle.py.This is LIGO Document P1700007.
REFERENCES
Aasi J., et al., 2015, Class. Quantum Grav., 32, 074001Aasi J., et al., 2016, Living Rev. Rel., 19, 1Abbott B. P., et al., 2016a, Phys. Rev. X, 6, 041014Abbott B. P., et al., 2016b, Phys. Rev. X, 6, 041015Abbott B. P., et al., 2016c, Phys. Rev. D, 93, 122003Abbott B. P., et al., 2016d, Phys. Rev. Lett., 116, 061102Abbott B. P., et al., 2016e, Phys. Rev. Lett., 116, 241102Abbott B. P., et al., 2016f, Phys. Rev. Lett., 116, 241103Abbott B. P., et al., 2016g, ApJ, 818, L22Abbott B. P., et al., 2016h, ApJ, 833, L1Abbott B. P., et al., 2017, Phys. Rev. Lett., 118, 221101Abt H. A., 1983, ARA&A, 21, 343Acernese F., et al., 2015, Class. Quantum Grav., 32, 024001Ajith P., et al., 2011, Phys. Rev. Lett., 106, 241101Albrecht S., Reffert S., Snellen I., Quirrenbach A., Mitchell D. S.,2007, A&AAlbrecht S., Reffert S., Snellen I. A. G., Winn J. N., 2009, Nature,461, 373Albrecht S., et al., 2014, ApJ, 785, 83Apostolatos T. A., Cutler C., Sussman G. J., Thorne K. S., 1994,Phys. Rev. D, 49, 6274Aso Y., Michimura Y., Somiya K., Ando M., Miyakawa O.,Sekiguchi T., Tatsumi D., Yamamoto H., 2013, Phys. Rev.D, 88, 043007Baird E., Fairhurst S., Hannam M., Murphy P., 2013, Phys. Rev.D, 87, 024035Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., VinkJ. S., Hurley J. R., 2010, ApJ, 714, 1217Belczynski K., Bulik T., Mandel I., Sathyaprakash B., ZdziarskiA. A., et al., 2013, ApJ, 764, 96Berry C. P. L., et al., 2015, ApJ, 804, 114Blaauw A., 1961, Bulletin of the Astronomical Institutes of theNetherlands, 15, 265Blanchet L., 2014, Living Rev. Relat., 17, 2Boss A. P., 1988, Comments on Astrophysics, 12, 169Brandt N., Podsiadlowski P., 1995, MNRAS, 274, 461 Bulik T., Belczynski K., Prestwich A., 2011, ApJ, 730, 140Crowther P. A., Barnard R., Carpano S., Clark J. S., DhillonV. S., Pollock A. M. T., 2010, MNRAS, 403, L41Cutler C., Flanagan E. E., 1994, Phys. Rev. D, 49, 2658de Mink S. E., Belczynski K., 2015, ApJ, 814, 58de Mink S. E., Langer N., Izzard R. G., Sana H., de Koter A.,2013, ApJ, 764, 166Dominik M., Belczynski K., Fryer C., Holz D., Berti E., et al.,2012, ApJ, 759, 52Farr W. M., Kremer K., Lyutikov M., Kalogera V., 2011, Astro-phys. J., 742, 81Farr W. M., Mandel I., Aldridge C., Stroud K., 2014a, preprint,( arXiv:1412.4849 )Farr B., Ochsner E., Farr W. M., O’Shaughnessy R., 2014b, Phys.Rev. D, 90, 024018Farr W. M., Gair J. R., Mandel I., Cutler C., 2015, Phys. Rev.D, 91, 023005Farr B., et al., 2016, ApJ, 825, 116Farr W. M., Stevenson S., Miller M. C., Mandel I., Farr B., Vec-chio A., 2017, preprint, ( arXiv:1706.01385 )Finn L. S., 1992, Phys. Rev., D46, 5236Finn L. S., Chernoff D. F., 1993, Phys. Rev. D, 47, 2198Fishbach M., Holz D., Farr B., 2017, preprint,( arXiv:1703.06869 )Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013,PASP, 125, 306Foreman-Mackey D., Hogg D. W., Morton T. D., 2014, ApJ, 795,64Fragos T., McClintock J. E., 2015, ApJ, 800, 17Fragos T., Tremmel M., Rantsiou E., Belczynski K., 2010, ApJ,719, L79Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M.,Kalogera V., et al., 2012, ApJ, 749, 91Gerosa D., Berti E., 2017, preprint, ( arXiv:1703.06223 )Gerosa D., Kesden M., Berti E., O’Shaughnessy R., Sperhake U.,2013, Phys. Rev. D, 87, 104028Gerosa D., O’Shaughnessy R., Kesden M., Berti E., Sperhake U.,2014, Phys. Rev. D, 89, 124025Gerosa D., Kesden M., Sperhake U., Berti E., O’Shaughnessy R.,2015, Phys. Rev. D, 92, 064016Gerosa D., Sperhake U., Voˇsmera J., 2017, Classical and QuantumGravity, 34, 064004Goodman J., Weare J., 2010, Communications in Applied Math-ematics and Computational Science, 5, 65Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS,360, 974Hogg D. W., Myers A. D., Bovy J., 2010, ApJ, 725, 2166Hut P., 1981, A&A, 99, 126Ivanova N., et al., 2013, A&ARv, 21, 59Kalogera V., 1996, ApJ, 471, 352Kalogera V., 2000, ApJ, 541, 319Kelley L. Z., Ramirez-Ruiz E., Zemp M., Diemand J., Mandel I.,2010, ApJ, 725, L91Kesden M., Sperhake U., Berti E., 2010, Phys. Rev. D, 81, 084054Kesden M., Gerosa D., O’Shaughnessy R., Berti E., Sperhake U.,2015, Phys. Rev. Lett., 114, 081103King A. R., Lubow S. H., Ogilvie G. I., Pringle J. E., 2005, MN-RAS, 363, 49Kroupa P., 2001, MNRAS, 322, 231Kruckow M. U., Tauris T. M., Langer N., Sz´ecsi D., Marchant P.,Podsiadlowski P., 2016, A&A, 596, A58Kushnir D., Zaldarriaga M., Kollmeier J. A., Waldman R., 2016,MNRAS, 462, 844Lieu M., Farr W. M., Betancourt M., Smith G. P., Sereno M.,McCarthy I. G., 2017, MNRAS, 468, 4872Lipunov V. M., Postnov K. A., Prokhorov M. E., BogomazovA. I., 2009, Astronomy Reports, 53, 915MNRAS , 1–12 (2017) S. Stevenson et al.
Littenberg T. B., Farr B., Coughlin S., Kalogera V., Holz D. E.,2015, ApJ, 807, L24Loredo T. J., 2004, in Fischer R., Preuss R., Toussaint U. V.,eds, American Institute of Physics Conference Series Vol. 735,American Institute of Physics Conference Series. pp 195–206( arXiv:astro-ph/0409387 ), doi:10.1063/1.1835214MacKay D. J. C., 2003, Information Theory, Inference and Learn-ing Algorithms. Cambridge University Press, CambridgeMandel I., 2010, Phys. Rev. D, 81, 084029Mandel I., 2016, MNRAS, 456, 578Mandel I., O’Shaughnessy R., 2010, Class. Quantum Grav., 27,114007Mandel I., Haster C.-J., Dominik M., Belczynski K., 2015, MN-RAS, 450, L85Mandel I., Farr W. M., Gair J. R., 2016, Extracting dis-tribution parameters from multiple uncertain observationswith selection biases. In prep., https://dcc.ligo.org/LIGO-P1600187/publicMandel I., Farr W. M., Colonna A., Stevenson S., Tiˇno P., VeitchJ., 2017, MNRAS, 465, 3254Marchant P., Langer N., Podsiadlowski P., Tauris T. M., MoriyaT. J., 2016, A&A, 588, A50Martin R. G., Reis R. C., Pringle J. E., 2008, MNRAS, 391, L15McClintock J. E., et al., 2011, Class. Quantum Grav., 28, 114009Miller M. C., Miller J. M., 2015, Phys. Rep., 548, 1Miller B., O’Shaughnessy R., Littenberg T. B., Farr B., 2015,Phys. Rev. D, 92, 044056Mirabel I. F., 2016, preprint, ( arXiv:1609.08411 )Moore C. J., Cole R. H., Berry C. P. L., 2015, Class. QuantumGrav., 32, 015014Naoz S., Farr W. M., Rasio F. A., 2012, ApJ, 754, L36O’Shaughnessy R. W., Kim C., Kalogera V., Belczynski K., 2008,ApJ, 672, 479Orosz J. A., et al., 2001, ApJ, 555, 489Peters P. C., 1964, Phys. Rev., 136, 1224Peters P. C., Mathews J., 1963, Phys. Rev., 131, 435Poisson E., Will C. M., 1995, Phys. Rev. D, 52, 848Postnov K. A., Yungelson L. R., 2014, Living Rev. Relat., 17, 3Ram´ırez-Agudelo O. H., et al., 2015, A&A, 580, A92Reisswig C., Husa S., Rezzolla L., Dorband E. N., Pollney D.,Seiler J., 2009, Phys. Rev. D, 80, 124026Repetto S., Davies M. B., Sigurdsson S., 2012, MNRAS, 425, 2799Repetto S., Igoshev A. P., Nelemans G., 2017, MNRAS,Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S.,Haster C.-J., Rasio F. A., 2015, Phys. Rev. Lett., 115, 051101Rodriguez C. L., Zevin M., Pankow C., Kalogera V., Rasio F. A.,2016, ApJ, 832, L2Sadowski A., Belczynski K., Bulik T., Ivanova N., Rasio F. A.,O’Shaughnessy R., 2008, ApJ, 676, 1162Salpeter E. E., 1955, ApJ, 121, 161Schmidt P., Ohme F., Hannam M., 2015, Phys. Rev. D, 91, 024043Schnittman J. D., 2004, Phys. Rev. D, 70, 124020Schutz B. F., 2011, Class. Quantum Grav., 28, 125023Sigurdsson S., Hernquist L., 1993, Nature, 364, 423Stevenson S., Ohme F., Fairhurst S., 2015, ApJ, 810, 58Trifir`o D., O’Shaughnessy R., Gerosa D., Berti E., Kesden M.,Littenberg T., Sperhake U., 2016, Phys. Rev. D, 93, 044071Vallisneri M., 2011, Phys. Rev. Lett., 107, 191104Veitch J., et al., 2015a, Phys. Rev. D, 91, 042003Veitch J., P¨urrer M., Mandel I., 2015b, Phys. Rev. Lett., 115,141101Vitale S., Lynch R., Sturani R., Graff P., 2017, Class. QuantumGrav., 34, 03LT01Voss R., Tauris T. M., 2003, MNRAS, 342, 1169Woosley S. E., Heger A., Weaver T. A., 2002, Rev. Mod. Phys.,74, 1015Zahn J.-P., 1977, A&A, 57, 383 This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000