<figure><img src="image_1.png" /></figure>

HAL
open science

# Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo

B.P. Abbott, R. Abbott, T.D. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R.X. Adhikari, V.B. Adya, C. Affeldt, et al.

► To cite this version:

B.P. Abbott, R. Abbott, T.D. Abbott, S. Abraham, F. Acernese, et al., Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo. Astrophys.J.Lett., 2019, 882 (2), pp.L24. 10.3847/2041-8213/ab3800. hal-02051700

HAL Id: hal-02051700

https://hal.science/hal-02051700v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

CrossMark

Binary Black Hole Population Properties Inferred from the First and Second Observing
Runs of Advanced LIGO and Advanced Virgo

B. P. Abbott¹, R. Abbott¹, T. D. Abbott², S. Abraham³, F. Acernese⁴,⁵, K. Ackley⁶, C. Adams⁷, R. X. Adhikari¹, V. B. Adya⁸,⁹,
C. Affeldt⁸,⁹, M. Agathos¹⁰, K. Agatsuma¹¹, N. Aggarwal¹², O. D. Aguiar¹³, L. Aiello¹⁴,¹⁵, A. Ain³, P. Ajith¹⁶, G. Allen¹⁷,
A. Allocca¹⁸,¹⁹, M. A. Aloy²⁰, P. A. Altin²¹, A. Amato²², A. Ananyeva¹, S. B. Anderson¹, W. G. Anderson²³, S. V. Angelova²⁴,
S. Antier²⁵, S. Appert¹, K. Arai¹, M. C. Araya¹, J. S. Areeda²⁶, M. Arène²⁷, N. Arnaud²⁵,²⁸, K. G. Arun²⁹, S. Ascenzi³⁰,³¹,
G. Ashton⁶, S. M. Aston⁷, P. Astone³², F. Aubin³³, P. Aufmuth⁹, K. AultONeal³⁴, C. Austin², V. Avendano³⁵, A. Avila-Alvarez²⁶,
S. Babak²⁷,³⁶, P. Bacon²⁷, F. Badaracco¹⁴,¹⁵, M. K. M. Bader³⁷, S. Bae³⁸, P. T. Baker³⁹, F. Baldaccini⁴⁰,⁴¹, G. Ballardin²⁸,
S. W. Ballmer⁴², S. Banagiri⁴³, J. C. Barayoga¹, S. E. Barclay⁴⁴, B. C. Barish¹, D. Barker⁴⁵, K. Barkett⁴⁶, S. Barnum¹², F. Barone⁴,⁵,
B. Barr⁴⁴, L. Barsotti¹², M. Barsuglia²⁷, D. Barta⁴⁷, J. Bartlett⁴⁵, I. Bartos⁴⁸, R. Bassiri⁴⁹, A. Basti¹⁸,¹⁹, M. Bawaj⁴¹,⁵⁰,
J. C. Bayley⁴⁴, M. Bazzan⁵¹,⁵², B. Bécsy⁵³, M. Bejger²⁷,⁵⁴, I. Belahcene²⁵, A. S. Bell⁴⁴, D. Beniwal⁵⁵, B. K. Berger⁴⁹,
G. Bergmann⁸,⁹, S. Bernuzzi⁵⁶,⁵⁷, J. J. Bero⁵⁸, C. P. L. Berry⁵⁹, D. Bersanetti⁶⁰, A. Bertolini³⁷, J. Betzwieser⁷, R. Bhandare⁶¹,
J. Bidler²⁶, I. A. Bilenko⁶², S. A. Bilgili³⁹, G. Billingsley¹, J. Birch⁷, I. A. Birney²⁴, O. Birnholtz⁵⁸, S. Biscans¹,¹², S. Biscoveanu⁶,
A. Bisht⁹, M. Bitossi¹⁹,²⁸, M. A. Bizouard²⁵, J. K. Blackburn¹, C. D. Blair⁷, D. G. Blair⁶³, R. M. Blair⁴⁵, S. Bloemen⁶⁴, N. Bode⁸,⁹,
M. Boer⁶⁵, Y. Boetzel⁶⁶, G. Bogaert⁶⁵, F. Bondu⁶⁷, E. Bonilla⁴⁹, R. Bonnand³³, P. Booker⁸,⁹, B. A. Boom³⁷, C. D. Booth⁶⁸,
R. Bork¹, V. Boschi²⁸, S. Bose³,⁶⁹, K. Bossie⁷, V. Bossilkov⁶³, J. Bosveld⁶³, Y. Bouffanais²⁷, A. Bozzi²⁸, C. Bradaschia¹⁹,
P. R. Brady²³, A. Bramley⁷, M. Branchesi¹⁴,¹⁵, J. E. Brau⁷⁰, T. Briant⁷¹, J. H. Briggs⁴⁴, F. Brighenti⁷²,⁷³, A. Brillet⁶⁵,
M. Brinkmann⁸,⁹, V. Brisson²⁵,¹⁸¹, P. Brockill²³, A. F. Brooks¹, D. D. Brown⁵⁵, S. Brunett¹, A. Buikema¹², T. Bulik⁷⁴,
H. J. Bulten³⁷,⁷⁵, A. Buonanno³⁶,⁷⁶, R. Buscicchio¹¹, D. Buskulic³³, C. Buy²⁷, R. L. Byer⁴⁹, M. Cabero⁸,⁹, L. Cadonati⁷⁷,
G. Cagnoli²²,⁷⁸, C. Cahillane¹, J. Calderón Bustillo⁶, T. A. Callister¹, E. Calloni⁵,⁷⁹, J. B. Camp⁸⁰, W. A. Campbell⁶,
M. Canepa⁶⁰,⁸¹, K. C. Cannon⁸², H. Cao⁵⁵, J. Cao⁸³, E. Capocasa²⁷, F. Carbognani²⁸, S. Caride⁸⁴, M. F. Carney⁵⁹, G. Carullo¹⁸,
J. Casanueva Diaz¹⁹, C. Casentini³⁰,³¹, S. Caudill³⁷, M. Cavaglià⁸⁵, F. Cavalier²⁵, R. Cavalieri²⁸, G. Cella¹⁹, P. Cerdá-Durán²⁰,
G. Cerretani¹⁸,¹⁹, E. Cesarini³¹,⁸⁶, O. Chaibi⁶⁵, K. Chakravarti³, S. J. Chamberlin⁸⁷, M. Chan⁴⁴, S. Chao⁸⁸, P. Charlton⁸⁹,
E. A. Chase⁵⁹, E. Chassande-Mottin²⁷, D. Chatterjee²³, M. Chaturvedi⁶¹, K. Chatziioannou⁹⁰, B. D. Cheeseboro³⁹, H. Y. Chen⁹¹,
X. Chen⁶³, Y. Chen⁴⁶, H.-P. Cheng⁴⁸, C. K. Cheong⁹², H. Y. Chia⁴⁸, A. Chincarini⁶⁰, A. Chiummo²⁸, G. Cho⁹³, H. S. Cho⁹⁴,
M. Cho⁷⁶, N. Christensen⁶⁵,⁹⁵, Q. Chu⁶³, S. Chua⁷¹, K. W. Chung⁹², S. Chung⁶³, G. Ciani⁵¹,⁵², A. A. Ciobanu⁵⁵, R. Ciolfi⁹⁶,⁹⁷
F. Cipriano⁶⁵, A. Cirone⁶⁰,⁸¹, F. Clara⁴⁵, J. A. Clark⁷⁷, P. Clearwater⁹⁸, F. Cleva⁶⁵, C. Cocchieri⁸⁵, E. Coccia¹⁴,¹⁵, P.-F. Cohadon⁷¹,
D. Cohen²⁵, R. Colgan⁹⁹, M. Colleoni¹⁰⁰, C. G. Collette¹⁰¹, C. Collins¹¹, L. R. Cominsky¹⁰², M. Constancio Jr.¹³, L. Conti⁵²,
S. J. Cooper¹¹, P. Corban⁷, T. R. Corbitt², I. Cordero-Carrión¹⁰³, K. R. Corley⁹⁹, N. Cornish⁵³, A. Corsi⁸⁴, S. Cortese²⁸,
C. A. Costa¹³, R. Cotesta³⁶, M. W. Coughlin¹, S. B. Coughlin⁵⁹,⁶⁸, J.-P. Coulon⁶⁵, S. T. Countryman⁹⁹, P. Couvares¹,
P. B. Covas¹⁰⁰, E. E. Cowan⁷⁷, D. M. Coward⁶³, M. J. Cowart⁷, D. C. Coyne¹, R. Coyne¹⁰⁴, J. D. E. Creighton²³,
T. D. Creighton¹⁰⁵, J. Cripe², M. Croquette⁷¹, S. G. Crowder¹⁰⁶, T. J. Cullen², A. Cumming⁴⁴, L. Cunningham⁴⁴, E. Cuoco²⁸,
T. Dal Canton⁸⁰, G. Dálya¹⁰⁷, S. L. Danilishin⁸,⁹, S. D'Antonio³¹, K. Danzmann⁸,⁹, A. Dasgupta¹⁰⁸, C. F. Da Silva Costa⁴⁸,
L. E. H. Datrier⁴⁴, V. Dattilo²⁸, I. Dave⁶¹, M. Davier²⁵, D. Davis⁴², E. J. Daw¹⁰⁹, D. DeBra⁴⁹, M. Deenadayalan³, J. Degallaix²²,
M. De Laurentis⁵,⁷⁹, S. Deléglise⁷¹, W. Del Pozzo¹⁸,¹⁹, L. M. DeMarchi⁵⁹, N. Demos¹², T. Dent⁸,⁹,¹¹⁰, R. De Pietri⁵⁷,¹¹¹, J. Derby²⁶,
R. De Rosa⁵,⁷⁹, C. De Rossi²²,²⁸, R. DeSalvo¹¹², O. de Varona⁸,⁹, S. Dhurandhar³, M. C. Díaz¹⁰⁵, T. Dietrich³⁷, L. Di Fiore⁵,
M. Di Giovanni⁹⁷,¹¹³, T. Di Girolamo⁵,⁷⁹, A. Di Lieto¹⁸,¹⁹, B. Ding¹⁰¹, S. Di Pace³²,¹¹⁴, I. Di Palma³²,¹¹⁴, F. Di Renzo¹⁸,¹⁹,
A. Dmitriev¹¹, Z. Doctor⁹¹, F. Donovan¹², K. L. Dooley⁶⁸,⁸⁵, S. Doravari⁸,⁹, I. Dorrington⁶⁸, T. P. Downes²³, M. Drago¹⁴,¹⁵,
J. C. Driggers⁴⁵, Z. Du⁸³, J.-G. Ducoin²⁵, P. Dupej⁴⁴, S. E. Dwyer⁴⁵, P. J. Easter⁶, T. B. Edo¹⁰⁹, M. C. Edwards⁹⁵, A. Effler⁷,
P. Ehrens¹, J. Eichholz¹, S. S. Eikenberry⁴⁸, M. Eisenmann³³, R. A. Eisenstein¹², R. C. Essick⁹¹, H. Estelles¹⁰⁰, D. Estevez³³,
Z. B. Etienne³⁹, T. Etzel¹, M. Evans¹², T. M. Evans⁷, V. Fafone¹⁴,³⁰,³¹, H. Fair⁴², S. Fairhurst⁶⁸, X. Fan⁸³, S. Farinon⁶⁰, B. Farr⁷⁰,
W. M. Farr¹¹, E. J. Fauchon-Jones⁶⁸, M. Favata³⁵, M. Fays¹⁰⁹, M. Fazio¹¹⁵, C. Fee¹¹⁶, J. Feicht¹, M. M. Fejer⁴⁹, F. Feng²⁷,
A. Fernandez-Galiana¹², I. Ferrante¹⁸,¹⁹, E. C. Ferreira¹³, T. A. Ferreira¹³, F. Ferrini²⁸, F. Fidecaro¹⁸,¹⁹, I. Fiori²⁸, D. Fiorucci²⁷,
M. Fishbach⁹¹, R. P. Fisher⁴²,¹¹⁷, J. M. Fishner¹², M. Fitz-Axen⁴³, R. Flaminio³³,¹¹⁸, M. Fletcher⁴⁴, E. Flynn²⁶, H. Fong⁹⁰,
J. A. Font²⁰,¹¹⁹, P. W. F. Forsyth²¹, J.-D. Fournier⁶⁵, S. Frasca³²,¹¹⁴, F. Frasconi¹⁹, Z. Frei¹⁰⁷, A. Freise¹¹, R. Frey⁷⁰, V. Frey²⁵,
P. Fritschel¹², V. V. Frolov⁷, P. Fulda⁴⁸, M. Fyffe⁷, H. A. Gabbard⁴⁴, B. U. Gadre³, S. M. Gaebel¹¹, J. R. Gair¹²⁰, L. Gammaitoni⁴⁰,
M. R. Ganija⁵⁵, S. G. Gaonkar³, A. Garcia²⁶, C. García-Quirós¹⁰⁰, F. Garufi⁵,⁷⁹, B. Gateley⁴⁵, S. Gaudio³⁴, G. Gaur¹²¹,
V. Gayathri¹²², G. Gemme⁶⁰, E. Genin²⁸, A. Gennai¹⁹, D. George¹⁷, J. George⁶¹, L. Gergely¹²³, V. Germain³³, S. Ghonge⁷⁷,
Abhirup Ghosh¹⁶, Archisman Ghosh³⁷, S. Ghosh²³, B. Giacomazzo⁹⁷,¹¹³, J. A. Giaime²,⁷, K. D. Giardina⁷, A. Giazotto¹⁹,¹⁸²,
K. Gill³⁴, G. Giordano⁴,⁵, L. Glover¹¹², P. Godwin⁸⁷, E. Goetz⁴⁵, R. Goetz⁴⁸, B. Goncharov⁶, G. González²,
J. M. Gonzalez Castro¹⁸,¹⁹, A. Gopakumar¹²⁴, M. L. Gorodetsky⁶², S. E. Gossan¹, M. Gosselin²⁸, R. Gouaty³³, A. Grado⁵,¹²⁵,
C. Graef⁴⁴, M. Granata²², A. Grant⁴⁴, S. Gras¹², P. Grassia¹, C. Gray⁴⁵, R. Gray⁴⁴, G. Greco⁷²,⁷³, A. C. Green¹¹,⁴⁸, R. Green⁶⁸,
E. M. Gretarsson³⁴, P. Groot⁶⁴, H. Grote⁶⁸, S. Grunewald³⁶, P. Gruning²⁵, G. M. Guidi⁷²,⁷³, H. K. Gulati¹⁰⁸, Y. Guo³⁷, A. Gupta⁸⁷,
M. K. Gupta¹⁰⁸, E. K. Gustafson¹, R. Gustafson¹²⁶, L. Haegel¹⁰⁰, O. Halim¹⁴,¹⁵, B. R. Hall⁶⁹, E. D. Hall¹², E. Z. Hamilton⁶⁸,
G. Hammond<sup>44</sup>, M. Haney<sup>66</sup>, M. M. Hanke<sup>8,9</sup>, J. Hanks<sup>45</sup>, C. Hanna<sup>87</sup>, M. D. Hannam<sup>68</sup>, O. A. Hannuksela<sup>92</sup>, J. Hanson<sup>7</sup>,
T. Hardwick<sup>2</sup>, K. Haris<sup>16</sup>, J. Harms<sup>14,15</sup>, G. M. Harry<sup>127</sup>, I. W. Harry<sup>36</sup>, C.-J. Haster<sup>90</sup>, K. Haughian<sup>44</sup>, F. J. Hayes<sup>44</sup>, J. Healy<sup>58</sup>,
A. Heidmann<sup>71</sup>, M. C. Heintze<sup>7</sup>, H. Heitmann<sup>65</sup>, P. Hello<sup>25</sup>, G. Hemming<sup>28</sup>, M. Hendry<sup>44</sup>, I. S. Heng<sup>44</sup>, J. Hennig<sup>8,9</sup>,
A. W. Heptonstall<sup>1</sup>, Francisco Hernandez Vivanco<sup>6</sup>, M. Heurs<sup>8,9</sup>, S. Hild<sup>44</sup>, T. Hinderer<sup>37,128,129</sup>, D. Hoak<sup>28</sup>, S. Hochheim<sup>8,9</sup>,
D. Hofman<sup>22</sup>, A. M. Holgado<sup>17</sup>, N. A. Holland<sup>21</sup>, K. Holt<sup>7</sup>, D. E. Holz<sup>91</sup>, P. Hopkins<sup>68</sup>, C. Horst<sup>23</sup>, J. Hough<sup>44</sup>, E. J. Howell<sup>63</sup>,
C. G. Hoy<sup>68</sup>, A. Hreibi<sup>65</sup>, E. A. Huerta<sup>17</sup>, D. Huet<sup>25</sup>, B. Hughey<sup>34</sup>, M. Hulko<sup>1</sup>, S. Husa<sup>100</sup>, S. H. Huttner<sup>44</sup>, T. Huynh-Dinh<sup>7</sup>,
B. Idzkowski<sup>74</sup>, A. Iess<sup>30,31</sup>, C. Ingram<sup>55</sup>, R. Inta<sup>84</sup>, G. Intini<sup>32,114</sup>, B. Irwin<sup>116</sup>, H. N. Isa<sup>44</sup>, J.-M. Isac<sup>71</sup>, M. Isi<sup>1</sup>, B. R. Iyer<sup>16</sup>,
K. Izumi<sup>45</sup>, T. Jacqmin<sup>71</sup>, S. J. Jadhav<sup>130</sup>, K. Jani<sup>77</sup>, N. N. Janthalur<sup>130</sup>, P. Jaranowski<sup>131</sup>, A. C. Jenkins<sup>132</sup>, J. Jiang<sup>48</sup>,
D. S. Johnson<sup>17</sup>, A. W. Jones<sup>11</sup>, D. I. Jones<sup>133</sup>, R. Jones<sup>44</sup>, R. J. G. Jonker<sup>37</sup>, L. Ju<sup>63</sup>, J. Junker<sup>8,9</sup>, C. V. Kalaghatgi<sup>68</sup>, V. Kalogera<sup>59</sup>,
B. Kamai<sup>1</sup>, S. Kandhasamy<sup>85</sup>, G. Kang<sup>38</sup>, J. B. Kanner<sup>1</sup>, S. J. Kapadia<sup>23</sup>, S. Karki<sup>70</sup>, K. S. Karvinen<sup>8,9</sup>, R. Kashyap<sup>16</sup>,
M. Kasprzack<sup>1</sup>, S. Katsanevas<sup>28</sup>, E. Katsavounidis<sup>12</sup>, W. Katzman<sup>7</sup>, S. Kaufer<sup>9</sup>, K. Kawabe<sup>45</sup>, N. V. Keerthana<sup>3</sup>, F. Kéfélian<sup>65</sup>,
D. Keitel<sup>44</sup>, R. Kennedy<sup>109</sup>, J. S. Key<sup>134</sup>, F. Y. Khalili<sup>62</sup>, H. Khan<sup>26</sup>, I. Khan<sup>14,31</sup>, S. Khan<sup>8,9</sup>, Z. Khan<sup>108</sup>, E. A. Khazanov<sup>135</sup>,
M. Khursheed<sup>61</sup>, N. Kijbunchoo<sup>21</sup>, Chunglee Kim<sup>136</sup>, J. C. Kim<sup>137</sup>, K. Kim<sup>92</sup>, W. Kim<sup>55</sup>, W. S. Kim<sup>138</sup>, Y.-M. Kim<sup>139</sup>,
C. Kimball<sup>59</sup>, E. J. King<sup>55</sup>, P. J. King<sup>45</sup>, M. Kinley-Hanlon<sup>127</sup>, R. Kirchhoff<sup>8,9</sup>, J. S. Kissel<sup>45</sup>, L. Kleybolte<sup>140</sup>, J. H. Klika<sup>23</sup>,
S. Klimenko<sup>48</sup>, T. D. Knowles<sup>39</sup>, P. Koch<sup>8,9</sup>, S. M. Koehlenbeck<sup>8,9</sup>, G. Koekoek<sup>37,141</sup>, S. Koley<sup>37</sup>, V. Kondrashov<sup>1</sup>, A. Kontos<sup>12</sup>,
N. Koper<sup>8,9</sup>, M. Korobko<sup>140</sup>, W. Z. Korth<sup>1</sup>, I. Kowalska<sup>74</sup>, D. B. Kozak<sup>1</sup>, V. Kringel<sup>8,9</sup>, N. Krishnendu<sup>29</sup>, A. Królak<sup>142,143</sup>,
G. Kuehn<sup>8,9</sup>, A. Kumar<sup>130</sup>, P. Kumar<sup>144</sup>, R. Kumar<sup>108</sup>, S. Kumar<sup>16</sup>, L. Kuo<sup>88</sup>, A. Kutynia<sup>142</sup>, S. Kwang<sup>23</sup>, B. D. Lackey<sup>36</sup>,
K. H. Lai<sup>92</sup>, T. L. Lam<sup>92</sup>, M. Landry<sup>45</sup>, B. B. Lane<sup>12</sup>, R. N. Lang<sup>145</sup>, J. Lange<sup>58</sup>, B. Lantz<sup>49</sup>, R. K. Lanza<sup>12</sup>, A. Larтаux-Vollard<sup>25</sup>,
P. D. Lasky<sup>6</sup>, M. Laxen<sup>7</sup>, A. Lazzarini<sup>1</sup>, C. Lazzaro<sup>52</sup>, P. Leaci<sup>32,114</sup>, S. Leavey<sup>8,9</sup>, Y. K. Lecoeuche<sup>45</sup>, C. H. Lee<sup>94</sup>, H. K. Lee<sup>146</sup>,
H. M. Lee<sup>147</sup>, H. W. Lee<sup>137</sup>, J. Lee<sup>93</sup>, K. Lee<sup>44</sup>, J. Lehmann<sup>8,9</sup>, A. Lenon<sup>39</sup>, N. Leroy<sup>25</sup>, N. Letendre<sup>33</sup>, Y. Levin<sup>6,99</sup>, J. Li<sup>83</sup>,
K. J. L. Li<sup>92</sup>, T. G. F. Li<sup>92</sup>, X. Li<sup>46</sup>, F. Lin<sup>6</sup>, F. Linde<sup>37</sup>, S. D. Linker<sup>112</sup>, T. B. Littenberg<sup>148</sup>, J. Liu<sup>63</sup>, X. Liu<sup>23</sup>, R. K. L. Lo<sup>1,92</sup>,
N. A. Lockerbie<sup>24</sup>, L. T. London<sup>68</sup>, A. Longo<sup>149,150</sup>, M. Lorenzini<sup>14,15</sup>, V. Loriette<sup>151</sup>, M. Lormand<sup>7</sup>, G. Losurdo<sup>19</sup>, J. D. Lough<sup>8,9</sup>,
C. O. Lousto<sup>58</sup>, G. Lovelace<sup>26</sup>, M. E. Lower<sup>152</sup>, H. Lück<sup>8,9</sup>, D. Lumaca<sup>30,31</sup>, A. P. Lundgren<sup>153</sup>, R. Lynch<sup>12</sup>, Y. Ma<sup>46</sup>, R. Macas<sup>68</sup>,
S. Macfoy<sup>24</sup>, M. MacInnis<sup>12</sup>, D. M. Macleod<sup>68</sup>, A. Macquet<sup>65</sup>, F. Magaña-Sandoval<sup>42</sup>, L. Magaña Zertuche<sup>85</sup>, R. M. Magee<sup>87</sup>,
E. Majorana<sup>32</sup>, I. Maksimovic<sup>151</sup>, A. Malik<sup>61</sup>, N. Man<sup>65</sup>, V. Mandic<sup>43</sup>, V. Mangano<sup>44</sup>, G. L. Mansell<sup>12,45</sup>, M. Manske<sup>21,23</sup>,
M. Mantovani<sup>28</sup>, M. Mapelli<sup>51,52,96,154</sup>, F. Marchesoni<sup>41,50</sup>, F. Marion<sup>33</sup>, S. Márka<sup>99</sup>, Z. Márka<sup>99</sup>, C. Markakis<sup>10,17</sup>,
A. S. Markosyan<sup>49</sup>, A. Markowitz<sup>1</sup>, E. Maros<sup>1</sup>, A. Marquina<sup>103</sup>, S. Marsat<sup>36</sup>, F. Martelli<sup>72,73</sup>, I. W. Martin<sup>44</sup>, R. M. Martin<sup>35</sup>,
D. V. Martynov<sup>11</sup>, K. Mason<sup>12</sup>, E. Massera<sup>109</sup>, A. Masserot<sup>33</sup>, T. J. Massinger<sup>1</sup>, M. Masso-Reid<sup>44</sup>, S. Mastrogiovanni<sup>32,114</sup>,
A. Matas<sup>36,43</sup>, F. Matichard<sup>1,12</sup>, L. Matone<sup>99</sup>, N. Mavalvala<sup>12</sup>, N. Mazumder<sup>69</sup>, J. J. McCann<sup>63</sup>, R. McCarthy<sup>45</sup>, D. E. McClelland<sup>21</sup>,
S. McCormick<sup>7</sup>, L. McCuller<sup>12</sup>, S. C. McGuire<sup>155</sup>, J. McIver<sup>1</sup>, D. J. McManus<sup>21</sup>, T. McRae<sup>21</sup>, S. T. McWilliams<sup>39</sup>, D. Meacher<sup>87</sup>,
G. D. Meadors<sup>6</sup>, M. Mehmet<sup>8,9</sup>, A. K. Mehta<sup>16</sup>, J. Meidam<sup>37</sup>, A. Melatos<sup>98</sup>, G. Mendell<sup>45</sup>, R. A. Mercer<sup>23</sup>, L. Mereni<sup>22</sup>,
E. L. Merilh<sup>45</sup>, M. Merzougui<sup>65</sup>, S. Meshkov<sup>1</sup>, C. Messenger<sup>44</sup>, C. Messick<sup>87</sup>, R. Metzdorff<sup>71</sup>, P. M. Meyers<sup>98</sup>, H. Miao<sup>11</sup>,
C. Michel<sup>22</sup>, H. Middleton<sup>98</sup>, E. E. Mikhailov<sup>156</sup>, L. Milano<sup>5,79</sup>, A. L. Miller<sup>48</sup>, A. Miller<sup>32,114</sup>, M. Millhouse<sup>53</sup>, J. C. Mills<sup>68</sup>,
M. C. Milovich-Goff<sup>112</sup>, O. Minazzoli<sup>65,157</sup>, Y. Minenkov<sup>31</sup>, A. Mishkin<sup>48</sup>, C. Mishra<sup>158</sup>, T. Mistry<sup>109</sup>, S. Mitra<sup>3</sup>,
V. P. Mitrofanov<sup>62</sup>, G. Mitselmakher<sup>48</sup>, R. Mittleman<sup>12</sup>, G. Mo<sup>95</sup>, D. Moffa<sup>116</sup>, K. Mogushi<sup>85</sup>, S. R. P. Mohapatra<sup>12</sup>,
M. Montani<sup>72,73</sup>, C. J. Moore<sup>10</sup>, D. Moraru<sup>45</sup>, G. Moreno<sup>45</sup>, S. Morisaki<sup>82</sup>, B. Mours<sup>33</sup>, C. M. Mow-Lowry<sup>11</sup>,
Arunava Mukherjee<sup>8,9</sup>, D. Mukherjee<sup>23</sup>, S. Mukherjee<sup>105</sup>, N. Mukund<sup>3</sup>, A. Mullavey<sup>7</sup>, J. Munch<sup>55</sup>, E. A. Muñiz<sup>42</sup>, M. Muratore<sup>34</sup>,
P. G. Murray<sup>44</sup>, A. Nagar<sup>86,159,160</sup>, I. Nardecchia<sup>30,31</sup>, L. Naticchioni<sup>32,114</sup>, R. K. Nayak<sup>161</sup>, J. Neilson<sup>112</sup>, G. Nelemans<sup>37,64</sup>,
T. J. N. Nelson<sup>7</sup>, M. Nery<sup>8,9</sup>, A. Neunzert<sup>126</sup>, K. Y. Ng<sup>12</sup>, S. Ng<sup>55</sup>, P. Nguyen<sup>70</sup>, D. Nichols<sup>37,128</sup>, S. Nissanke<sup>37,128</sup>, F. Nocera<sup>28</sup>,
C. North<sup>68</sup>, L. K. Nuttall<sup>153</sup>, M. Obergaulinger<sup>20</sup>, J. Oberling<sup>45</sup>, B. D. O'Brien<sup>48</sup>, G. D. O'Dea<sup>112</sup>, G. H. Ogin<sup>162</sup>, J. J. Oh<sup>138</sup>,
S. H. Oh<sup>138</sup>, F. Ohme<sup>8,9</sup>, H. Ohta<sup>82</sup>, M. A. Okada<sup>13</sup>, M. Oliver<sup>100</sup>, P. Oppermann<sup>8,9</sup>, Richard J. Oram<sup>7</sup>, B. O'Reilly<sup>7</sup>,
R. G. Ormiston<sup>43</sup>, L. F. Ortega<sup>48</sup>, R. O'Shaughnessy<sup>58</sup>, S. Ossokine<sup>36</sup>, D. J. Ottaway<sup>55</sup>, H. Overmier<sup>7</sup>, B. J. Owen<sup>84</sup>, A. E. Pace<sup>87</sup>,
G. Pagano<sup>18,19</sup>, M. A. Page<sup>63</sup>, A. Pai<sup>122</sup>, S. A. Pai<sup>61</sup>, J. R. Palamos<sup>70</sup>, O. Palashov<sup>135</sup>, C. Palomba<sup>32</sup>, A. Pal-Singh<sup>140</sup>,
Huang-Wei Pan<sup>88</sup>, B. Pang<sup>46</sup>, P. T. H. Pang<sup>92</sup>, C. Pankow<sup>59</sup>, F. Pannarale<sup>32,114</sup>, B. C. Pant<sup>61</sup>, F. Paoletti<sup>19</sup>, A. Paoli<sup>28</sup>, A. Parida<sup>3</sup>,
W. Parker<sup>7,155</sup>, D. Pascucci<sup>44</sup>, A. Pasqualetti<sup>28</sup>, R. Passaquieti<sup>18,19</sup>, D. Passuello<sup>19</sup>, M. Patil<sup>143</sup>, B. Patricelli<sup>18,19</sup>, B. L. Pearlstone<sup>44</sup>,
C. Pedersen<sup>68</sup>, M. Pedraza<sup>1</sup>, R. Pedurand<sup>22,163</sup>, A. Pele<sup>7</sup>, S. Penn<sup>164</sup>, C. J. Perez<sup>45</sup>, A. Perreca<sup>97,113</sup>, H. P. Pfeiffer<sup>36,90</sup>, M. Phelps<sup>8,9</sup>,
K. S. Phukon<sup>3</sup>, O. J. Piccinni<sup>32,114</sup>, M. Pichot<sup>65</sup>, F. Piergiovanni<sup>72,73</sup>, G. Pillant<sup>28</sup>, L. Pinard<sup>22</sup>, M. Pirello<sup>45</sup>, M. Pitkin<sup>44</sup>,
R. Poggiani<sup>18,19</sup>, D. Y. T. Pong<sup>92</sup>, S. Ponrathnam<sup>3</sup>, P. Popolizio<sup>28</sup>, E. K. Porter<sup>27</sup>, J. Powell<sup>152</sup>, A. K. Prajapati<sup>108</sup>, J. Prasad<sup>3</sup>,
K. Prasai<sup>49</sup>, R. Prasanna<sup>130</sup>, G. Pratten<sup>100</sup>, T. Prestegard<sup>23</sup>, S. Privitera<sup>36</sup>, G. A. Prodi<sup>97,113</sup>, L. G. Prokhorov<sup>62</sup>, O. Puncken<sup>8,9</sup>,
M. Punturo<sup>41</sup>, P. Puppo<sup>32</sup>, M. Pürrer<sup>36</sup>, H. Qi<sup>23</sup>, V. Quetschke<sup>105</sup>, P. J. Quinonez<sup>34</sup>, E. A. Quintero<sup>1</sup>, R. Quitzow-James<sup>70</sup>,
F. J. Raab<sup>45</sup>, H. Radkins<sup>45</sup>, N. Radulescu<sup>65</sup>, P. Raffai<sup>107</sup>, S. Raja<sup>61</sup>, C. Rajan<sup>61</sup>, B. Rajbhandari<sup>84</sup>, M. Rakhmanov<sup>105</sup>,
K. E. Ramirez<sup>105</sup>, A. Ramos-Buades<sup>100</sup>, Javed Rana<sup>3</sup>, K. Rao<sup>59</sup>, P. Rapagnani<sup>32,114</sup>, V. Raymond<sup>68</sup>, M. Razzano<sup>18,19</sup>, J. Read<sup>26</sup>,
T. Regimbau<sup>33</sup>, L. Rei<sup>60</sup>, S. Reid<sup>24</sup>, D. H. Reitz<sup>1,48</sup>, W. Ren<sup>17</sup>, F. Ricci<sup>32,114</sup>, C. J. Richardson<sup>34</sup>, J. W. Richardson<sup>1</sup>, P. M. Ricker<sup>17</sup>,
K. Riles<sup>126</sup>, M. Rizzo<sup>59</sup>, N. A. Robertson<sup>1,44</sup>, R. Robie<sup>44</sup>, F. Robinet<sup>25</sup>, A. Rocchi<sup>31</sup>, L. Rolland<sup>33</sup>, J. G. Rollins<sup>1</sup>, V. J. Roma<sup>70</sup>,
M. Romanelli<sup>67</sup>, R. Romano<sup>4,5</sup>, C. L. Romel<sup>45</sup>, J. H. Romie<sup>7</sup>, K. Rose<sup>116</sup>, D. Rosińska<sup>54,165</sup>, S. G. Rosofsky<sup>17</sup>, M. P. Ross<sup>166</sup>,
S. Rowan<sup>44</sup>, A. Rüdiger<sup>8,9,183</sup>, P. Ruggi<sup>28</sup>, G. Rutins<sup>167</sup>, K. Ryan<sup>45</sup>, S. Sachdev<sup>1</sup>, T. Sadecki<sup>45</sup>, M. Sakellariadou<sup>132</sup>,
O. S. Salafia<sup>168,169,170</sup>, L. Salconi<sup>28</sup>, M. Saleem<sup>29</sup>, A. Samajdar<sup>37</sup>, L. Sammut<sup>6</sup>, E. J. Sanchez<sup>1</sup>, L. E. Sanchez<sup>1</sup>, N. Sanchis-Gual<sup>20</sup>,
V. Sandberg<sup>45</sup>, J. R. Sanders<sup>42</sup>, K. A. Santiago<sup>35</sup>, N. Sarin<sup>6</sup>, B. Sassolas<sup>22</sup>, B. S. Sathyaprakash<sup>68,87</sup>, P. R. Saulson<sup>42</sup>, O. Sauter<sup>126</sup>, R. L. Savage<sup>45</sup>, P. Schale<sup>70</sup>, M. Scheel<sup>46</sup>, J. Scheuer<sup>59</sup>, P. Schmidt<sup>64</sup>, R. Schnabel<sup>140</sup>, R. M. S. Schofield<sup>70</sup>, A. Schönbeck<sup>140</sup>, E. Schreiber<sup>8,9</sup>, B. W. Schulte<sup>8,9</sup>, B. F. Schutz<sup>68</sup>, S. G. Schwalbe<sup>34</sup>, J. Scott<sup>44</sup>, S. M. Scott<sup>21</sup>, E. Seidel<sup>17</sup>, D. Sellers<sup>7</sup>, A. S. Sengupta<sup>171</sup>, N. Sennett<sup>36</sup>, D. Sentenac<sup>28</sup>, V. Sequino<sup>14,30,31</sup>, A. Sergeev<sup>135</sup>, Y. Setyawati<sup>8,9</sup>, D. A. Shaddock<sup>21</sup>, T. Shaffer<sup>45</sup>, M. S. Shahriar<sup>59</sup>, M. B. Shaner<sup>112</sup>, L. Shao<sup>36</sup>, P. Sharma<sup>61</sup>, P. Shawhan<sup>76</sup>, H. Shen<sup>17</sup>, R. Shink<sup>172</sup>, D. H. Shoemaker<sup>12</sup>, D. M. Shoemaker<sup>77</sup>, S. ShyamSundar<sup>61</sup>, K. Siellez<sup>77</sup>, M. Sieniawska<sup>54</sup>, D. Sigg<sup>15</sup>, A. D. Silva<sup>13</sup>, L. P. Singer<sup>80</sup>, N. Singh<sup>74</sup>, A. Singhal<sup>14,32</sup>, A. M. Sintes<sup>100</sup>, S. Sitmukhambetov<sup>105</sup>, V. Skliris<sup>68</sup>, B. J. J. Slagmolen<sup>21</sup>, T. J. Slaven-Blair<sup>63</sup>, J. R. Smith<sup>26</sup>, R. J. E. Smith<sup>6</sup>, S. Somala<sup>173</sup>, E. J. Son<sup>138</sup>, B. Sorazu<sup>44</sup>, F. Sorrentino<sup>60</sup>, T. Souradeep<sup>3</sup>, E. Sowell<sup>84</sup>, A. P. Spencer<sup>44</sup>, M. Spera<sup>51,52,59,96,154</sup>, A. K. Srivastava<sup>108</sup>, V. Srivastava<sup>42</sup>, K. Staats<sup>59</sup>, C. Stachie<sup>65</sup>, M. Standke<sup>8,9</sup>, D. A. Steer<sup>27</sup>, M. Steinke<sup>8,9</sup>, J. Steinlechner<sup>44,140</sup>, S. Steinlechner<sup>140</sup>, D. Steinmeyer<sup>8,9</sup>, S. P. Stevenson<sup>152</sup>, D. Stocks<sup>49</sup>, R. Stone<sup>105</sup>, D. J. Stops<sup>11</sup>, K. A. Strain<sup>44</sup>, G. Stratta<sup>72,73</sup>, S. E. Strigin<sup>62</sup>, A. Strunk<sup>45</sup>, R. Sturani<sup>174</sup>, A. L. Stuver<sup>175</sup>, V. Sudhir<sup>12</sup>, T. Z. Summerscales<sup>176</sup>, L. Sun<sup>1</sup>, S. Sunil<sup>108</sup>, J. Suresh<sup>3</sup>, P. J. Sutton<sup>68</sup>, B. L. Swinkels<sup>37</sup>, M. J. Szczepańczyk<sup>34</sup>, M. Tacca<sup>37</sup>, S. C. Tait<sup>44</sup>, C. Talbot<sup>6</sup>, D. Talukder<sup>70</sup>, D. B. Tanner<sup>48</sup>, M. Tápai<sup>123</sup>, A. Taracchini<sup>36</sup>, J. D. Tasson<sup>95</sup>, R. Taylor<sup>1</sup>, F. Thies<sup>8,9</sup>, M. Thomas<sup>7</sup>, P. Thomas<sup>45</sup>, S. R. Thondapu<sup>61</sup>, K. A. Thorne<sup>7</sup>, E. Thrane<sup>6</sup>, Shubhanshu Tiwari<sup>97,113</sup>, Srishti Tiwari<sup>124</sup>, V. Tiwari<sup>68</sup>, K. Toland<sup>44</sup>, M. Tonelli<sup>18,19</sup>, Z. Tornasi<sup>44</sup>, A. Torres-Forné<sup>177</sup>, C. I. Torrie<sup>1</sup>, D. Töyrä<sup>11</sup>, F. Travasso<sup>28,41</sup>, G. Traylor<sup>7</sup>, M. C. Tringali<sup>74</sup>, A. Trovato<sup>27</sup>, L. Trozzo<sup>19,178</sup>, R. Trudeau<sup>1</sup>, K. W. Tsang<sup>37</sup>, M. Tse<sup>12</sup>, R. Tso<sup>46</sup>, L. Tsukada<sup>82</sup>, D. Tsuna<sup>82</sup>, D. Tuyenbayev<sup>105</sup>, K. Ueno<sup>82</sup>, D. Ugolini<sup>179</sup>, C. S. Unnikrishnan<sup>124</sup>, A. L. Urban<sup>2</sup>, S. A. Usman<sup>68</sup>, H. Vahlbruch<sup>9</sup>, G. Vajente<sup>1</sup>, G. Valdes<sup>2</sup>, N. van Bakel<sup>37</sup>, M. van Beuzekom<sup>37</sup>, J. F. J. van den Brand<sup>37,75</sup>, C. Van Den Broeck<sup>37,180</sup>, D. C. Vander-Hyde<sup>42</sup>, L. van der Schaaf<sup>37</sup>, J. V. van Heijningen<sup>63</sup>, A. A. van Veggel<sup>44</sup>, M. Vardaro<sup>51,52</sup>, V. Varma<sup>46</sup>, S. Vass<sup>1</sup>, M. Vasúth<sup>47</sup>, A. Vecchio<sup>11</sup>, G. Vedovato<sup>52</sup>, J. Veitch<sup>44</sup>, P. J. Veitch<sup>55</sup>, K. Venkateswara<sup>166</sup>, G. Venugopalan<sup>1</sup>, D. Verkindt<sup>33</sup>, F. Vetrano<sup>72,73</sup>, A. Viceré<sup>72,73</sup>, A. D. Viets<sup>23</sup>, D. J. Vine<sup>167</sup>, J.-Y. Vinet<sup>65</sup>, S. Vitale<sup>12</sup>, T. Vo<sup>42</sup>, H. Vocca<sup>40,41</sup>, C. Vorvick<sup>45</sup>, S. P. Vyatchanin<sup>62</sup>, A. R. Wade<sup>1</sup>, L. E. Wade<sup>116</sup>, M. Wade<sup>116</sup>, R. Walet<sup>37</sup>, M. Walker<sup>26</sup>, L. Wallace<sup>1</sup>, S. Walsh<sup>23</sup>, G. Wang<sup>14,19</sup>, H. Wang<sup>11</sup>, J. Z. Wang<sup>126</sup>, W. H. Wang<sup>105</sup>, Y. F. Wang<sup>92</sup>, R. L. Ward<sup>21</sup>, Z. A. Warden<sup>34</sup>, J. Warner<sup>45</sup>, M. Was<sup>33</sup>, J. Watchi<sup>101</sup>, B. Weaver<sup>45</sup>, L.-W. Wei<sup>8,9</sup>, M. Weinert<sup>8,9</sup>, A. J. Weinstein<sup>1</sup>, R. Weiss<sup>12</sup>, F. Wellmann<sup>8,9</sup>, L. Wen<sup>63</sup>, E. K. Wessel<sup>17</sup>, P. Weßels<sup>8,9</sup>, J. W. Westhouse<sup>34</sup>, K. Wette<sup>21</sup>, J. T. Whelan<sup>58</sup>, B. F. Whiting<sup>48</sup>, C. Whittle<sup>12</sup>, D. M. Wilken<sup>8,9</sup>, D. Williams<sup>44</sup>, A. R. Williamson<sup>37,128</sup>, J. L. Willis<sup>1</sup>, B. Willke<sup>8,9</sup>, M. H. Wimmer<sup>8,9</sup>, W. Winkler<sup>8,9</sup>, C. C. Wipf<sup>1</sup>, H. Wittel<sup>8,9</sup>, G. Woan<sup>44</sup>, J. Woehler<sup>8,9</sup>, J. K. Wofford<sup>58</sup>, J. Worden<sup>45</sup>, J. L. Wright<sup>44</sup>, D. S. Wu<sup>8,9</sup>, D. M. Wysocki<sup>58</sup>, L. Xiao<sup>1</sup>, H. Yamamoto<sup>1</sup>, C. C. Yancey<sup>76</sup>, L. Yang<sup>115</sup>, M. J. Yap<sup>21</sup>, M. Yazback<sup>48</sup>, D. W. Yeeles<sup>68</sup>, Hang Yu<sup>12</sup>, Haocun Yu<sup>12</sup>, S. H. R. Yuen<sup>92</sup>, M. Yvert<sup>33</sup>, A. K. Zadrożny<sup>105,142</sup>, M. Zanolin<sup>34</sup>, T. Zelenova<sup>28</sup>, J.-P. Zendri<sup>52</sup>, M. Zevin<sup>59</sup>, J. Zhang<sup>63</sup>, L. Zhang<sup>1</sup>, T. Zhang<sup>44</sup>, C. Zhao<sup>63</sup>, M. Zhou<sup>59</sup>, Z. Zhou<sup>59</sup>, X. J. Zhu<sup>6</sup>, A. B. Zimmerman<sup>90</sup>, Y. Zlochower<sup>58</sup>, M. E. Zucker<sup>1,12</sup>, and J. Zweizig<sup>1</sup>

The LIGO Scientific Collaboration and the Virgo Collaboration

<footnote n="1"><sup>1</sup> LIGO, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="2"><sup>2</sup>Louisiana State University, Baton Rouge, LA 70803, USA</footnote>

<footnote n="3">³ Inter-University Centre for Astronomy and Astrophysics, Pune 411007, India</footnote>

<footnote n="4">⁴ Università di Salerno, Fisciano, I-84084 Salerno, Italy</footnote>

<footnote n="5">⁵ INFN, Sezione di Napoli, Complesso Universitario di Monte S. Angelo, I-80126 Napoli, Italy</footnote>

<footnote n="6">⁶ OzGrav, School of Physics & Astronomy, Monash University, Clayton 3800, Victoria, Australia</footnote>

<footnote n="7">⁷ LIGO Livingston Observatory, Livingston, LA 70754, USA</footnote>

<footnote n="8">⁸ Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany</footnote>

<footnote n="9">⁹ Leibniz Universität Hannover, D-30167 Hannover, Germany</footnote>

<footnote n="10">¹⁰ University of Cambridge, Cambridge CB2 1TN, UK</footnote>

<footnote n="11">¹¹ University of Birmingham, Birmingham B15 2TT, UK</footnote>

<footnote n="12">¹² LIGO, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="13">¹³ Instituto Nacional de Pesquisas Espaciais, 12227-010 São José dos Campos, São Paulo, Brazil</footnote>

<footnote n="14">¹⁴ Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy</footnote>

<footnote n="15">¹⁵ INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy</footnote>

<footnote n="16">¹⁶ International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560089, India</footnote>

<footnote n="17">¹⁷ NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA</footnote>

<footnote n="18">¹⁸ Università di Pisa, I-56127 Pisa, Italy</footnote>

<footnote n="19">¹⁹ INFN, Sezione di Pisa, I-56127 Pisa, Italy</footnote>

<footnote n="20">²⁰ Departamento de Astronomía y Astrofísica, Universitat de València, E-46100 Burjassot, València, Spain</footnote>

<footnote n="21">²¹ OzGrav, Australian National University, Canberra, Australian Capital Territory 0200, Australia</footnote>

<footnote n="22">²² Laboratoire des Matériaux Avancés (LMA), CNRS/IN2P3, F-69622 Villeurbanne, France</footnote>

<footnote n="23">²³ University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA</footnote>

<footnote n="24">²⁴ SUPA, University of Strathclyde, Glasgow G1 XQ, UK</footnote>

<footnote n="25">²⁵ LAL, Univ. Paris-Sud, CNRS/IN2P3, Université Paris-Saclay, F-91898 Orsay, France</footnote>

<footnote n="26">²⁶ California State University Fullerton, Fullerton, CA 92831, USA</footnote>

<footnote n="27">²⁷ APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cité,
F-75205 Paris Cedex 13, France</footnote>

<footnote n="28">²⁸ European Gravitational Observatory (EGO), I-56021 Cascina, Pisa, Italy</footnote>

<footnote n="29">²⁹Chennai Mathematical Institute, Chennai 603103, India</footnote>

<footnote n="30">³⁰ Università di Roma Tor Vergata, I-00133 Roma, Italy</footnote>

<footnote n="31">³¹ INFN, Sezione di Roma Tor Vergata, I-00133 Roma, Italy</footnote>

<footnote n="32">³²INFN, Sezione di Roma, I-00185 Roma, Italy</footnote>

<footnote n="33">³³Laboratoire d'Annecy de Physique des Particules (LAPP), Univ. Grenoble Alpes, Université Savoie Mont Blanc, CNRS/IN2P3, F-74941 Annecy, France</footnote>

<footnote n="34">³⁴Embry-Riddle Aeronautical University, Prescott, AZ 86301, USA</footnote>

<footnote n="35">³⁵ Montclair State University, Montclair, NJ 07043, USA</footnote>
³⁶ Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam-Golm, Germany

³⁷ Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands

³⁸ Korea Institute of Science and Technology Information, Daejeon 34141, Republic of Korea

³⁹ West Virginia University, Morgantown, WV 26506, USA

⁴⁰ Università di Perugia, I-06123 Perugia, Italy

⁴¹ INFN, Sezione di Perugia, I-06123 Perugia, Italy

⁴² Syracuse University, Syracuse, NY 13244, USA

⁴³ University of Minnesota, Minneapolis, MN 55455, USA

⁴⁴ SUPA, University of Glasgow, Glasgow G12 8QQ, UK

⁴⁵ LIGO Hanford Observatory, Richland, WA 99352, USA

⁴⁶ Caltech CaRT, Pasadena, CA 91125, USA

⁴⁷ Wigner RCP, RMKI, H-1121 Budapest, Konkoly Thege Miklós út 29-33, Hungary

⁴⁸ University of Florida, Gainesville, FL 32611, USA

⁴⁹ Stanford University, Stanford, CA 94305, USA

⁵⁰ Università di Camerino, Dipartimento di Fisica, I-62032 Camerino, Italy

⁵¹ Università di Padova, Dipartimento di Fisica e Astronomia, I-35131 Padova, Italy

⁵² INFN, Sezione di Padova, I-35131 Padova, Italy

⁵³ Montana State University, Bozeman, MT 59717, USA

⁵⁴ Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, 00-716, Warsaw, Poland

⁵⁵ OzGrav, University of Adelaide, Adelaide, South Australia 5005, Australia

⁵⁶ Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität Jena, D-07743 Jena, Germany

⁵⁷ INFN, Sezione di Milano Bicocca, Gruppo Collegato di Parma, I-43124 Parma, Italy

⁵⁸ Rochester Institute of Technology, Rochester, NY 14623, USA

⁵⁹ Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA

⁶⁰ INFN, Sezione di Genova, I-16146 Genova, Italy

⁶¹ RRCAT, Indore, Madhya Pradesh 452013, India

⁶² Faculty of Physics, Lomonosov Moscow State University, Moscow 119991, Russia

⁶³ OzGrav, University of Western Australia, Crawley, Western Australia 6009, Australia

⁶⁴ Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands

⁶⁵ Artemis, Université Côte d'Azur, Observatoire Côte d'Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France

⁶⁶ Physik-Institut, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

⁶⁷ Univ Rennes, CNRS, Institut FOTON—UMR6082, F-3500 Rennes, France

⁶⁸ Cardiff University, Cardiff CF24 3AA, UK

⁶⁹ Washington State University, Pullman, WA 99164, USA

⁷⁰ University of Oregon, Eugene, OR 97403, USA

⁷¹ Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-Université PSL, Collège de France, F-75005 Paris, France

⁷² Università degli Studi di Urbino “Carlo Bo,” I-61029 Urbino, Italy

⁷³ INFN, Sezione di Firenze, I-50019 Sesto Fiorentino, Firenze, Italy

⁷⁴ Astronomical Observatory Warsaw University, 00-478 Warsaw, Poland

⁷⁵ VU University Amsterdam, 1081 HV Amsterdam, The Netherlands

⁷⁶ University of Maryland, College Park, MD 20742, USA

⁷⁷ School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA

⁷⁸ Université Claude Bernard Lyon 1, F-69622 Villeurbanne, France

⁷⁹ Università di Napoli “Federico II,” Complesso Universitario di Monte S.Angelo, I-80126 Napoli, Italy

⁸⁰ NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

⁸¹ Dipartimento di Fisica, Università degli Studi di Genova, I-16146 Genova, Italy

⁸² RESCEU, University of Tokyo, Tokyo, 113-0033, Japan

⁸³ Tsinghua University, Beijing 100084, People's Republic of China

⁸⁴ Texas Tech University, Lubbock, TX 79409, USA

⁸⁵ The University of Mississippi, University, MS 38677, USA

⁸⁶ Museo Storico della Fisica e Centro Studi and Ricerche “Enrico Fermi,” I-00184 Roma, Italy

⁸⁷ The Pennsylvania State University, University Park, PA 16802, USA

⁸⁸ National Tsing Hua University, Hsinchu City, 30013 Taiwan, Republic of China

⁸⁹ Charles Sturt University, Wagga Wagga, New South Wales 2678, Australia

⁹⁰ Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, ON M5S 3H8, Canada

⁹¹ University of Chicago, Chicago, IL 60637, USA

⁹² The Chinese University of Hong Kong, Shatin, NT, Hong Kong

⁹³ Seoul National University, Seoul 08826, Republic of Korea

⁹⁴ Pusan National University, Busan 46241, Republic of Korea

⁹⁵ Carleton College, Northfield, MN 55057, USA

⁹⁶ INAF, Osservatorio Astronomico di Padova, I-35122 Padova, Italy

⁹⁷ INFN, Trento Institute for Fundamental Physics and Applications, I-38123 Povo, Trento, Italy

⁹⁸ OzGrav, University of Melbourne, Parkville, Victoria 3010, Australia

⁹⁹ Columbia University, New York, NY 10027, USA

¹⁰⁰ Universitat de les Illes Balears, IAC3—IEEC, E-07122 Palma de Mallorca, Spain

¹⁰¹ Université Libre de Bruxelles, Brussels B-1050, Belgium

¹⁰² Sonoma State University, Rohnert Park, CA 94928, USA

¹⁰³ Departamento de Matemáticas, Universitat de València, E-46100 Burjassot, València, Spain

¹⁰⁴ University of Rhode Island, Kingston, RI 02881, USA

¹⁰⁵ The University of Texas Rio Grande Valley, Brownsville, TX 78520, USA

¹⁰⁶ Bellevue College, Bellevue, WA 98007, USA

¹⁰⁷ MTA-ELTE Astrophysics Research Group, Institute of Physics, Eötvös University, Budapest 1117, Hungary

¹⁰⁸ Institute for Plasma Research, Bhat, Gandhinagar 382428, India

¹⁰⁹ The University of Sheffield, Sheffield S10 2TN, UK

¹¹⁰ IGFAE, Campus Sur, Universidade de Santiago de Compostela, E-15782, Spain

¹¹¹ Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Università di Parma, I-43124 Parma, Italy
112 California State University, Los Angeles, 5151 State University Dr., Los Angeles, CA 90032, USA

113 Università di Trento, Dipartimento di Fisica, I-38123 Povo, Trento, Italy

114 Università di Roma “La Sapienza,” I-00185 Roma, Italy

115 Colorado State University, Fort Collins, CO 80523, USA

116 Kenyon College, Gambier, OH 43022, USA

117 Christopher Newport University, Newport News, VA 23606, USA

118 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

119 Observatori Astronòmic, Universitat de València, E-46980 Paterna, València, Spain

120 School of Mathematics, University of Edinburgh, Edinburgh EH9 3FD, UK

121 Institute of Advanced Research, Gandhinagar 382426, India

122 Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India

123 University of Szeged, Döm tér 9, Szeged 6720, Hungary

124 Tata Institute of Fundamental Research, Mumbai 400005, India

125 INAF, Osservatorio Astronomico di Capodimonte, I-80131, Napoli, Italy

126 University of Michigan, Ann Arbor, MI 48109, USA

127 American University, Washington, DC 20016, USA

128 GRAPPA, Anton Pannekoek Institute for Astronomy and Institute of High-Energy Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands

129 Delta Institute for Theoretical Physics, Science Park 904, 1090 GL Amsterdam, The Netherlands

130 Directorate of Construction, Services & Estate Management, Mumbai 400094, India

131 University of Bialystok, 15-424 Bialystok, Poland

132 King's College London, University of London, London WC2R 2LS, UK

133 University of Southampton, Southampton SO17 1BJ, UK

134 University of Washington Bothell, Bothell, WA 98011, USA

135 Institute of Applied Physics, Nizhny Novgorod, 603950, Russia

136 Ewha Womans University, Seoul 03760, Republic of Korea

137 Inje University Gimhae, South Gyeongsang 50834, Republic of Korea

138 National Institute for Mathematical Sciences, Daejeon 34047, Republic of Korea

139 Ulsan National Institute of Science and Technology, Ulsan 44919, Republic of Korea

140 Universität Hamburg, D-22761 Hamburg, Germany

141 Maastricht University, P.O. Box 616, 6200 MD Maastricht, The Netherlands

142 NCBJ, 05-400 Świerk-Otwock, Poland

143 Institute of Mathematics, Polish Academy of Sciences, 00656 Warsaw, Poland

144 Cornell University, Ithaca, NY 14850, USA

145 Hillsdale College, Hillsdale, MI 49242, USA

146 Hanyang University, Seoul 04763, Republic of Korea

147 Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea

148 NASA Marshall Space Flight Center, Huntsville, AL 35811, USA

149 Dipartimento di Matematica e Fisica, Università degli Studi Roma Tre, I-00146 Roma, Italy

150 INFN, Sezione di Roma Tre, I-00146 Roma, Italy

151 ESPCI, CNRS, F-75005 Paris, France

152 OzGrav, Swinburne University of Technology, Hawthorn VIC 3122, Australia

153 University of Portsmouth, Portsmouth, PO1 3FX, UK

154 Institut für Astro- und Teilchenphysik, Universität Innsbruck, Technikerstrasse 25/8, A-6020, Innsbruck, Austria

155 Southern University and A&M College, Baton Rouge, LA 70813, USA

156 College of William and Mary, Williamsburg, VA 23187, USA

157 Centre Scientifique de Monaco, 8 quai Antoine Ier, MC-98000, Monaco

158 Indian Institute of Technology Madras, Chennai 600036, India

159 INFN Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy

160 Institut des Hautes Etudes Scientifiques, F-91440 Bures-sur-Yvette, France

161 IISER-Kolkata, Mohanpur, West Bengal 741252, India

162 Whitman College, 345 Boyer Avenue, Walla Walla, WA 99362, USA

163 Université de Lyon, F-69361 Lyon, France

164 Hobart and William Smith Colleges, Geneva, NY 14456, USA

165 Janusz Gil Institute of Astronomy, University of Zielona Góra, 65-265 Zielona Góra, Poland

166 University of Washington, Seattle, WA 98195, USA

167 SUPA, University of the West of Scotland, Paisley PA1 2BE, UK

168 NAF, Osservatorio Astronomico di Brera sede di Merate, I-23807 Merate, Lecco, Italy

169 Università degli Studi di Milano-Bicocca, I-20126 Milano, Italy

170 INFN, Sezione di Milano-Bicocca, I-20126 Milano, Italy

171 Indian Institute of Technology, Gandhinagar Ahmedabad Gujarat 382424, India

172 Université de Montréal/Polytechnique, Montreal, QC H3T 1J4, Canada

173 Indian Institute of Technology Hyderabad, Sangareddy, Khandi, Telangana 502285, India

174 International Institute of Physics, Universidade Federal do Rio Grande do Norte, Natal RN 59078-970, Brazil

175 Villanova University, 800 Lancaster Ave., Villanova, PA 19085, USA

176 Andrews University, Berrien Springs, MI 49104, USA

177 Max Planck Institute for Gravitationalphysik (Albert Einstein Institute), D-14476 Potsdam-Golm, Germany

178 Università di Siena, I-53100 Siena, Italy

179 Trinity University, San Antonio, TX 78212, USA

180 Van Swinderen Institute for Particle Physics and Gravity, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands
*Received 2018 December 15; revised 2019 July 11; accepted 2019 July 21; published 2019 September 9*

<footnote n="181">181 Deceased, 2018 February.</footnote>

<footnote n="182">182 Deceased, 2017 November.</footnote>

<footnote n="183">183 Deceased, 2018 July.</footnote>
Abstract

We present results on the mass, spin, and redshift distributions with phenomenological population models using the 10 binary black hole (BBH) mergers detected in the first and second observing runs completed by Advanced LIGO and Advanced Virgo. We constrain properties of the BBH mass spectrum using models with a range of parameterizations of the BBH mass and spin distributions. We find that the mass distribution of the more massive BH in such binaries is well approximated by models with no more than 1% of BHs more massive than 45 $M_\odot$ and a power-law index of $\alpha = 1.3_{-1.7}^{+1.4}$ (90% credibility). We also show that BBHs are unlikely to be composed of BHs with large spins aligned to the orbital angular momentum. Modeling the evolution of the BBH merger rate with redshift, we show that it is flat or increasing with redshift with 93% probability. Marginalizing over uncertainties in the BBH population, we find robust estimates of the BBH merger rate density of $R = 53.2_{-28.2}^{+55.8} \text{ Gpc}^{-3} \text{ yr}^{-1}$ (90% credibility). As the BBH catalog grows in future observing runs, we expect that uncertainties in the population model parameters will shrink, potentially providing insights into the formation of BHs via supernovae, binary interactions of massive stars, stellar cluster dynamics, and the formation history of BHs across cosmic time.

Key words: black holes – gravitational waves – statistical

1. Introduction

The second LIGO/Virgo observing run (O2) spanned 9 months between 2016 November and 2017 August, building on the first, 4-month run (O1) in 2015. The LIGO/Virgo gravitational-wave (GW) interferometer network is composed of two instruments in the United States (LIGO; LIGO Scientific Collaboration et al. 2015; Abbott et al. 2016a) and a third in Europe (Virgo; Acernese et al. 2015), the latter joining the run in the summer of 2017. In total, 10 binary black hole (BBH) mergers have been detected to date (Abbott et al. 2018a). The BBHs detected possess a wide range of physical properties. The lightest so far is GW170608 (Abbott et al. 2017a), with an inferred total mass of $18.7_{-0.7}^{+3.3} M_\odot$. GW170729 (Abbott et al. 2018a)—exceptional in several ways—is likely to be the heaviest BBH to date, having total mass $85.2_{-11.2}^{+15.4} M_\odot$, as well as the most distant, at redshift $0.48_{-0.20}^{+0.19}$. Both GW151226 and GW170729 show evidence for at least one BH with a spin greater than zero (Abbott et al. 2016b, 2018a).

By measuring the distributions of mass, spin, and merger redshift in the BBH population, we may make inferences about the physics of binary mergers and better understand the origin of these systems. We employ Bayesian inference and modeling (Gelman et al. 2004; Mandel 2010; Foreman-Mackey et al. 2014; Hilbe et al. 2017; Asensio Ramos 2018), which, when applied to parameterized models of the population, is able to infer population-level parameters—sometimes called hyperparameters to distinguish them from the event-level parameters—while properly accounting for the uncertainty in the measurements of each event's parameters (Hogg et al. 2010; Mandel 2010).

The structure and parameterization of BBH population models are guided by the physical processes and evolutionary environments in which BBHs are expected to form and merge. Several BBH formation channels have been proposed in the literature, each of them involving a specific environment and a number of physical processes. For example, BBHs might form from isolated massive binaries in the galactic field through common-envelope evolution (Bethe & Brown 1998; Portegies Zwart & Yungelson 1998; Belczynski et al. 2002, 2007, 2008, 2014; Voss & Tauris 2003; Dewi et al. 2006; Dominik et al. 2013; Mennekens & Vanbeveren 2014; Spera et al. 2015; Eldridge & Stanway 2016; Mapelli et al. 2017; Stevenson et al. 2017b; Tauris et al. 2017; Chruslinska et al. 2018; Giacobbo & Mapelli 2018; Giacobbo et al. 2018; Kruckow et al. 2018; Mapelli &

Giacobbo 2018) or via chemically homogeneous evolution (de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016). Alternatively, BBHs might form via dynamical processes in stellar clusters (Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Grindlay et al. 2006; O'Leary et al. 2006; Ivanova et al. 2008; Sadowski et al. 2008; Downing et al. 2010, 2011; Clausen et al. 2013; Ziosi et al. 2014; Rodriguez et al. 2015, 2016a; Mapelli 2016; Askar et al. 2017; Banerjee 2017; Chatterjee et al. 2017) and galactic nuclei (Antonini & Perets 2012; Antonini & Rasio 2016; Petrovich & Antonini 2017), evolution of hierarchical triple systems (Antonini et al. 2014, 2017; Kimpson et al. 2016; Liu & Lai 2018), gas drag and stellar scattering in accretion disks surrounding supermassive BHs (McKernan et al. 2012; Bartos et al. 2017; Stone et al. 2017). Finally, BBHs might originate as part of a primordial BH population in the early universe (Carr & Hawking 1974; Bird et al. 2016; Carr et al. 2016; Inayoshi et al. 2016; Sasaki et al. 2016; Ali-Haïmoud et al. 2017; Clesse & García-Bellido 2017; Inomata et al. 2017; Ando et al. 2018; Chen & Huang 2018), where their mass spectrum is typically proposed as having power-law behavior, but spanning a much wider range of masses than stellar-mass BHs. Each channel contributes differently to the distributions of the mass, spin, distance, and orbital characteristics of BBHs.

There are several processes common to most pathways through stellar evolution that affect the properties of the resultant BBH system. Examples include mass loss (Vink et al. 2001; Vink & de Koter 2005; Gräfener & Hamann 2008) and supernovae (O'Connor & Ott 2011; Fryer et al. 2012; Janka 2012; Ugliano et al. 2012; Ertl et al. 2016; Sukhbold et al. 2016). The mass of the compact object left after the supernova is directly related to its pre-supernova mass and the supernova mechanism itself. Metallicity has been shown (Kudritzki & Puls 2000; Vink et al. 2001; Brott et al. 2011) to have important effects on stellar mass loss through winds—line-driven winds are quenched in metal-poor progenitors, enabling large BHs to form through direct collapse or post-supernova mass fallback (Heger et al. 2003; Mapelli et al. 2009; Belczynski et al. 2010; Spera et al. 2015). This also, in turn, might suppress supernova kicks (Fryer et al. 2012) and hence enhance the number of binaries that are not disrupted.

Theoretical and phenomenological models of BBH formation are explored by population synthesis. This requires modeling of not only stellar evolution but also the influence
of their evolutionary environments. For instance, isolated evolution in galactic fields requires prescriptions for binary interactions, such as common-envelope physics, as well as mass transfer episodes (see reviews in Kalogera et al. 2007; Vanbeveren 2009; Postnov & Yungelson 2014) and, more recently, the effects of rapid rotation (de Mink et al. 2009; Mandel & de Mink 2016; Marchant et al. 2016). Meanwhile, BBH formation in dense stellar clusters (Ziosi et al. 2014; Rodriguez et al. 2015, 2016a; Mapelli 2016; Askar et al. 2017; Banerjee 2017) is impacted primarily by dynamical interactions within the cluster (Fregeau 2004; Morscher et al. 2013), but also by cluster size and initial mass functions (Scheepmaker et al. 2007; Portegies Zwart et al. 2010; Kremer et al. 2019). GW observations provide an alternative to sharpen our understanding of those processes.

Electromagnetic observations and modeling of systems containing BHs have led to speculation about the existence of potential gaps in the BH mass spectrum. Both gaps may be probed using data from current ground-based GW interferometers and as such have been the target of parametric studies. At low masses, observations of X-ray binaries (XRBs) combined via Bayesian population modeling (Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011b) suggest a minimum BH mass well above the largest neutron star masses. While the existence and nature of this gap are still uncertain (Kreidberg et al. 2012), it is proposed to exist between the most massive neutron stars (Freire et al. 2008; Özel & Freire 2016; Margalit & Metzger 2017; 2.1–2.5 $M_\odot$) and the lightest BHs ($\sim 5 M_\odot$). It is possible to constrain the existence of this lower-mass gap with GW observations (Litenberg et al. 2015; Mandel et al. 2015, 2017; Kovetz et al. 2017). In Section 3, we find that our current GW observations do not inform the upper edge of this gap, inferring a minimum mass on the primary BH at $m_{\min} \lesssim 9 M_\odot$. Our volumetric sensitivity to BBH systems with masses less than 5 $M_\odot$ is small enough that we expect (and observe) no events in the lower gap region. Thus, our ability to place constraints in this region is severely limited.

Recently, there have been claims of an upper cutoff in the BBH mass spectrum based on the first few LIGO detections (Fishbach & Holz 2017; Bai et al. 2018; Talbot & Thrane 2018; Wysocki et al. 2018; Roulet & Zaldarriaga 2019). This might be expected as a consequence of a different supernova type, called the (pulsational) pair-instability supernova (Heger & Woosley 2002; Belczynski et al. 2016a; Spera & Mapelli 2017; Woosley 2017; Marchant et al. 2018). Evolved stars with a helium core mass $\gtrsim 30 M_\odot$ are expected to become unstable because efficient pair production softens their equation of state. For helium core mass $\sim 30–64 M_\odot$, the star undergoes a sequence of pulsations, losing mass until stability is reestablished (Woosley et al. 2007). The enhanced mass loss during pulsational pair instability is expected to affect the final collapse of the star, leading to smaller BH masses. The fate of a star with He core mass $\sim 64–135 M_\odot$ is more dramatic: the entire star is disrupted by a pair-instability supernova, leaving no remnant (Fowler & Hoyle 1964; Barkat et al. 1967; Rakavy & Shaviv 1967). From the combination of pair instability and pulsational pair instability, it is expected that pair-instability supernovae should leave no BH remnants between $\sim 50$ and 150 $M_\odot$ because the progenitor star is partially or entirely disrupted by the explosion. It is also possible that contributions from the merger of previous merger products—

second-generation mergers (O’Leary et al. 2016; Fishbach et al. 2017; Gerosa & Berti 2017; Rodriguez et al. 2018b)—could occupy this gap. Primordial BHs could also span numerous decades of the mass spectrum (Georg & Watson 2017), but their number density in either mass gap is dependent on the behavior of fluctuations in the early universe (Byrnes et al. 2018). Nonetheless, consistent with prior work, we find that all our mass models have almost no merging BHs above $\sim 45 M_\odot$.

Observational constraints on the BBH merger rate (Abbott et al. 2016c, 2018a) generally assume a rate density that is uniform in the comoving volume. As first shown in Fishbach et al. (2018), it is also possible to search for redshift evolution in the rate density using current data. Different redshift-dependent evolutionary behavior is possible (Dominik et al. 2013; Mandel & de Mink 2016; Rodriguez et al. 2016a; Mapelli et al. 2017; Rodriguez & Loeb 2018) with different environments and stellar evolution scenarios (O’Shaughnessy et al. 2010; Belczynski et al. 2016b). For instance, theoretical models of isolated evolution through common envelope lead to a distribution of times to merger $p(t_{GW}) \propto t_{GW}^{-1}$ (Dominik et al. 2012; Belczynski et al. 2016b). This would imply that many isolated binaries will coalesce near their formation redshift and produce a BBH merger rate that approximately tracks the star formation rate, peaking near $z \sim 2$. We find in Section 4 that the current sample of BBH mergers does not provide enough information to confidently constrain any but the most extreme models. While we place more posterior mass on merger rates that increase with increasing redshift than those that decrease, the scenario of a uniform rate in comoving volume is comfortably within our constraints.

BH spin measurements also provide a powerful tool to discriminate between different channels of BBH formation (Mandel & O’Shaughnessy 2010; Abbott et al. 2016d; Rodriguez et al. 2016c; Farr et al. 2017, 2018; Gerosa & Berti 2017; Vitale et al. 2017; Gerosa et al. 2018). For example, BBHs formed in a dynamic environment will have no preferred direction for alignment, producing isotropically oriented spins (Sigurdsson & Hernquist 1993; Portegies Zwart & McMillan 2000; Mandel & O’Shaughnessy 2010; Rodriguez et al. 2015, 2016c; Stone et al. 2017). However, some evidence has been presented for correlation in spin direction due to the natal environment of the progenitor stars within the cluster (Corsaro et al. 2017). In contrast, isolated binaries are expected to preferentially produce mergers with alignment between the spins of the constituent BHs and the orbital angular momentum of the system (Tutukov & Yungelson 1993; Kalogera 2000; Grandclément et al. 2004; Belczynski et al. 2016b; Mandel & de Mink 2016; Marchant et al. 2016; Rodriguez et al. 2016c; O’Shaughnessy et al. 2017; Stevenson et al. 2017b; Gerosa et al. 2018). Other effects occurring in stellar systems like hierarchical triples could also produce a weak preference for certain spin–orbit misalignments (Rodriguez & Antonini 2018). All of our parameterized models point to preferences against high spin magnitudes when the spin tilts are aligned with the orbital angular momentum. In Section 5, we find that the dimensionless spin magnitude inference prefers distributions that decline as the spin magnitude increases from zero, but our ability to distinguish between assumed distributions of spin orientation is very limited.

GW170817, the first binary neutron star merger observed through GW emission (Abbott et al. 2017b), was detected by
GW observatories and associated with a short GRB (Abbott et al. 2017c) in 2017 August. A subsequent post-merger transient (AT 2017gfo) was observed across the electromagnetic spectrum, from radio (Alexander et al. 2017), to near-IR/optical (Chornock et al. 2017; Coulter et al. 2017; Cowperthwaite et al. 2017; Nicholl et al. 2017; Pian et al. 2017; Soares-Santos et al. 2017), to X-ray (Margutti et al. 2017; Troja et al. 2017), to γ-ray (Abbott et al. 2017c; Goldstein et al. 2017; Savchenko et al. 2017). Unfortunately, with only one confident detection, it is not yet possible to infer details of binary neutron star populations more than to note that the GW measurement is mostly compatible with the observed Galactic population (Özel et al. 2012). However, if GW170817 did form a BH, it would also occupy the lower-mass gap described previously.

We structure the paper as follows. First, notation and models are established in Section 2. Section 3 describes our modeling of the BH mass distribution, followed by rate distributions and evolution in Section 4. The BH spin magnitude and orientation distributions are discussed in Section 5. We conclude in Section 6. Studies of various systematics are presented in Appendix A. In Appendix B we present additional studies of spin distributions with model selection for a number of zero-parameter spin models and mixtures of spin orientations. To motivate and enable more detailed studies, we have established a repository of our samples and other derived products.<sup>184</sup>

## 2. Data, Notation, and Models

In this work, we analyze the population of 10 BBH merger events confidently identified in the first and second observing runs (O1 and O2; Abbott et al. 2018a). We do not include marginal detections, but these likely have a minimal impact on our conclusions here (Gaebel et al. 2019). Ordered roughly from smallest to most massive by source-frame chirp mass, the mergers considered in this paper are GW170608, GW151226, GW151012, GW170104, GW170814, GW170809, GW170818, GW150914, GW170823, and GW170729.

The individual properties of those 10 sources were inferred using a Bayesian framework, with results summarized in Abbott et al. (2018a). For BBH systems, two waveform models have been used, both calibrated to numerical relativity simulations and incorporating spin effects, albeit differently: IMRPhenomPv2 (Hannam et al. 2014; Husa et al. 2016; Khan et al. 2016), which includes an effective representation (Schmidt et al. 2015) of precession effects, and SEOBNRv3 (Pan et al. 2014; Taracchini et al. 2014; Babak et al. 2017), which incorporates all spin degrees of freedom. The results presented in this work use IMRPhenomPv2; we discuss potential systematic biases in our inference in Appendix A. We also refer to Appendix B in Abbott et al. (2018a) for more details on comparisons between those two waveform families.

To assess the stability of our results to statistical effects and systematic error, we focus on one modestly exceptional event. Both GW151226 and GW170729 exhibit evidence for measurable BH spin, but GW170729 in particular is an outlier by several other metrics as well. In addition to spins, it is also more massive and more distant than any of the other events in the catalog. All events used in the population analysis have

confident probabilities of astrophysical origin, but GW170729 is the least significant, having the smallest odds ratio of astrophysical versus noise origin (Abbott et al. 2018a). As we describe in Sections 3 and 4, this event has an impact on our inferred merger rate versus both mass and redshift. To demonstrate the robustness of our result, we present these analyses twice: once using every event, and again omitting GW170729.

### 2.1. Binary Parameters

A coalescing compact binary in a quasi-circular orbit can be completely characterized by its eight intrinsic parameters, namely, its component masses $m_i$ and spins $S_i$, and its seven extrinsic parameters: R.A., decl., luminosity distance, coalescence time, and three Euler angles characterizing its orientation (e.g., inclination, orbital phase, and polarization). Binary eccentricity is also a potentially observable quantity in BBH mergers, with several channels having imprints on eccentricity distributions (e.g., Quinlan & Shapiro 1987; Kocsis & Levin 2012; Samsing et al. 2014; Fragione et al. 2018; Rodriguez et al. 2018a). However, our ability to parameterize (Huerta et al. 2014, 2017; Hinder et al. 2018; Klein et al. 2018) and measure (Coughlin et al. 2015; Abbott et al. 2016d, 2017d; Lower et al. 2018) eccentricity is an area of active development. For low to moderate eccentricity at formation, binaries are expected to circularize (Peters 1964; Hinder et al. 2008) before entering the bandwidth of ground-based GW interferometers. We therefore assume zero eccentricity in our models.

In this work, we define the mass ratio as $q = m_2/m_1$, where $m_1 \ge m_2$. The frequency of GW emission is directly related to the component masses. However, due to the expansion of spacetime as the GW is propagating, the frequencies measured by the instrument are redshifted relative to those emitted at the source (Thorne 1983). We capture these effects by distinguishing between masses as they would be measured in the source frame, denoted as above, and the redshifted masses, $(1+z)m_i$, which are measured in the detector frame. Meanwhile, the amplitude of the wave scales inversely with the luminosity distance (Misner et al. 1973). We use the GW measurement of the luminosity distance to obtain the cosmological redshift and therefore convert between detector-frame and source-frame masses. We assume a fixed Planck 2015 (Planck Collaboration et al. 2016) cosmology throughout to convert between a source's luminosity distance and its redshift (Hogg 1999).

We characterize BH spins using the dimensionless spin parameter $\chi_i = S_i/m_i^2$. Of particular interest are the magnitude of the dimensionless spin, $a_i = |\chi_i|$, and the tilt angle with respect to the orbital angular momentum, $\hat{L}$, given by $\cos t_i = \hat{L} \cdot \hat{\chi}_i$. We also define an overall effective spin, $\chi_{\text{eff}}$ (Damour 2001; Racine 2008; Ajith et al. 2011), which is a combination of the individual spin components along the orbital angular momentum:

$$ \chi_{\text{eff}} = \frac{(\chi_1 + q\chi_2) \cdot \hat{L}}{1 + q}. \qquad (1) $$

$\chi_{\text{eff}}$ is approximately proportional to the lowest-order contribution to the GW waveform phase that contains spin for systems with similar masses. Additionally, $\chi_{\text{eff}}$ is conserved throughout

<footnote n="184">184 The data release for this work can be found at https://dcc.ligo.org/LIGO-P1800324/public.</footnote>
<table><caption>**Table 1**
Parameters Describing the Binary Black Hole Population</caption><tr><td>α</td><td>Spectral index of m<sub>1</sub> for the power-law distributed component of the mass spectrum</td></tr><tr><td>m<sub>max</sub></td><td>Maximum mass of the power-law distributed component of the mass spectrum</td></tr><tr><td>m<sub>min</sub></td><td>Minimum BH mass</td></tr><tr><td>β<sub>q</sub></td><td>Spectral index of the mass ratio distribution</td></tr><tr><td>λ<sub>m</sub></td><td>Fraction of binary BHs in the Gaussian component</td></tr><tr><td>μ<sub>m</sub></td><td>Mean mass of BHs in the Gaussian component</td></tr><tr><td>σ<sub>m</sub></td><td>Standard deviation of masses of BHs in the Gaussian component</td></tr><tr><td>δm</td><td>Mass range over which BH mass spectrum turns on</td></tr><tr><td>ζ</td><td>Fraction of binaries with isotropic spin orientations</td></tr><tr><td>σ<sub>i</sub></td><td>Width of the preferentially aligned component of the distribution of BH spin orientations</td></tr><tr><td>E[a]</td><td>Mean of the beta distribution of spin magnitudes</td></tr><tr><td>Var[a]</td><td>Variance of the beta distribution of spin magnitudes</td></tr><tr><td>λ</td><td>How the merger rate evolves with redshift</td></tr></table>

Note. See the text for a more thorough discussion and the functional forms of the models.

the binary evolution to high accuracy (Racine 2008; Gerosa et al. 2015).

## 2.2. Model Features

The current sample is not sufficient to allow for a high-fidelity comparison with models (e.g., population synthesis) that include more detailed descriptions of stellar evolution and environmental influences. As such, we adopt the union of the parameterizations presented in Talbot & Thrane (2017), Fishbach & Holz (2017), Wysocki et al. (2018), Talbot & Thrane (2018), and Fishbach et al. (2018). This allows for better facilitation of comparison between models and the ability to vary the subsets of parameters influencing the mass and spin distributions while leaving others fixed.

The general model family has eight parameters to characterize the mass model, three parameters to characterize each BH's spin distribution, one parameter describing the local merger rate $R_0$, and one parameter characterizing redshift dependence. We refer to the set of these population parameters as $\theta$. All of the population parameters introduced in this section are summarized in Table 1.

## 2.3. Parameterized Mass Models

The power-law distribution considered previously (Abbott et al. 2016c, 2017e) modeled the BBH primary mass distribution as a one-parameter power law, with fixed limits on the minimum and maximum allowed BH mass. With our sample of 10 binaries, we extend this analysis by considering three increasingly complex models for the distribution of BH masses. The first extension, Model A (derived from Fishbach & Holz 2017; Wysocki et al. 2018), allows the maximum BH mass $m_{max}$ and the power-law index $\alpha$ to vary. In Model B (derived from Fishbach & Holz 2017; Kovetz et al. 2017; Talbot & Thrane 2018) the minimum BH mass $m_{min}$ and the mass ratio power-law index $\beta_q$ are also free parameters. However, the priors on Models B and C enforce a minimum of $5 M_\odot$ on $m_{min}$—see Table 2. Explicitly, the mass distribution in

Models A and B takes the form

$$p(m_1, m_2|m_{\text{min}}, m_{\text{max}}, \alpha, \beta_q) \propto \begin{cases} C(m_1)m_1^{-\alpha}q^{\beta_q} & \text{if } m_{\text{min}} \le m_2 \le m_1 \le m_{\text{max}}, \\ 0 & \text{otherwise} \end{cases}, \quad (2)$$

where $C(m_1)$ is chosen so that the marginal distribution is a power law in $m_1$: $p(m_1|m_{\text{min}}, m_{\text{max}}, \alpha, \beta_q) = m_1^{-\alpha}$.

Model A fixes $m_{\text{min}} = 5 M_\odot$ and $\beta_q = 0$, whereas Model B fits for all four parameters. Equation (2) implies that the conditional mass ratio distribution is a power law with $p(q|m_1) \propto q^{\beta_q}$. When $\beta_q = 0$, $C(m_1) \propto 1/(m_1 - m_{\text{min}})$, as assumed in Abbott et al. (2016c, 2017e).

Model C (from Talbot & Thrane 2018) further builds on the mass distribution in Equation (2) by allowing for a second, Gaussian component at high mass, as well as introducing smoothing scales $\delta m$, which taper the hard edges of the low- and high-mass cutoffs of the primary- and secondary-mass power law. The second Gaussian component is designed to capture a possible buildup of high-mass BHs created from pulsational pair-instability supernovae. The tapered low-mass smoothing reflects the fact that parameters such as metallicity probably blur the edge of the lower-mass gap, if it exists. Model C therefore introduces four additional model parameters, the mean, $\mu_m$, and standard deviation, $\sigma_m$, of the Gaussian component; $\lambda_m$, the fraction of primary BHs in this Gaussian component; and $\delta m$, the smoothing scale at the low-mass end of the distribution.

The full form of this distribution is

$$p(m_1|\theta) = \left[ (1 - \lambda_m)A(\theta)m_1^{-\alpha}\Theta(m_{\text{max}} - m_1) + \lambda_m B(\theta)\exp\left(-\frac{(m_1 - \mu_m)^2}{2\sigma_m^2}\right) \right] S(m_1, m_{\text{min}}, \delta m)$$

$$p(q|m_1, \theta) = C(m_1, \theta)q^{\beta_q}S(m_2, m_{\text{min}}, \delta m). \quad (3)$$

The factors A, B, and C ensure that the power-law component, Gaussian component, and mass ratio distributions are correctly normalized. S is a smoothing function that rises from 0 at $m_{\text{min}}$ to 1 at $m_{\text{min}} + \delta m$ as defined in Talbot & Thrane (2018). $\Theta$ is the Heaviside step function. Models A, B, and C are displayed with a selection of parameters for demonstration purposes in the left panels of Figure 1.

## 2.4. Parameterized Spin Models

The BH spin distribution is decomposed into independent models of spin magnitudes, $a$, and orientations, $t$. For simplicity and lacking compelling evidence to the contrary, we assume that both BH spin magnitudes in a binary, $a_i$, are drawn from a beta distribution (Wysocki et al. 2018):

$$p(a_i|\alpha_a, \beta_a) = \frac{a_i^{\alpha_a-1}(1-a_i)^{\beta_a-1}}{B(\alpha_a, \beta_a)}. \quad (4)$$

This distribution is a convenient and flexible parameterization for describing values on the unit interval (Ferrari & Cribari-Neto 2004). Two examples of this distribution are shown in the top right panel of Figure 1. We choose to model the moments of the beta distribution using the mean ($E[a]$) and variance (Var
<table><caption>**Table 2**
Summary of Models Used in Sections 3–5, with the Prior Ranges for the Population Parameters</caption><thead><tr><th rowspan="2">Model</th><th colspan="8">Mass Parameters</th><th colspan="4">Spin Parameters</th></tr><tr><th>α</th><th>m<sub>max</sub></th><th>m<sub>min</sub></th><th>β<sub>q</sub></th><th>λ<sub>m</sub></th><th>μ<sub>m</sub></th><th>σ<sub>m</sub></th><th>δ<sub>m</sub></th><th>E[a]</th><th>Var[a]</th><th>ζ</th><th>σ<sub>i</sub></th></tr></thead><tbody><tr><td>A</td><td>[-4, 12]</td><td>[30, 100]</td><td><b>5</b></td><td><b>0</b></td><td><b>0</b></td><td>N/A</td><td>N/A</td><td>N/A</td><td>[0, 1]</td><td>[0, 0.25]</td><td><b>1</b></td><td>[0, 4]</td></tr><tr><td>B</td><td>[-4, 12]</td><td>[30, 100]</td><td>[5, 10]</td><td>[-4, 12]</td><td><b>0</b></td><td>N/A</td><td>N/A</td><td>N/A</td><td>[0, 1]</td><td>[0, 0.25]</td><td><b>1</b></td><td>[0, 4]</td></tr><tr><td>C</td><td>[-4, 12]</td><td>[30, 100]</td><td>[5, 10]</td><td>[-4, 12]</td><td>[0, 1]</td><td>[20, 50]</td><td>(0, 10)</td><td>[0, 10]</td><td>[0, 1]</td><td>[0, 0.25]</td><td>[0, 1]</td><td>[0, 4]</td></tr></tbody></table>

**Note.** The fixed parameters are in bold. Each of these distributions is uniform over the stated range. All models in this section assume rates that are uniform in the comoving volume ($\lambda = 0$). The lower limit on $m_{min}$ is chosen to be consistent with Abbott et al. (2018a).

[a]), given by

$$ 
\begin{aligned}
\mathbb{E}[a] &= \frac{\alpha_a}{\alpha_a + \beta_a}; \\
\mathrm{Var}[a] &= \frac{\alpha_a \beta_a}{(\alpha_a + \beta_a)^2 (\alpha_a + \beta_a + 1)}.
\end{aligned}
\qquad (5) 
$$

We adopt a prior on the spin magnitude model parameters, which are uniform over the values of $\mathbb{E}[a]$ and $\mathrm{Var}[a]$, which satisfy $\alpha_a, \beta_a \ge 1$. This choice of which values to sample avoids numerically challenging singular spin distributions.

To describe the spin orientation, we assume that the tilt angles between each BH spin and the orbital angular momentum, $t_i$, are drawn from a mixture of two distributions: an isotropic component and a preferentially aligned component, represented by a truncated Gaussian distribution in $\cos t_i$ peaked at $\cos t_i = 1$ (Talbot & Thrane 2017),

$$ p(\cos t_1, \cos t_2 | \sigma_1, \sigma_2, \zeta) = \frac{(1 - \zeta)}{4} + \frac{2\zeta}{\pi} \prod_{i \in \{1, 2\}} \frac{\exp(-(1 - \cos t_i)^2 / (2\sigma_i^2))}{\sigma_i \operatorname{erf}(\sqrt{2}/\sigma_i)}. \qquad (6) $$

We choose to parameterize the cosine of the tilt angles, rather than the angles themselves. This choice prompts the selection of a Gaussian (or uniform) model, rather than a wrapped distribution, which would be more appropriate for an angular variable. An example of the Mixture distribution is displayed in the bottom right panel of Figure 1.

The parameter $\zeta$ denotes the fraction of binaries that are preferentially aligned with the orbital angular momentum; $\zeta = 1$ implies that all BH spins are preferentially aligned, and $\zeta = 0$ is an isotropic distribution of spin orientations. The typical degree of spin misalignment is represented by the $\sigma_i$. For spin orientations we explore two parameterized families of models:

1. Gaussian (G): $\zeta = 1$.

2. Mixture (M): $0 \le \zeta \le 1$.

The Gaussian model is motivated by formation in isolated binary evolution, with significant natal misalignment, while the mixture scenarios allow for an arbitrary combination of this scenario and randomly oriented spins, which arise naturally in dynamical formation.

## 2.5. Redshift Evolution Models

The previous two subsections described the probability distributions of intrinsic parameters $p(\xi)$ (i.e., masses and spins) that characterize the population of BBHs. In addition, we also

measure the value of one extrinsic parameter of the population: the overall merger rate density $R$. The models described in the previous two subsections assume that the distribution of intrinsic parameters is independent of cosmological redshift $z$, at least over the redshift range accessible to the LIGO and Virgo interferometers during the first two observing runs ($z \lesssim 1$). However, we consider an additional model in which the overall event rate evolves with redshift. We follow Fishbach et al. (2018) by parameterizing the evolving merger rate density $R(z)$ in the comoving frame by

$$ R(z|\lambda) = R_0(1 + z)^{\lambda}, \qquad (7) $$

where $R_0$ is the rate density at $z = 0$. In this model, $\lambda = 0$ corresponds to a merger rate density that is uniform in comoving volume and source-frame time, while $\lambda \sim 3$ corresponds to a merger rate that approximately follows the star formation rate in the redshift range relevant to the detections in O1 and O2 (Madau & Dickinson 2014). Various BBH formation channels predict different merger rate histories, ranging from rate densities that will peak in the future ($\lambda < 0$) to rate densities that peak earlier than the star formation rate ($\lambda \gtrsim 3$). These depend on the formation rate history and the distribution of delay times between formation and redshift. In cases where we do not explicitly write the event rate density as $R(z)$, it is assumed that the rate density $R$ is constant in comoving volume and source-frame time.

The general model family, including the distributions of masses, spins, and merger redshift, is therefore given by the distribution

$$ \frac{dN}{d\xi dz}(\theta) = R(z) \left[ \frac{dV_c}{dz}(z) \right] \frac{T_{\text{obs}}}{1+z} p(\xi|\theta), \qquad (8) $$

where $N$ is the total number of mergers that occur within the detection horizon (i.e., the maximum redshift considered) over the total observing time, $T_{\text{obs}}$, as measured in the detector frame, $\theta$ is the collection of all hyperparameters that characterize the distribution, and $dV_c/dz$ is the differential comoving volume per unit redshift. The merger rate density $R(z)$ is related to $N$ by

$$ R(z) = \frac{dN}{dV_c dt}(z), \qquad (9) $$

where $t$ is the time in the source frame, so that Equation (8) can be written equivalently in terms of the merger rate density:

$$ \frac{dR}{d\xi}(z|\theta) = R_0 p(\xi|\theta)(1+z)^{\lambda}. \qquad (10) $$
<figure><img src="image_4.png" /><figcaption>Figure 1. Probability distributions for models encoded by Equations (3), (4), and (6) are shown in the left panels, top right panel, and bottom right panel, respectively. In each, the legend indicates the parameter values for the models plotted. In the case of the lower left (q) distribution, we condition the value of $m_1 = 40 M_\odot$ rather than marginalizing for simplicity.</figcaption></figure>

## 2.6. Hierarchical Population Model

We perform a hierarchical Bayesian analysis, accounting for measurement uncertainty and selection effects (Loredo 2004; Abbott et al. 2016c; Fishbach et al. 2018; Wysocki et al. 2018; Mandel et al. 2019; Mortlock et al. 2018). We model the occurrence rate of events through a Poisson process with a mean dependent on the parameter distribution of the binaries.<fnref n="185" /> The likelihood of the observed GW data given the population hyperparameters $\theta$ that describe the general astrophysical distribution, $dN/d\xi dz$, is given by the inhomogeneous Poisson likelihood:

$$ \mathcal{L}(\{d_n\}|\theta) \propto e^{-\mu(\theta)} \prod_{n=1}^{N_{\text{obs}}} \int \mathcal{L}(d_n|\xi, z) \frac{dN}{d\xi dz}(\theta) d\xi dz, \quad (11) $$

where $\mu(\theta)$ is the rate constant describing the mean number of events as a function of the population hyperparameters, $N_{\text{obs}}$ is the number of detections, and $\mathcal{L}(d_n|\xi, z)$ is the individual-event likelihood for the nth detection having parameters $\xi, z$.

In order to calculate the expected number of detections $\mu(\theta)$, we must understand the selection effects of our detectors. The sensitivity of GW detectors is a strong function of the binary masses and distance and also varies with spin. For any binary, we define the sensitive spacetime volume $VT(\xi)$ of a network with a given sensitivity to be

$$ VT(\xi) = T_{\text{obs}} \int_{0}^{\infty} f(z|\xi) \frac{dV_c}{dz} \frac{1}{1+z} dz, \quad (12) $$

where the sensitivity is assumed to be constant over the observing time, $T_{\text{obs}}$, as measured in the detector frame and

$f(z|\xi)$ is the detection probability of a BBH with the given parameter set $\xi$ at redshift $z$ (O'Shaughnessy et al. 2010), averaged over the extrinsic binary orientation parameters (Finn & Chernoff 1993). The factor of $1/(1+z)$ arises from the difference in clocks timed between the source frame and the detector frame. For a given population with hyperparameters $\theta$, we can calculate the total observed spacetime volume

$$ \langle VT \rangle_{\theta} = \int_{\xi} p(\xi|\theta)VT(\xi)d\xi, \quad (13) $$

where $p(\xi|\theta)$ describes the underlying distribution of the intrinsic parameters. We performed large-scale simulation runs wherein the spacetime volume in the above equation is estimated by Monte Carlo integration (Tiwari 2018)—these runs are restricted to have no BH less massive than $5 M_\odot$. We then use a semianalytic prescription, calibrated to the simulation results, to derive the $\langle VT \rangle_\theta$ for specific hyperparameters.

Allowing the merger rate to evolve with redshift, the expected number of detections is given by

$$ \mu(\theta) = T_{\text{obs}} \int_{\xi} \int_{0}^{\infty} p(\xi|\theta)f(z|\xi)R(z) \frac{dV_c}{dz} \frac{1}{1+z} dz d\xi. \quad (14) $$

If the merger rate does not evolve with redshift, i.e., $R(z) = R_0$, this reduces to $\mu(\theta) = R_0 \langle VT \rangle_\theta$.

We note that the hyperparameter likelihood given by Equation (11) reduces to the likelihood used in the O1 mass distribution analysis (Equation (D10) of Abbott et al. 2016c), which fit only for the shape, not the rate/normalization of the mass distribution, if one marginalizes over the rate parameter with a flat-in-log prior $p(R_0) \propto 1/R_0$ (Fishbach et al. 2018; Mandel et al. 2019). For consistency with previous analyses, we adopt a flat-in-log prior on the rate parameter throughout this work.

<footnote n="185">¹⁸⁵ While this assumption is embedded (Farr et al. 2015b) in the selection of events used in this work, studies of event count per time do not show significant evidence for deviations from Poissonian statistics (Abbott et al. 2018a).</footnote>
## 2.7. Statistical Framework and Prior Choices

In practice, we sample the likelihood $L(d_n|\xi, z)$ using the parameter estimation pipeline LALInference (Veitch et al. 2015). Since LALInference gives us a set of posterior samples for each event, we first divide out the priors used in the individual-event analyses before applying Equation (11) (Hogg et al. 2010; Mandel 2010; see Appendix C).

Where not fixed, we adopt uniform priors on population parameters describing the models. Unless otherwise noted, for the event rate distribution we use a log-uniform distribution in $R_0$, bounded between $[10^{-1}, 10^3]$. While this is a different form than the priors adopted in Abbott et al. (2018a), we note that similar results are obtained on the rates (see Section 4), indicating that the choice of prior does not strongly influence the posterior distributions. We provide specific limits on all priors when the priors for a given model are introduced. Unless otherwise stated, all posterior credible intervals are 90% intervals, symmetric in the quantiles around the median. The MCMC-based analyses presented in this work have approximately 10⁴ effective samples, after thinning by their autocorrelation time.

The normalization factor of the posterior density in Bayes's theorem is the evidence—it is the probability of the data given the model. We are interested in the preferences of the data for one model versus another. This preference is encoded in the Bayes factor, or the ratio of evidences. The odds ratio is the Bayes factor multiplied by their ratio of the model prior probabilities. In all cases presented here, the prior model probabilities are assumed to be equal, and odds ratios are equivalent to Bayes factors.

We often present the posterior population distribution (PPD) of various quantities. The PPD is the expected distribution of new mergers conditioned on previously obtained observations. It integrates the distribution of values (e.g., $\xi$, such as the masses and spins) conditioned on the model parameters (e.g., the power-law index) over the posteriors obtained for the model parameters:

$$p(\xi_{\text{new}}|\xi_{\text{observed}}) = \int p(\xi_{\text{new}}|\theta)p(\theta|\xi_{\text{observed}})d\theta. \quad (15)$$

It is a predictor for future merger values $\xi_{new}$ given observed data $\xi_{observed}$ and factors in the uncertainties imposed by the posterior on the model parameters. Note that the PPD does not incorporate the detector sensitivity and therefore is not a straightforward predictor of the properties of future observed mergers.

## 3. The Mass Distribution

For context, Figure 4 in Abbott et al. (2018a) illustrates the inferred masses for all of the significant BBH observations identified in our GW surveys in O1 and O2. Despite at least moderate sensitivity to total masses between 0.1 and 500 $M_\odot$, current observations occupy only a portion of the binary mass parameter space. Notably, we have not yet observed a pair of very massive (e.g., 100 $M_\odot$) BHs, a binary that is bounded away from equal mass in its posterior, or a binary with a component mass confidently below 5 $M_\odot$. In our survey, we also find a preponderance of observations at higher masses: six with significant posterior support above 30 $M_\odot$. In this section, we attempt to reconstruct the BBH merger rate as a function of the component masses using parameterized models. Table 2

summarizes the mass models adopted from Section 2.3 and the prior distributions for each of the parameters in those models. We present results for three increasingly general mass and spin models, the most complex of which ranges over the full set of model parameters in Section 2, with the exception of dependence of rate on redshift. The interdependence of the mass and redshift distribution is explored more fully in Section 4.

### 3.1. Parameterized Modeling Results

Figure 2 shows our updated inference for the compact binary primary mass $m_1$ and mass ratio $q$ distributions for several increasingly general population models. In addition to inferring the mass distribution, all of these calculations self-consistently marginalize over the parameterized spin distribution presented in Section 5 and the merger rate. Figures 3 and 4 show the posterior distribution on selected model hyperparameters.

If we assume that the BH masses follow a power-law distribution and fix the minimum BH mass to be $m_{min} = 5 M_\odot$ (Model A), we find $\alpha = 0.4^{+1.4}_{-1.9}$ and $m_{max} = 41.6^{+9.6}_{-4.3} M_\odot$. In Model B we infer the power-law index of the primary mass to be $\alpha = 1.3^{+1.4}_{-1.7}$ with corresponding limits $m_{min} = 7.8^{+1.2}_{-2.5} M_\odot$ and $m_{max} = 40.8^{+11.8}_{-4.4} M_\odot$.

Figure 4, shows the posterior over the population parameters present in Models A and B, as well as a second, Gaussian population parameterized with $m_{max}$ and $\sigma_m$. $\lambda_m$ is the mixing fraction of binaries in the Gaussian population versus the power law, with $\lambda_m = 0$ indicating only the power-law component. The Gaussian component is centered at $\mu_m = 29.8^{+5.8}_{-7.3} M_\odot$, has a width $\sigma_m = 6.4^{+3.2}_{-4.2} M_\odot$, and is consistent with the parameters of the seven highest-mass events in our sample as seen in Figure 5. Also as a consequence of this mixture, the second component can account for many of the high-mass events. Without needing to accommodate higher-mass events, the inferred power law is much steeper ($\alpha = 7.1^{+4.4}_{-4.8}$) than Model A or Model B; however, the posterior distribution for Model C is less informative for $\alpha \gtrsim 4$. This in turn means that we cannot constrain the parameter $m_{max}$ in Model C since the power-law component has negligible support above ~45 $M_\odot$ (see the top panel of Figure 2). In the intermediate regime, ~15–25 $M_\odot$, Model C infers a smaller rate than Model A or Model B as a consequence of the steeper power-law behavior. The low-mass smoothing allowed in this model also weakens constraints we can place on the minimum BH mass; in this model we find $m_{min} = 6.9^{+1.7}_{-2.8} M_\odot$. All three models produce consistent results for the marginal merger rate distribution, as is further discussed in Section 4.

All models feature a parameter $m_{max}$, which defines a cutoff of the power law. However, the interpretation of that parameter within Model C is not a straightforward comparison with Models A and B, due to the presence of the Gaussian component at high mass and the large value of the power-law spectral index. Instead, to compare those two features, we compute the 99th percentile of the mass distribution inferred from the model PPDs (see Equation (15)). Model A obtains 44.0 $M_\odot$, Model B obtains 41.8 $M_\odot$, and Model C obtains 41.8 $M_\odot$. Therefore, all models self-consistently infer a dearth of BHs above ~45 $M_\odot$. This is determined by the lower limit for the mass of the most massive BH in the sample because $m_{max}$ can be no smaller than this value. Similarly, the models that allow $m_{min}$ to vary (B and C) disfavor populations with
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Inferred differential merger rate as a function of primary mass, $m_1$, and mass ratio, $q$, for three different assumptions. For each of the three increasingly complex assumptions A, B, C described in the text we show the PPD (dashed) and median (solid), as well as the 90% symmetric credible intervals (shaded regions), for the differential rate. The results shown marginalize over the spin distribution model. The falloff at small masses in Models B and C is driven by our choice of the prior limits on the $m_{\min}$ parameter (see Table 2). All three models give consistent mass distributions within their 90% credible intervals over a broad range of masses, consistent with their near-unity evidence ratios (Table 3); in particular, the peaks and trough seen in Model C, while suggestive, are not identified at high credibility in the mass distribution.</figcaption></figure>

$m_{\min}$ above $\simeq 9 M_{\odot}$. This parameter is close to the largest allowed mass for the least massive BH in the sample, for similar reasons.

The lower limits we place on $m_{\min}$ are dominated by our prior choices that constrain $m_{\min} \in [5, 10] M_{\odot}$ (see Table 2). For example, in Figure 3, the posterior on $m_{\min}$ becomes flat as $m_{\min}$ approaches the prior boundary at $5 M_{\odot}$. Given current sensitivities, this is to be expected (Littenberg et al. 2015; Mandel et al. 2015). In the inspiral-dominated regime, the sensitive time volume scales as $VT \sim m^{15/6}$ (Finn & Chernoff 1993); extending our inferred mass distributions and merger rates into the possible lower BH mass gap from 3–5 $M_{\odot}$ (Özel et al. 2010; Farr et al. 2011b; Kreidberg et al. 2012) yields an expected number of detected BBH mergers $\lesssim 1$. Thus, we are unable to place meaningful constraints on the presence or absence of a mass gap at low BH mass.

Models B and C also allow the distribution of mass ratios to vary according to $\beta_q$. In these cases the inferred mass ratio distribution favors comparable-mass binaries (i.e., distributions with most support near $q \simeq 1$); see the bottom panel of Figure 2. Within the context of our parameterization, we find $\beta_q = 6.9_{-5.7}^{+4.6}$ for Model B and $\beta_q = 4.5_{-5.2}^{+6.6}$ for Model C. These values are consistent with each other and are bounded above zero at 95% confidence, thus implying that the mass ratio distribution is nearly flat or declining with more extreme mass ratios. The posterior on $\beta_q$ returns the prior for $\beta_q \gtrsim 4$. Thus, we cannot say much about the relative likelihood of asymmetric binaries, beyond their overall rarity.

The distribution of the parameter controlling the fraction of the power law versus the Gaussian component in Model C is $\lambda_m = 0.3_{-0.2}^{+0.4}$, which peaks away from zero, implying that this

model prefers a contribution to the mass distribution from the Gaussian population in addition to the power laws modeled in A and B. To determine preference among the three models presented in this section, we compute the Bayes factors comparing the mass models using a nested sampler (Skill- ing 2004), cPNest (Veitch et al. 2017). These are shown in Table 3. Model B, which allows $m_{\min}$ and $\beta_q$ to vary, is preferred over Model A (ln $BF_B^A = -1.42$). To isolate the contributions of the Gaussian component and low-mass smoothing in Model C, we compute the Savage-Dickey density ratio, $p(\theta = 0)/p_{\text{prior}}(\theta = 0)$, equivalent to the Bayes factor comparing without and with the feature. The model including a Gaussian component in addition to the power-law distribution is preferred over the pure power-law models (ln $BF_C^{\lambda=0} = -1.92$); nevertheless, all models infer mass distribu- tions that agree within their 90% credible bounds (see Figure 2). We caution that the mild preferences in Table 3 are influenced by our choices of the range and shape of the priors we apply to the parameters, particularly for models where the number of parameters is comparable to the number of events. Moreover, the credible intervals on the distributions of the primary mass overlap, indicating that the model predictions agree to within the individual model uncertainties. We are unable to distinguish between a gradual and a sharp cutoff at low mass (ln $BF_C^{\delta_m=0} = 0.14$). This is unsurprising, since we are less sensitive to structure in the mass distribution at low masses (Talbot & Thrane 2018).

The analysis above includes all 10 BBH detections, though not all events have the same statistical detection confidence (Gaebel et al. 2019). To assess the stability of our results against systematics in the estimated significance, we have
<figure><img src="image_6.png" /><figcaption>**Figure 3.** 1D and 2D posterior distributions for the hyperparameters describing Models A and B. Large values of $\alpha$ correspond to a mass distribution that rapidly decays with increasing mass. Large values of $\beta$ correspond to a mass ratio distribution that prefers equal-mass binaries. Also shown is the 1D posterior distribution for the merger rate discussed in Abbott et al. (2018a) and the stability of Model A to the removal of the GW170729 event.</figcaption></figure>

repeated these analyses after omitting the least significant detection. For our sample, the least significant detection, GW170729, is also the most massive binary. Most features we derive from our observations remain unchanged, with one exception shown in Figure 3: since we have omitted the most massive binary, the maximum BH mass $m_{max}$ reported in Models A and B is decreased by about 5 $M_{\odot}$. Without GW170729, the $m_{max}$ distribution is $38.3_{-3.6}^{+7.3} M_{\odot}$ for Model A and $37.3_{-4.5}^{+8.5} M_{\odot}$ for Model B. This is consistent with the difference between GW170729 and the next-highest-mass binary, GW170823, when comparing the less massive end of their primary mass posteriors.

## 3.2. Comparison with Theoretical and Observational Models

Previous modeling of the primary mass distribution with a power-law distribution (Abbott et al. 2016c) was last updated with the discovery of GW170104 (Abbott et al. 2017e). This analysis measured spectral index of the power law to be $\alpha = 2.3_{-1.4}^{+1.3}$ at 90% confidence assuming a minimum BH mass of 5 $M_{\odot}$ and maximum total mass of 100 $M_{\odot}$. None of our models directly emulate this one, but Model A is the closest analog. When allowing $m_{max}$ to vary, 100 $M_{\odot}$ is strongly disfavored. As a consequence of the lower $m_{max}$, the power-law index inferred is also shallower than previously obtained (Fishbach & Holz 2017) but remains consistent with the previous distribution.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** 1D and 2D posterior distributions for the hyperparameters describing Model C. This model consists of the power-law distribution in Model B with an additional Gaussian component at high mass. The parameters $\alpha$, $\beta$, $m_{\text{max}}$, and $m_{\text{min}}$ describe the power-law component. The Gaussian has mean $\mu_m$ and standard deviation $\sigma_m$. The fraction of BHs in the Gaussian component is $\lambda_m$. This model also allows for a gradual turn-on at low masses over a mass range $\delta_m$.</figcaption></figure>

In Figure 5, we highlight the two mass gaps predicted by models of stellar evolution: the first gap between ~2 and ~5 $M_{\odot}$ and the second between ~50 and ~150 $M_{\odot}$, compared against the observed BHs. A set of tracks (Spera & Mapelli 2017) relating the progenitor mass and compact object is also shown for reference purposes. The tracks are subject to many uncertainties in stellar and binary evolution and only serve as representative examples. We discuss some of those uncertainties in the context of our results below.

The minimum mass of a BH and the existence of a mass gap between neutron stars and BHs (lower gray shaded area,

right panel of Figure 5) are currently debated. Claims (Özel et al. 2010; Farr et al. 2011b) of the existence of a mass gap between the heaviest neutron stars (~2 $M_{\odot}$) and the lightest BHs (~5 $M_{\odot}$) are based on the sample of about a dozen XRBs with dynamical mass measurements. However, Kreidberg et al. (2012) suggested that the dearth of observed BH masses in the gap could be due to a systematic offset in mass measurements. Moreover, only a subset of theoretical models (e.g., the “rapid” model in Fryer et al. 2012) reproduce this gap in stellar modeling. We can see in Figure 5 that none of the observed binaries sit in this gap, but the sample is not
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left panel: compact-object masses ($m_{CO}$) from GW detections in O1 and O2, with the black squares and error bars representing the component masses of the merging BHs and their uncertainties, and red triangles representing the mass and associated uncertainties of the merger products. The horizontal green line shows the 99th percentile of the mass distribution inferred from the Model B PPD. Right panel: predicted compact-object mass as a function of the zero-age main-sequence mass of the progenitor star ($m_{ZAMS}$) and for four different metallicities of the progenitor star (ranging from $Z = 10^{-4}$ to $Z = 2 \times 10^{-2}$; Spera & Mapelli 2017). This model accounts for single stellar evolution from the PARSEC stellar evolution code (Bressan et al. 2012), for core-collapse supernovae (Fryer et al. 2012), and for pulsational pair-instability and pair-instability supernovae (Woosley 2017). The shaded areas represent the lower- and upper-mass gaps. There is uncertainty as to the final product of GW170817. It is shown in the left panel to emphasize that BNS mergers might fill the lower gap.</figcaption></figure>

<table><caption>**Table 3**
The Log Bayes Factor Comparing Each of the Models Described in Table 2 to the Most Complex Model, Model C</caption><thead><tr><th>Model</th><th>A</th><th>B</th><th>C, λ<sub>m</sub> = 0</th><th>C, δ<sub>m</sub> = 0</th></tr></thead><tbody><tr><td>ln BF<sub>C</sub><sup>i</sup></td><td>-2.28</td><td>-0.86</td><td>-1.92</td><td>0.14</td></tr></tbody></table>

**Note.** The evidence for the three mass models is computed using nested sampling, while the limits $\lambda_m = 0$ and $\delta_m = 0$ of Model C are computed using the Savage-Dickey density ratio.

sufficient to definitively confirm or refute the existence of this mass gap.

From the first six announced BBH detections, Fishbach & Holz (2017) argued that there is evidence for missing BHs with mass greater than $\ge 40 M_{\odot}$. The existence of this second mass gap—see the upper gray shaded area in the right panel of Figure 5 between ~50 and ~150 $M_{\odot}$—has been further explored by Talbot & Thrane (2018), Wysocki et al. (2018), Bai et al. (2018), and Roulet & Zaldarriaga (2019). This gap might arise from the combined effect of pulsational pair-instability (Barkat et al. 1967; Heger et al. 2003; Woosley et al. 2007; Woosley 2017) and pair-instability (Fowler & Hoyle 1964; Ober et al. 1983; Bond et al. 1984) supernovae. Uncertainties in stellar evolution models (e.g., stellar winds, rotation) and in the treatment of the final outcomes of (pulsational) pair instability lead to a range of possible low-mass edges for the upper-mass gap, as well as the shape and abundance in a putative buildup. Predictions for the maximum

mass of BHs born after pulsational pair-instability supernovae are ~50 $M_{\odot}$ (Belczynski et al. 2016a; Spera & Mapelli 2017). Our inferred maximum mass is consistent with these predictions.

## 4. Merger Rates and Evolution with Redshift

As illustrated in previous work (Abbott et al. 2016e, 2018a; Fishbach & Holz 2017; Fishbach et al. 2018; Wysocki et al. 2018), the inferred BBH merger rate depends on and correlates with our assumptions about their intrinsic mass (and to a lesser extent, spin) distribution. In the most recent catalog of GW BBH events (Abbott et al. 2018a), we infer the overall BBH merger rate for two fixed-parameter populations. The first of these populations follows the power-law model given by Equation (2) with $\alpha = 2.3$, $\beta_q = 0$, $m_{min} = 5 M_{\odot}$, and $m_{max} = 50 M_{\odot}$. The second population follows a distribution in which both BH masses are independently drawn from a flat-in-log distribution:

$$p(m_1, m_2) \propto \frac{1}{m_1 m_2}, \qquad (16)$$

subject to the same mass cutoffs $5 M_{\odot} < m_2 < m_1 < 50 M_{\odot}$ as the fixed power-law population. Both the power-law and flat-in-log populations assume an isotropic and uniform-magnitude spin distribution ($\alpha_a = \beta_a = 1$). These two fixed-parameter populations are used to estimate the population-averaged sensitive volume $\langle VT \rangle$ with a Monte Carlo injection campaign as described in Abbott et al. (2018a), with each population
<table><caption>**Table 4**
Summary of Models in Section 4, with Prior Ranges for the Population Parameters Determining the Rate Models</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Mass Model</th><th rowspan="2">Rate Parameters<br>λ</th><th colspan="4">Spin Parameters</th></tr><tr><th>α<sub>a</sub></th><th>β<sub>a</sub></th><th>E[a]</th><th>Var[a]</th></tr></thead><tbody><tr><td>Fixed parameter (power law)</td><td>A, with α = 2.3,<br>m<sub>max</sub> = 50 M<sub>☉</sub></td><td><b>0</b></td><td><b>1</b></td><td><b>1</b></td><td>N/A</td><td>N/A</td></tr><tr><td>Fixed parameter (flat-in-log)</td><td>Equation (16)</td><td><b>0</b></td><td><b>1</b></td><td><b>1</b></td><td>N/A</td><td>N/A</td></tr><tr><td>Nonevolving</td><td>A, B, C</td><td><b>0</b></td><td>N/A</td><td>N/A</td><td>[0, 1]</td><td>[0, 0.25]</td></tr><tr><td>Evolving<sup>a</sup></td><td>A</td><td>[-25, 25]</td><td>N/A</td><td>N/A</td><td><b>0</b></td><td><b>0</b></td></tr></tbody></table>

**Notes.** The fixed-parameter models are drawn from Abbott et al. (2018a). The fixed parameters are in bold. Each of these distributions is uniform over the stated range; as previously, we require $α_a, β_a ≥ 1$. Details of the mass models listed here are described in Table 2.

<sup>a</sup> This model assumes that the BHs have zero spin.

corresponding to a different $\langle VT \rangle$ because of the strong correlation between the mass spectrum and the sensitive volume. Under the assumption of a constant-in-redshift rate density, these $\langle VT \rangle$ estimates yield two different estimates of the rate: $57_{-25}^{+40}$ Gpc$^{-3}$ yr$^{-1}$ for the $α = 2.3$ population and $19_{-8.2}^{+13}$ Gpc$^{-3}$ yr$^{-1}$ for the flat-in-log population (90% credibility; combining the rate posteriors from the two analysis pipelines).

The two fixed-parameter distributions do not incorporate all information about the mass, mass ratio, spin distribution, and redshift evolution suggested by our observations in O1 and O2. In this section, rather than fixing the mass and spin distribution, we estimate the rate by marginalizing over the uncertainty in the underlying population, which we parameterize with the mass and spin models employed in Sections 3 and 5. When carrying out these analyses, it is computationally infeasible to determine $VT(\xi)$ for each point in parameter space with the full Monte Carlo injection campaign described in Abbott et al. (2018a), so we employ the semianalytic methods described in Appendix A. Furthermore, while the rate calculations in Abbott et al. (2018a) incorporate all triggers down to a very low threshold and fit the number of detections by modeling the signal and background distributions in the detection pipelines (Farr et al. 2015b; Abbott et al. 2016e), in this work we fix a high detection threshold Abbott et al. (2018a), which sets the number of detections to $N_{obs} = 10$. In principle, our results are sensitive to the choice of threshold, but this effect has been shown to be much smaller than the statistical uncertainties (Gaebel et al. 2019). The choice of detection threshold is further discussed in Appendix A. The full set of models used in this section is enumerated in Table 4.

In these calculations, we first maintain the assumption in Abbott et al. (2018a) that the merger rate is uniform in comoving volume and source-frame time, as discussed in Section 2. We then relax this assumption and consider a merger rate that evolves in redshift according to Equation (7), fitting the mass distribution jointly with the rate density as a function of redshift.

### 4.1. Nonevolving Merger Rate

We first consider the case of a uniform-in-volume merger rate and examine the effects of fitting the rate jointly with the distribution of masses and spins. The first column in Figures 3 and 4 shows the results of self-consistently determining the rate using the models for the mass and spin distribution described in the previous two sections.

<table><caption>**Table 5**
BBH Merger Rate Intervals for Each of the Mass Models Tested</caption><thead><tr><th>Model</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>R<sub>0</sub> (Gpc<sup>-3</sup> yr<sup>-1</sup>)</td><td>64.0<sup>+73.5</sup><sub>-33.0</sub></td><td>53.2<sup>+55.8</sup><sub>-28.2</sub></td><td>58.3<sup>+72.3</sup><sub>-32.2</sub></td></tr></tbody></table>

**Note.** These rates assume no evolution in redshift but otherwise marginalize over all other population parameters.

Table 5 contains the intervals on the distribution of $R_0$ for all three models. For Models B and C we deduce a merger rate in the range of $R_0 = 24.9–109.0$ Gpc$^{-3}$ yr$^{-1}$. Adopting Model A for the mass distribution yields a slightly higher rate estimate, $R_0 = 31.0–137.5$ Gpc$^{-3}$ yr$^{-1}$, as this model fixes $m_{min} = 5 M_{☉}$, whereas Models B and C favor a higher minimum mass and therefore larger population-averaged sensitive volumes. The rate estimates are consistent between all mass models considered, including the results presented for the fixed-parameter power-law model in Abbott et al. (2018a). However, the fixed-parameter models in Abbott et al. (2018a) are disfavored by our full fit to the mass distribution, particularly with respect to the maximum mass. Our results favor maximum masses $\lesssim 45 M_{☉}$, rather than $50 M_{☉}$ as used in Abbott et al. (2018a), and power-law slopes closer to $α \sim 1$. For this reason, although we infer a mass distribution slope that is similar to the flat-in-log population from Abbott et al. (2018a), we infer a rate that is closer to the rate inferred for the fixed-parameter power-law model.186 While $\langle VT \rangle$ gets larger (implying a smaller rate estimate) as $α$ is decreased, decreasing $m_{max}$ has the opposite effect, and so the $\langle VT \rangle$ for the fixed-parameter power-law model is similar to the $\langle VT \rangle$ values for our best-fit mass distributions, which favor smaller $α$ and smaller $m_{max}$.

We note that while our analysis differs from the rate calculations in Abbott et al. (2018a) by the choice of prior on the rate parameter (log-uniform in this work compared to a Jeffreys prior $p(R_0) \propto R_0^{-0.5}$ in Abbott et al. 2018a), adopting a Jeffreys prior has a negligible effect on our rate posteriors. For example, under a log-uniform prior, we recover a rate for Model A of $62.8_{-33.3}^{+74.0}$ Gpc$^{-3}$ yr$^{-1}$, whereas under a Jeffreys prior this shifts by only $\sim 10\%$ to $56.7_{-30.4}^{+65.4}$ Gpc$^{-3}$ yr$^{-1}$.

<footnote n="186">186 The flat-in-log population (Equation (16)) cannot be parameterized by the mass Models A, B, and C used in this work because the mass ratio distribution takes a different form. However, it is very close to Model A with $α = 1$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Constraints on evolution of the BBH merger rate density as a function of redshift. Including the 10 BBHs from O1 and O2 in our analysis, we find a preference for a merger rate that increases with increasing redshift. The solid blue line gives the posterior median merger rate density, and dark and light bands give 50% and 90% credible intervals, respectively. In orange, the solid line and shaded region show the median and 90% credible interval, respectively, of the rate inferred for Model B as discussed in Section 4.1, assuming a nonevolving merger rate.</figcaption></figure>

## 4.2. Evolution of the Merger Rate with Redshift

As discussed in the introduction, most formation channels predict some evolution of the merger rate with redshift, due to factors such as the star formation rate, time-delay distribution, metallicity evolution, and globular cluster formation rate (-Dominik et al. 2013; Belczynski et al. 2016b; Mandel & de Mink 2016; Rodriguez & Loeb 2018). Therefore, in this section we allow the merger rate to evolve with redshift, and we infer the redshift evolution jointly with the mass distribution. For simplicity, we adopt the two-parameter Model A for the mass distribution and fix spins to zero for this analysis. As discussed in Section 3, the additional mass and spin degrees of freedom have only a weak effect on the inferred merger rate. We assume the redshift evolution model given by Equation (7). Because massive binaries are detectable at higher redshifts, the observed redshift evolution correlates with the observed mass distribution of the population, and so we must fit them simultaneously. However, as in Fishbach et al. (2018), we assume that the underlying mass distribution does not vary with redshift. We therefore fit the joint mass-redshift distribution according to the model:

$$ \frac{dR}{dm_1 dm_2}(z) = R_0 p(m_1, m_2 | \alpha, m_{\text{max}})(1+z)^{\lambda}. \quad (17) $$

Note that this model assumes that the merger rate density increases or decreases monotonically with redshift over the sensitive range $z < 1$. If the merger rate follows the star formation rate, we expect the rate to peak around $z \sim 2$, which is currently far beyond the horizon redshift for BBH detections.

Figure 6 shows the merger rate density as a function of redshift (blue band), compared to the rate inferred in Section 4.1 for the nonevolving Model B (orange band). The joint posterior probability density function (pdf) on $\lambda, \alpha$, and $m_{\text{max}}$, marginalized over the local rate parameter $R_0$, is shown in Figure 7. There is a strong correlation between the mass power-law slope and the redshift evolution parameter. This is due to the fact that higher-mass BBHs are detectable at higher redshifts, and so, for the same underlying mass distribution, an increasing rate evolution with redshift implies that a greater

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Posterior pdf on the redshift evolution parameter $\lambda$, mass power-law slope $\alpha$, and maximum mass $m_{\text{max}}$, marginalized over the local rate parameter $R_0$, and assuming a flat prior on $\lambda, \alpha$, and $m_{\text{max}}$ and a flat-in-log prior on $R_0$. In order to analyze the stability of the model against outliers, we repeat the analysis once with the sample of 10 BBHs (results shown in blue) and once excluding the most distant and massive event in our sample, GW170729 (results shown in red). The contours show 50% and 90% credible intervals. The dashed black lines show the values of hyperparameters assumed for the fixed-parameter power-law model. We infer a redshift evolution that is consistent with a flat in comoving volume and source-frame time merger rate ($\lambda = 0$) with a preference for $\lambda \ge 0$ at 0.93 credibility when considering all 10 events. This preference becomes less significant with the exclusion of GW170729 from the analysis. The inferred power-law slope and maximum mass are consistent with the values inferred in Section 3. This analysis recovers a broader posterior on the mass power-law slope because of the correlation with the redshift evolution parameter, but the maximum mass remains well constrained at $\lesssim 45 M_{\odot}$.</figcaption></figure>

fraction of detected BBHs will be massive. This effect is hard to disentangle from a shallower mass distribution, which will also produce comparatively more massive BBH detections. Note that the constraints on $\alpha$ and $m_{\text{max}}$ in Section 3 are correlated by the same effect. Compared to the constraints on $\alpha$ and $m_{\text{max}}$ discussed in Section 3, which assume a constant-in-redshift merger rate density, allowing for additional freedom in the redshift distribution of BBHs relaxes the constraints on the mass distribution parameters, especially the power-law slope $\alpha$ ($m_{\text{max}}$ is sufficiently well measured that the correlation with $\lambda$ is not as noticeable). Under the assumption of a constant merger rate density, Model A in Section 3 finds $\alpha = 0.4_{-1.9}^{+1.4}$, $m_{\text{max}} = 41.6_{-4.3}^{+9.6} M_{\odot}$, whereas allowing for redshift evolution yields $\alpha = 1.8_{-2.0}^{+1.7}$, $m_{\text{max}} = 41_{-5}^{+11} M_{\odot}$ when analyzing the sample of 10 BBHs from O1 and O2. As in Section 3, we carry out a leave-one-out analysis, excluding the most massive and distant BBH, GW170729, from the sample (red curves in Figure 7). Without GW170729, the marginalized mass distribution posteriors become $\alpha = 0.9_{-2.2}^{+1.8}$, $m_{\text{max}} = 38_{-4}^{+10} M_{\odot}$.

Marginalizing over the two mass distribution parameters and the redshift evolution parameter, the merger rate density is consistent with being constant in redshift ($\lambda = 0$), and in particular, it is consistent with the rate estimates recovered under the different mass distribution models in Section 4.1
above. However, we find a preference for a merger rate density that increases at higher redshift ($\lambda \ge 0$) with probability 0.93. This implies that models that predict a constant or slightly decreasing merger rate with redshift, such as certain models of primordial BHs (Mandic et al. 2016), are disfavored. This preference for a merger rate that increases with increasing redshift becomes less significant when GW170729 is excluded from the analysis because this event likely merged at redshift $z \gtrsim 0.5$, close to the O1–O2 detection horizon. Although GW170729 shifts the posterior toward larger values of $\lambda$, implying a stronger redshift evolution of the merger rate, the posterior remains well within the uncertainties inferred from the remaining nine BBHs. When including GW170729 in the analysis, we find $\lambda = 8.4_{-9.5}^{+9.6}$ at 90% credibility, compared to $\lambda = 2.3_{-10.9}^{+9.9}$ when excluding GW170729 from the analysis. With only 10 BBH detections so far, the wide range of possible values for $\lambda$ is consistent with most astrophysical formation channels. The precision of this measurement will improve as we accumulate more detections in future observing runs and may enable us to discriminate between different formation rate histories or time-delay distributions (Sathyaprakash et al. 2012; Van Den Broeck 2014; Fishbach et al. 2018).

## 5. The Spin Distribution

The GW signal depends on spins in a complicated way, but at leading order, and in the regime we are interested in here, some combinations of parameters have more impact on our inferences than others and thus are measurable. One such parameter is $\chi_{\text{eff}}$. For binaries that are near equal mass, we can see from Equation (1) that only when BH spins are high and aligned with the orbital angular momentum will $\chi_{\text{eff}}$ be measurably greater than zero. Figure 5 in Abbott et al. (2018a) illustrates the inferred $\chi_{\text{eff}}$ spin distributions for all of the BBHs identified in our GW surveys in O1 and O2. Only GW170729 and GW151226 show significant evidence for positive $\chi_{\text{eff}}$; the rest of the posteriors cluster around $\chi_{\text{eff}} = 0$.

Despite these degeneracies, several tests have been proposed to use spins to constrain BBH formation channels (Farr et al. 2017, 2018; Gerosa & Berti 2017; Stevenson et al. 2017a; Talbot & Thrane 2017; Vitale et al. 2017; Gerosa et al. 2018; Wysocki et al. 2018). Drawing on these methods, we now seek to estimate the BH spin magnitude and misalignment distributions, under different assumptions regarding isotropy or alignment.

### 5.1. Spin Magnitude and Tilt Distributions

We examine here the individual spin magnitudes and tilt distributions. Throughout this section, when referring to the parametric models, we also allow the merger rate and population parameters describing the most general mass model to vary (Model C, see Table 2). Changing the parameterization of the mass model does not significantly change our inferences about the spin distribution. However, to account for degeneracies between mass and spin that grow increasingly significant for longer, low-mass signals (Baird et al. 2013), we must consistently model the mass and spin distributions together. See Table 6 for a summary of the models and priors used in this section.

The inferred distributions of spin magnitude are shown in Figure 8. The top panel shows the PPD and the median and associated uncertainties on the spin magnitude inferred from

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Inferred distribution of spin magnitude for a parametric (top) and nonparametric binned model (bottom). Both component magnitudes are included in these distributions. The solid lines show the median, and the dashed line shows the PPD. The shaded regions denote the 50% and 90% symmetric intervals. In the upper panel, the parametric model is presented with both singular (orange) and non-singular (blue) model configurations. For comparison purposes, the V (very low spin magnitude) model is plotted in a dash-dotted black line. In the bottom panel, the distribution of spin magnitude is inferred over five bins, assuming an either perfectly aligned (pink) or isotropic (green) population. The solid lines denote the median, and the shaded regions denote the central 90% posterior credible bounds. In both cases, the magnitude is consistent within the uncertainties with the parametric (singular and nonsingular) results. The number of bins in the model was chosen to balance resolution with the amount of information in the data; analyses with more bins do not indicate any additional features in the spin distributions.</figcaption></figure>

<table><caption>**Table 6**
Summary of Spin Distribution Models Examined in Section 5.1, with Prior Ranges for the Population Parameters Determining the Spin Models</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Mass<br>Model</th><th colspan="5">Spin Parameters</th></tr><tr><th>E[a]</th><th>Var[a]</th><th>α<sub>a</sub>, β<sub>a</sub></th><th>ζ</th><th>σ<sub>i</sub></th></tr></thead><tbody><tr><td>Gaussian (G)</td><td>C</td><td>[0, 1]</td><td>[0, 0.25]</td><td>≥1</td><td>1</td><td>[0, 4]</td></tr><tr><td>Mixture (M)</td><td>C</td><td>[0, 1]</td><td>[0, 0.25]</td><td>≥1</td><td>[0, 1]</td><td>[0, 4]</td></tr></tbody></table>

**Note.** The fixed parameters are in bold. Each of these distributions is uniform over the stated range, with boundary conditions such that the inferred parameters $\alpha_a, \beta_a$ must be $\ge 1$. Details of the mass model listed here are described in Table 2.

the parametric Mixture model defined in Section 2.4 and using prior distributions shown in Table 6. It marginalizes over all other parameters, including the mass parameters in Model C
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Inferred distribution of cosine spin tilt for the more massive BH for two choices of prior (see Section 2.4). The dashed-dotted line denotes a completely isotropic distribution (see Appendix B). The solid lines show the median. The shaded regions denote the 90% symmetric intervals, and the dashed line denotes the PPD.</figcaption></figure>

and the spin mixture fraction. We observe that spin distributions that decline with increasing magnitude are preferred. In terms of our Beta function parameterization—$E[a]$ and $\text{Var}[a]$, defined in Equation (5)—these have mean spin $E[a] < 1/2$ or equivalently have $\beta_a > \alpha_a$, at posterior probability 0.79. We find that 90% of BH spins in BBHs are less than $a \le 0.55$ from the PPD, and 50% of BH spins are less than $a \le 0.27$. We find similar conclusions if both BH spins are drawn from different distributions (i.e., 90% of BH spins on the more massive BH are less than 0.6). When avoiding singular values in the spin magnitude model distribution, the distribution exhibits a peak structure, i.e., $p(a=0) = p(a=1) = 0$. If allowed to capture the full range of model parameters, including “singular” configurations, the support for small values of $a$ is more pronounced. However, this scenario forces a small—and otherwise observationally unsupported—uptick of probability mass at $a$ near maximal spins. In both cases, the recovered spin distribution in the top panel of Figure 8 is driven by favoring declining spin distributions, which are more compatible with the observed population. This conclusion is also consistent with the preference in Appendix B for the very low spin magnitude model.

We also compute the posterior distribution for the magnitude of BH spins from $\chi_{\text{eff}}$ measurements by modeling the distribution of BH spin magnitudes nonparametrically with five bins, assuming either an isotropic or perfectly aligned population following Farr et al. (2018). We show in the bottom panel of Figure 8 that under the perfectly aligned scenario there is a preference for small BH spin, inferring 90% of BHs to have spin magnitudes below $0.6^{+0.24}_{-0.28}$. However, when spins are assumed to be isotropic, the distribution is relatively flat, with 90% of BH spin magnitudes below $0.8^{+0.15}_{-0.24}$. Thus, the nonparametric analysis produces conclusions consistent with our parametric analyses described above. These conclusions are also reinforced by computing the Bayes factor for a set of fixed-parameter models of spin magnitude and orientation in Appendix B. There we find that the very low spin magnitude model is preferred by a log Bayes factor of 1 or greater in most mass and spin orientation configurations tested (see Figure 13 and Table 8 for details).

Figure 9 shows the inferred distribution of the primary spin tilt for the more massive BH. These results were obtained

without including the effects of component spins on the detection probability; see Appendix A for further discussion. In the Gaussian model ($\zeta = 1$), all BH spin orientations are drawn from spin tilt distributions that are preferentially aligned and parameterized with $\sigma_i$. In that model, the $\sigma_i$ distributions do not differ appreciably from their flat priors. As such, the inferred spin tilt distribution is influenced by large $\sigma_i$, and the result resembles an isotropic distribution. The Mixture distribution does not return a decisive measurement of the mixture fraction, obtaining $\zeta = 0.6^{+0.4}_{-0.5}$. Since the Gaussian model is a subset of the Mixture model, we can compare preferences via the Savage-Dickey ratio. The log Bayes factor for $\zeta = 1$ is $\ln \text{BF} = 0.15$, indicating virtually no preference for any particular orientation distribution. While we allow both BHs to have different typical misalignment, the inference on the second tilt is less informative than the primary. The inferred distribution for $\cos t_2$ is similar to $\cos t_1$ but also closer to the prior.

The mixture fraction distribution is also modeled with the fixed-parameter models in Appendix B. The fixed magnitude distributions considered in Appendix B prefer isotropic to aligned, but the preference is weakened for distributions concentrated at lower spins. A few exceptions occur for the very low spin fixed mass ratio models, with aligned models being slightly preferred.

In general, we are not able to place strong constraints on the distribution of spin orientations. We elaborate in Appendix B.3 on how our BH spin measurements are not yet informative enough to discern between isotropic and aligned orientation distribution via $\chi_{\text{eff}}$.

## 5.2. Interpretation of Spin Distributions

The spins of BHs are affected by a number of uncertain processes that occur during the evolution of the binary. As a consequence, the magnitude distribution is difficult to predict from theoretical models of these processes alone. While the spin of a BH should be related to the rotation of the core of its progenitor star, the amount of spin that is lost during the final stages of the progenitor's life is still highly uncertain. While we have modeled the spins independently, correlations from binary evolution and stellar collapse are possible (Belczynski et al. 2017; Gerosa et al. 2018; Qin et al. 2018; Arca Sedda & Benacquista 2019; Postnov & Kuranov 2019). The core rotational angular momentum before the supernova can be changed from the birth spin of the progenitor by several processes (Langer 2012; de Mink et al. 2013; Amaro-Seoane & Chen 2016). Examples include mass transfer (Packet 1981; Shu & Lubow 1981) and tidal interactions (Petrovic et al. 2005), as well as internal mixing of the stellar layers across the core-envelope boundary via magnetic torquing (Spruit 2002; Maeder & Meynet 2003) and gravity waves (Talon & Charbonnel 2005, 2008; Fuller et al. 2015). In principle, an off-center supernova explosion could also impart significant angular momentum and tilt the spin of the remnant into the collapsing star (Farr et al. 2011a).

Once a BH is formed, however, changing the spin magnitude is more difficult owing to limitations on mass accretion rates affecting how much a BH can be spun up (Thorne 1974; Valsecchi et al. 2010; Wong et al. 2012; Qin et al. 2019). Once the BBH system is formed, the spin magnitudes do not change appreciably over the inspiral (Farr et al. 2014).

No BBH detected to date has a component with confidently high and aligned spin magnitude. The results in the previous
section imply that BHs tend to be born with spin less than our PPD bound of 0.55, or that another process (e.g., supernova kicks or dynamical processes involved in binary formation) induces tilts such that $\chi_{\text{eff}}$ is small.

The possibility of a spin magnitude distribution that peaks at low spins incurs a degeneracy between models that is not easily overcome: when the spin magnitudes are small enough, models produce features that cannot be distinguished within observational uncertainties.

## 6. Discussion and Conclusions

We have presented a variety of estimates for the mass, spin, and redshift distributions of BBH, based on the observed sample of 10 BBH and generic phenomenological population models motivated by electromagnetic observations and theory. Some model independent features are evident from the observations. Notably, no BBHs more massive than GW170729 have been observed to date, but several binaries have component masses likely between 20 and $40 M_{\odot}$. No highly asymmetric (small $q$) system has been observed. Only two systems (GW151226 and GW170729) produce a $\chi_{\text{eff}}$ distribution that is confidently different from zero; conversely, most BH binaries are consistent with $\chi_{\text{eff}}$ near zero. These features drive our inferences about the mass and spin distribution.

Despite exploring a wide range of mass and spin distributions, we find that the BBH merger rate density is $R = 64.0_{-33.0}^{+73.5} \text{ Gpc}^{-3} \text{yr}^{-1}$ for Model A and is within $R = 53.2_{-28.2}^{+55.8} \text{ Gpc}^{-3} \text{yr}^{-1}$ for Models B and C. This result is consistent with the fixed model assumptions reported in the combined O1 and O2 observational periods (Abbott et al. 2018a). We find a significant reduction in the merger rate for BBHs with primary masses larger than $\sim 45 M_{\odot}$. We do not have enough sensitivity to binaries with a BH mass less than $5 M_{\odot}$ to be able to place meaningful constraints on the minimum mass of BHs. We find mild evidence that the mass distribution of coalescing BHs may not be a pure power law, instead being slightly better fit by a model including a broad Gaussian distribution at high mass. We find that the best-fitting models preferentially produce comparable-mass binaries (i.e., $\beta_q > 0$ is preferred).

The mass models in this work supersede results from an older model from O1 that inferred only the power-law index (Abbott et al. 2016c, 2017e). That model found systematically larger values of $\alpha$ than its nearest counterpart in this work, Model A, because the older model used a fixed value for the minimum and maximum mass of 5 and $100 M_{\odot}$, respectively. This extreme $m_{\text{max}}$ is highly disfavored by our current results, and so the older model is also disfavored. Moreover, volumetric sensitivity grows as a strong function of mass. The lack of detections near the older $m_{\text{max}}$ drives a preference for a much smaller maximum BH mass in the new models (Fishbach & Holz 2017). A reduced maximum mass is associated with a shallower power-law fit.

Inferring the redshift distribution is difficult with only a small sample of local events (Fishbach et al. 2018). We have constrained models with extreme variation over redshift, favoring instead those that are uniform in the comoving volume or have increasing merger rates with higher redshift. Many potential formation channels in the literature (Belczynski et al. 2014; Antonini & Rasio 2016; Inayoshi et al. 2016; Mandel & de Mink 2016; Rodriguez et al. 2016b; Bartos et al.

2017; Mapelli et al. 2017; Kruckow et al. 2018) produce event rates that are compatible with those from the previous observing runs (Abbott et al. 2018a) and this work. It is, of course, plausible that several are contributing simultaneously, and no combination of mass, rate, or redshift dependence explored here rules out any of the channels proposed to date. The next generation of interferometers will allow for an exquisite probe into this dependence at large redshifts (Sathya-prakash et al. 2012; Van Den Broeck 2014; Vitale & Farr 2018).

We have modeled the spin distribution in several ways, forming inferences on the spin magnitude and tilt distributions. In all of our analysis, the evidence disfavors distributions with large spin components aligned (or nearly aligned) with the orbital angular momentum; specifically, we find that 90% of the spin magnitude PPD is smaller than 0.55. We cannot significantly constrain the degree of spin-orbit misalignment in the population. However, regardless of the mass or assumed spin tilt distribution, there is a preference (demonstrated in Figure 8 and Appendix B) for distributions that emphasize lower spin magnitudes. Our inferences suggest that 90% of coalescing BH binaries are formed with $\chi_{\text{eff}} < 0.3$. Low spins argue against so-called second-generation mergers, where at least one of the components of the binary is a BH formed from a previous merger (Berti et al. 2007; González et al. 2007) and possesses spins near 0.7 (Fishbach et al. 2017).

GW170729 is notable in several ways: it is the most massive, largest $\chi_{\text{eff}}$, and most distant redshift event detected so far. To quantify the impact it has on our results, where possible we have presented model posteriors that reflect its presence in or exclusion from the event set. Many of our predictions are robust despite its extreme values—by far, and not unexpectedly, its influence is most significant in the distribution of $m_{\text{max}}$. It also impacts our conclusions about redshift evolution, where its absence flattens the inferred redshift evolution.

Recent modeling using only the first six released events (Wysocki et al. 2018; Roulet & Zaldarriaga 2019) has come to similar conclusions about low spin magnitudes and the shape of the power-law distribution. The presence of an apparent upper limit to the merging BBH mass distribution was also observed after the first six released events (Fishbach & Holz 2017). An enhancement that will benefit these types of analyses in the future is a simultaneous fit of the astrophysical model and its parameters and noise background model (Gaebel et al. 2019).

Several studies have noted that population features (Mandel & O'Shaughnessy 2010; Stevenson et al. 2015, 2017a; Farr et al. 2017, 2018; Fishbach et al. 2017; Fishbach & Holz 2017; Gerosa & Berti 2017; Kovetz et al. 2017; Talbot & Thrane 2017, 2018; Zevin et al. 2017; Barrett et al. 2018; Gerosa et al. 2018; Wysocki et al. 2018) and complementary physics (Abbott et al. 2016f; Stevenson et al. 2017a; Zevin et al. 2017; Chen et al. 2018) will be increasingly accessible as observations accumulate. Additional events will also permit the enhancement of the simple phenomenological models used in this work and comparison with modeling of astrophysical processes. Given the event merger rates estimated here and anticipated improvements in sensitivity (Abbott et al. 2018b), hundreds of BBHs and tens of binary neutron stars are expected to be collected in the operational lifetime of second-generation GW instruments. Thus, the inventory of BBH in the coming
years will enable inquiries into astrophysics that were previously unobtainable.

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS), and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies, as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación, the Vicepresidència i Conselleria d'Innovació Recerca i Turisme and the Conselleria d'Educació i Universitat del Govern de les Illes Balears, the Conselleria d'Educació Investigació Cultura i Esport de la Generalitat Valenciana, the National Science Centre of Poland, the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, the Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, Innovations, and Communications, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan, and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS, and the State of Niedersachsen/Germany for provision of computational resources.

## Appendix A
Systematics

In this section, we discuss the systematic uncertainties that affect our analysis and show that they are subdominant to statistical uncertainties. We focus on two major sources of systematic uncertainty. The first of these is introduced by the waveform models that are used to extract the parameters of

<table><caption>**Table 7**
Summary of Intervals for Each of the Parameters Considered in the Models of Sections 3–5</caption><thead><tr><th>Parameter/Model</th><th>Reference</th><th>Spin &lang;VT&rang;</th><th>using SEOBNR</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center; font-weight: bold;">Mass</td></tr><tr><td>&alpha; (Model A)</td><td>0.4<sup>+1.9</sup><sub>-1.9</sub></td><td>0.4<sup>+1.3</sup><sub>-1.9</sub></td><td>0.3<sup>+1.4</sup><sub>-2.0</sub></td></tr><tr><td>&alpha; (Model B)</td><td>1.3<sup>+1.4</sup><sub>-1.7</sub></td><td>1.2<sup>+1.4</sup><sub>-1.7</sub></td><td>1.2<sup>+1.4</sup><sub>-1.8</sub></td></tr><tr><td>&alpha; (Model C)</td><td>7.1<sup>+4.4</sup><sub>-4.8</sub></td><td>...</td><td>7.3<sup>+4.2</sup><sub>-4.8</sub></td></tr><tr><td>&beta;<sub>q</sub> (Model B)</td><td>6.9<sup>+4.6</sup><sub>-5.7</sub></td><td>7.0<sup>+4.6</sup><sub>-5.6</sub></td><td>7.1<sup>+4.4</sup><sub>-6.0</sub></td></tr><tr><td>&beta;<sub>q</sub> (Model C)</td><td>4.5<sup>+6.6</sup><sub>-5.2</sub></td><td>...</td><td>5.0<sup>+6.3</sup><sub>-5.7</sub></td></tr><tr><td>m<sub>max</sub> (Model A)</td><td>41.6<sup>+9.6</sup><sub>-4.3</sub></td><td>41.6<sup>+9.9</sup><sub>-4.4</sub></td><td>41.3<sup>+9.2</sup><sub>-4.2</sub></td></tr><tr><td>m<sub>max</sub> (Model B)</td><td>40.8<sup>+11.8</sup><sub>-4.4</sub></td><td>40.6<sup>+10.8</sup><sub>-4.3</sub></td><td>40.5<sup>+12.5</sup><sub>-3.9</sub></td></tr><tr><td>m<sub>max</sub> (Model C)</td><td>62.0<sup>+34.0</sup><sub>-28.7</sub></td><td>...</td><td>62.5<sup>+33.8</sup><sub>-29.1</sub></td></tr><tr><td>m<sub>min</sub> (Model B)</td><td>7.8<sup>+1.2</sup><sub>-2.5</sub></td><td>7.7<sup>+1.3</sup><sub>-2.4</sub></td><td>7.7<sup>+1.3</sup><sub>-2.4</sub></td></tr><tr><td>m<sub>min</sub> (Model C)</td><td>6.9<sup>+1.7</sup><sub>-2.8</sub></td><td>...</td><td>6.9<sup>+1.7</sup><sub>-2.7</sub></td></tr><tr><td>&lambda;<sub>m</sub> (Model C)</td><td>0.3<sup>+0.4</sup><sub>-0.2</sub></td><td>...</td><td>0.3<sup>+0.4</sup><sub>-0.2</sub></td></tr><tr><td>&mu;<sub>m</sub> (Model C)</td><td>29.8<sup>+5.8</sup><sub>-7.3</sub></td><td>...</td><td>30.1<sup>+5.8</sup><sub>-7.4</sub></td></tr><tr><td>&sigma;<sub>m</sub> (Model C)</td><td>6.4<sup>+2.2</sup><sub>-4.2</sub></td><td>...</td><td>5.9<sup>+3.5</sup><sub>-4.0</sub></td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Rate</td></tr><tr><td>R (Model A)</td><td>64.0<sup>+73.5</sup><sub>-33.0</sub></td><td>62.8<sup>+74.0</sup><sub>-33.3</sub></td><td>62.4<sup>+74.0</sup><sub>-31.9</sub></td></tr><tr><td>R (Model B)</td><td>53.2<sup>+55.8</sup><sub>-28.2</sub></td><td>51.8<sup>+55.3</sup><sub>-26.9</sub></td><td>52.9<sup>+52.7</sup><sub>-28.2</sub></td></tr><tr><td>R (Model C)</td><td>58.3<sup>+72.3</sup><sub>-32.2</sub></td><td>...</td><td>58.0<sup>+70.3</sup><sub>-32.0</sub></td></tr><tr><td colspan="4" style="text-align: center; font-weight: bold;">Spin</td></tr><tr><td>cos t<sub>1</sub> (Model A)</td><td>2.0<sup>+1.8</sup><sub>-1.4</sub></td><td>2.2<sup>+1.6</sup><sub>-1.6</sub></td><td>1.2<sup>+2.4</sup><sub>-1.0</sub></td></tr><tr><td>cos t<sub>1</sub> (Model B)</td><td>2.1<sup>+1.7</sup><sub>-1.5</sub></td><td>2.2<sup>+1.6</sup><sub>-1.6</sub></td><td>1.4<sup>+2.2</sup><sub>-1.2</sub></td></tr><tr><td>cos t<sub>1</sub> (Model C)</td><td>1.8<sup>+1.9</sup><sub>-1.7</sub></td><td>...</td><td>1.2<sup>+2.4</sup><sub>-1.1</sub></td></tr><tr><td>cos t<sub>2</sub> (Model A)</td><td>2.3<sup>+1.5</sup><sub>-1.6</sub></td><td>2.4<sup>+1.5</sup><sub>-1.6</sub></td><td>2.1<sup>+1.7</sup><sub>-1.6</sub></td></tr><tr><td>cos t<sub>2</sub> (Model B)</td><td>2.3<sup>+1.5</sup><sub>-1.5</sub></td><td>2.4<sup>+1.5</sup><sub>-1.6</sub></td><td>2.0<sup>+1.8</sup><sub>-1.5</sub></td></tr><tr><td>cos t<sub>2</sub> (Model C)</td><td>2.0<sup>+1.8</sup><sub>-1.6</sub></td><td>...</td><td>2.0<sup>+1.8</sup><sub>-1.6</sub></td></tr></tbody></table>

**Note.** The reference uses the posteriors derived from the IMRPhenomPv2 waveform model and without spin effects included in $\langle VT \rangle$. The second column allows for spin effects in $\langle VT \rangle$ estimation. Finally, the third column shows the population model parameters inferred when the SEOBNRv3 waveform model is used to derive the event posteriors. Spin enabled $\langle VT \rangle$ is only available for Models A and B, but we expect that Model C would exhibit similar trends. Broadly, the mass and rate parameters are nearly the same and well within their respective uncertainties with and without spin effects in $\langle VT \rangle$, as well as considering the SEOBNRv3 waveform model. The most notable difference comes from the parameterized spin distribution. The differences are primarily related to the spin tilt distribution, and, for clarity, we suppress the spin magnitude distribution and mixture parameters, which are nearly identical.

individual events, and the second is in the estimation of the detection efficiency.

### A.1. Waveform Systematics

In Abbott et al. (2018a), two waveform families are used to extract the parameters of individual events: SEOBNRv3 (Pan et al. 2014; Babak et al. 2017) and IMRPhenomPv2 (Hannam et al. 2014; Husa et al. 2016; Khan et al. 2016). While both families capture a wide variety of physical effects, including simple precession and other spin effects, they do not match each other exactly over the whole of the parameter space. Differences between the waveforms can therefore lead to slight biases in the inference of individual events’ parameters and thereby impact the inferred population distributions. To directly assess the impact of these uncertainties on our results, we have repeated our calculations using parameter estimates based on
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Inferred distribution of cosine spin tilt for the more massive BH for two choices of prior (see Section 2.4) with the SEOBNRv3 waveform model, with the same definitions as in Figure 9.</figcaption></figure>

SEOBNRv3 (the results in the main text all use IMRPhenomPv2). See Table 7—we find that the two waveform models produce at most modestly different inferences about key parameters. For example, the standard Model B mass and spin distribution analysis with SEOBNRv3 leads us to infer that the 90% upper bound of $a_1$ is 0.5 and credible intervals on $m_{max}$ and $R$ are 36.7–53.0 $M_\odot$ and 24.8–105.7 Gpc$^{-3}$ yr$^{-1}$, which is consistent with the IMRPhenomPv2 Model B estimates of 0.6, 36.4–52.6 $M_\odot$, and 24.9–109.0 Gpc$^{-3}$ yr$^{-1}$ presented in the main text. Similarly, in the redshift evolution analysis, we infer that the redshift evolution parameter $\lambda = 8.4_{-9.7}^{+9.6}$ under the SEOBNRv3 waveform compared to $\lambda = 8.4_{-9.5}^{+9.6}$ under the IMRPhenomPv2 waveform. In summary, for the mass and rate part of the distributions, we find that there is no significant change whatsoever. The most significant change is to the parameter controlling the primary tilt angle. The SEOBNRv3 waveform predicts that this parameter is smaller by 50%, along with a smaller reduction in the secondary tilt angle parameter. Compare Figure 10, produced with SEOBNRv3 derived event posteriors, with Figure 9. SEOBNRv3 produces a distribution of tilts that is closer to isotropic than its IMRPhenomPv2 counterpart. However, the models are compatible to within their uncertainties over the distribution $p(\cos t_1)$.

## A.2. Selection Effects and Sensitive Volume

In this subsection we detail the various assumptions and possible systematics that enter into our calculation of the detection efficiency. The detectability of a BBH merger in GWs depends on the distance and orientation of the binary along with its intrinsic parameters, especially its component masses. In order to model the underlying population and determine the BBH merger rate, we must properly model the mass-, redshift-, and spin-dependent selection effects and incorporate them into our population analysis according to Equation (11). One way to infer the sensitivity of the detector network to a given population of BBH mergers is by carrying out large-scale simulations in which synthetic GW waveforms are injected into the detector data and subsequently searched for. The parameters of the injected waveforms can be drawn directly from the fixed population of interest, or alternatively, the injections can be placed to more broadly cover parameter space and reweighed to match the properties of the population (Tiwari 2018). Such injection campaigns were carried out in Abbott et al. (2018a) to measure the total sensitive spacetime

volume $\langle VT \rangle$ and the corresponding merger rate for two fixed-parameter populations (power law and flat-in-log). However, it is computationally expensive to carry out an injection campaign that sufficiently covers the multidimensional population hyperparameter space considered in this work. For this reason, for the parametric population studies in this work, we employ a semianalytic method to estimate the fraction of found detections as a function of masses, spins, and redshift (or equivalently, distance).

Our estimates of the network sensitivity are based on the semianalytic method that was used to infer the BBH mass distribution from the first four GW detections (Abbott et al. 2016c, 2017e). This method assumes that a BBH system is detectable if and only if it produces a signal-to-noise ratio (S/N) $\rho \ge \rho_{th}$ in a single detector, where the threshold S/N, $\rho_{th}$, is typically chosen to be 8. Given a BBH system with known component masses, spins, and cosmological redshift, and a detector with stationary Gaussian noise characterized by a given power spectral density (PSD), one can calculate the optimal S/N, $\rho_{opt}$, of the signal emitted by the BBH merger. The optimal S/N corresponds to the S/N of the signal produced by a face-on, directly overhead BBH merger with the same masses, spins, and redshift. Given $\rho_{opt}$, the distribution of single-detector S/Ns $\rho$—corresponding to sources with random orientations with respect to the detector—can be calculated using the analytic distribution of angular factors $\Theta \equiv \rho/\rho_{opt}$ (Finn & Chernoff 1993). Under these assumptions, the probability of detecting a system of given masses, spins, and redshift, $P_{det}(m_1, m_2, \chi_1, \chi_2, z)$, is given by the probability that $\rho \ge \rho_{th}$, or equivalently, that a randomly drawn $\Theta \ge \rho_{th}/\rho_{opt}(m_1, m_2, \chi_1, \chi_2, z)$. $P_{det}$ referred to in this section is equivalent to the $f(z|\xi)$ that appears in Equation (12) of Section 2.

The semianalytic calculation relies on two main simplifying assumptions: the detection threshold $\rho_{th}$, and the choice of PSD for characterizing the detector noise. When fitting the mass distribution to the first four BBH events in Abbott et al. (2017e), we assumed that the PSD in each LIGO interferometer could be approximated by the Early High Sensitivity curve in Abbott et al. (2018b) during O1 and the first few months of O2, and we fixed $\rho_{th} = 8$. We refer to the sensitivity estimate under these assumptions as the raw semianalytic calculation. In reality, the detector PSD fluctuates throughout the observing period. Additionally, the fixed detection threshold on S/N does not directly account for the empirical distributions of astrophysical and noise triggers and does not have a direct
correspondence with the detection statistic used by the GW searches to rank significance of triggers (Messick et al. 2017; Nitz et al. 2017; Abbott et al. 2018a). Consequently, the sensitive spacetime volume of a population estimated using an S/N threshold may differ from the one obtained using injections, which has a threshold on the pipeline-dependent detection statistic.

We therefore pursue two modifications to the raw semi-analytic calculation in order to reduce the bias in our sensitivity estimates and the resulting population estimates. We emphasize that these modifications do not noticeably affect the inferred shape of the population, e.g., the mass power-law slope, but do lead to different rate estimates, reflecting a systematic uncertainty in the inferred merger rate and its evolution with redshift that, given the small number of events and uncertainty in the phenomenological population models, remains subdominant to the statistical uncertainty. This is explicitly shown in the remainder of this section.

In the first modification, which we employ throughout the mass distribution analysis (Section 3), we calibrate the raw semianalytic method to the injection campaign in Abbott et al. (2018a). The calibration takes the form of mass-dependent calibration factors, calculated by least-squares regression as described below; see Wysocki & O'Shaughnessy (2018) for relevant data products. Specifically, we use injections to evaluate $\langle VT \rangle_i \equiv \int d\xi p(\xi|\theta_i) VT_{\text{true}}(\xi)$ for a set of reference hyperparameters $\theta_i$ (here, mass distribution models with different exponents $\alpha$ and maximum masses $m_{\text{max}}$), where $\xi$ denotes all binary parameters. To calculate $\langle VT \rangle_i$ from injections into the PyCBC detection pipeline, we consider injections to be “detected” if they have a detection statistic $\rho \ge 8$, where $\rho$ is the statistic used in the PyCBC analysis of O2 data (Nitz et al. 2017; Abbott et al. 2018a). This is comparable to the detection statistic $\rho = 8.7$ of the lowest-significance GW event included in our analysis, GW170729. Note that, as discussed in Section 4, because we adopt a fixed detection threshold, our analysis differs from the rate analysis in Abbott et al. (2018a), which does not fix a detection threshold, instead assigning to each trigger a probability of astrophysical origin (Farr et al. 2015b). Once we have computed $\langle VT \rangle_i$, we correct the raw semianalytic model $VT_{\text{raw}}$ described above by a factor $f(\xi)$, which is a low-order polynomial in $\xi$: $f(\xi) = \sum_{\alpha} \lambda_{\alpha} F_{\alpha}(\xi)$, with $F_{\alpha}$ the relevant basis polynomials. We minimize the mean-square difference between $\langle VT \rangle$ as computed by injections and $\int d\xi f(\xi) \langle VT \rangle_{\text{raw}}(\xi) p(\xi|\theta)$. If $H$ is the precomputed matrix of weight “moments” $H_{k,\alpha} = \int d\xi p(\xi|\theta_k) VT(\xi) F_{\alpha}(\xi)$, then the coefficients of this least-squares expression can be computed analytically as $\lambda = (H^T \gamma H)^{-1} H^T \gamma \langle VT \rangle$, where $\gamma$ is a diagonal inverse covariance matrix characterizing the Monte Carlo integration errors of each individual $\langle VT \rangle_i$. This procedure yields the mass-calibrated sensitive volume $\langle VT \rangle_{\text{cal}}$.

The top panel of Figure 11 shows the comparison between the raw semianalytic $\langle VT \rangle$, the calibrated $\langle VT \rangle$, and the injection $\langle VT \rangle$ across the 2D hyperparameter space of Model A for the mass distribution. We have repeated our mass distribution analysis with different choices of the $\langle VT \rangle$ calibration and found that the effect on the shape of the mass distribution and the overall merger rate $R$ are much smaller than the differences between Models A, B, and C and the statistical errors associated with a small sample of 10 events.

As shown in Figure 11, the main effect of this calibration is to decrease $\langle VT \rangle$ by a factor of $\sim 1.6$. Over the relevant part of parameter space (i.e., the regions of the $\alpha-m_{\text{max}}$ plane that have likelihood support), this factor remains fairly constant, implying that the inferred shape of the mass distribution is not affected by applying the $\langle VT \rangle$ calibration, although the overall rate is increased by about a factor of $\sim 1.6$ compared to the raw semianalytic calculation. We have verified this explicitly by repeating the analysis with and without calibrated $\langle VT \rangle$.

For the redshift evolution analysis (Section 4), it is not sufficient to calibrate the mass dependence of the detection probability; we must verify that the semianalytic calculation reproduces the proper redshift dependence. Therefore, we pursue an alternative modification to the raw semianalytic calculation. In this modification, we replace the single PSD of the raw semianalytic calculation with a different PSD calculated for the Livingston detector for each 5-day chunk of observing time in O1 and O2. We find that this assumption correctly reproduces the redshift-dependent sensitivity empirically determined by the injection campaigns into the GstLAL pipeline for two fixed mass distributions (see Figure 12), whereas adopting different assumptions, such as using the PSDs calculated for the Hanford detector instead of the Livingston detector, or changing the single-detector S/N threshold away from 8, yields curves in Figure 12 that deviate noticeably from the distribution of recovered injections. This modification to the sensitivity calculation is necessary in the redshift analysis because the detection probability can fluctuate significantly at high redshifts $z > 0.5$, where there is a very small probability of detection but considerable physical volume. Due to computational cost, the number of detections available at high redshift is insufficient to directly calibrate the redshift-dependent detection probability to injections as we did in the mass distribution section.

We find that between the two methods we use to estimate detection efficiency, the effect on the inferred mass distribution is negligible. However, the second time-varying approach employed in the redshift analysis underpredicts the overall merger rate by $\sim 70\%$ compared to the first calibrated approach (see the bottom right panel of Figure 11). This reflects a systematic uncertainty in the high-redshift detection efficiency and the implied merger rate. When additional detections lead to improved statistical constraints on the merger rate across redshift, it will become increasingly necessary to place a very large number of injections at high redshift and closely spaced in time in order to accurately estimate the high-redshift sensitivity.

Another difference between the mass distribution analysis presented in Section 3 and the redshift evolution analysis of Section 4 is in the treatment of BBH spins. Section 3 marginalizes over the spin distribution and includes first-order spin effects in the calculation of $\langle VT \rangle$, while the redshift analysis of Section 4 does not. From Table 7, we find that including first-order spin effects in the calculation of $P_{\text{det}}$ and the corresponding sensitive spacetime volume $\langle VT \rangle$ results in mostly indistinguishable population estimates compared to neglecting spin entirely. Similarly, fixing the spin distribution does not appreciably affect the inferred mass distribution. Therefore, for simplicity we neglect the effect of spin distribution in the redshift evolution analysis. Meanwhile, the effects of spin on the sensitive volume $\langle VT \rangle$ do have a moderate influence on inferences about the spin tilt angles,
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Ratio between the raw semianalytic computation of $\langle VT \rangle$ and the $\langle VT \rangle$ computed by injections into the PyCBC search pipeline (top left panel), and the same ratio for the mass-calibrated $\langle VT \rangle$ (top right panel), for different mass distributions described by the two-parameter Model A. The bottom panel shows the same ratios, but this time comparing the $\langle VT \rangle$ derived with the time-varying PSDs applied to the semianalytic calculation $VT_{psd}$ to the raw and the calibrated $\langle VT \rangle$. The $\langle VT \rangle$ for the injections is calculated for a threshold of $\varrho = 8.0$, where $\varrho$ is the signal-to-noise ratio model statistic used in the PyCBC analysis of O2 data. This threshold roughly matches the detection statistic of the lowest-significance detection, GW170729, which has $\varrho = 8.7$ in PyCBC. We use the mass-calibrated $VT_{cal}$ for the parametric mass and spin distribution analyses in Sections 3 and 5 in order to better match the injection results, while the redshift evolution analysis uses $VT_{psd}$ in order to carefully track the sensitivity at high redshift. The discrepancy between the methods may be due to the limited number of high-redshift injections. However, the difference between all three methods is relatively constant as a function of the mass population, particularly where posterior support for the mass distribution hyperparameters is high, indicating that systematic uncertainties in the $\langle VT \rangle$ estimation do not have a large impact on our results.</figcaption></figure>

presented in Section 5.1. When considering the effects of spin with $\langle VT \rangle$, there is about a 10% shift in the median spin tilt angle parameters inferred, but this is well within the much wider credible interval. Therefore, such effects do not change our overall astrophysical conclusions, and their influence on the results shown is comparable to what would result from

different priors on the population parameters (e.g., choosing a different prior range of $\sigma_i$ as compared to Table 6).

We also note that all our calculations of the detection efficiency are based on the IMRPhenomPv2 waveform. Differences between the phasing and, more importantly, the amplitude of the waveform can lead to different S/Ns and
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Redshift distribution of injections recovered with a false-alarm rate (FAR) less than 0.1 yr⁻¹ by the search pipeline GstLAL for the two fixed-parameter injection sets, power law (red) and flat-in-log (green), compared to the expectation from the semianalytic calculation used for the redshift evolution analysis, as described in the text. The underlying redshift distributions of the injected populations are assumed to follow a uniform in comoving volume and source-frame time distribution. The FAR threshold of 0.1 yr⁻¹ nearly matches the FAR of the lowest-significance GW event, GW170729, with an FAR of 0.18 yr⁻¹ in the GstLAL pipeline. The semianalytic calculation closely predicts the redshift distribution of the found injections.</figcaption></figure>

detection statistics for the same sets of physical parameters. To bound the significance of this effect, we carry out the injection-based $\langle VT \rangle$ estimation for both the IMRPhenomPv2 and SEOBNRv2 waveforms and find that for populations described by the two-parameter mass Model A, the waveforms produce $\langle VT \rangle$ estimates consistent to 10% across the relevant region of hyperparameter space with high posterior probability. Therefore, compared to the statistical uncertainties, the choice in waveform does not contribute a significant systematic uncertainty for the $\langle VT \rangle$ estimation.

Finally, an additional systematic uncertainty we have neglected in the $\langle VT \rangle$ and parametric rate calculations is the calibration uncertainty. While the event posterior samples have incorporated a marginalization over uncertainties on the calibration (Farr et al. 2015a) for both strain amplitude and phase, the $\langle VT \rangle$ estimation here does not. The amplitude calibration uncertainty results in an 18% volume uncertainty (Abbott et al. 2018a), which is currently below the level of statistical uncertainty in our population-averaged merger rate estimate.

## Appendix B
### Alternative Spin Models

We perform here a number of complementary analyses to reinforce the robustness of the results in Section 5 and gauge the effect of fixed-parameter choices on spin inferences. Instead of a parameterized model such as those used in Section 5, we focus on a few discrete choices of model parameters to reinforce the conclusions in that section. These choices provide a complementary view to the results presented earlier and also display our current ability (or inability) to measure features in differing parts of the mass and spin parameter space.

#### B.1. *Model Selection*

We choose a set of specific realizations of the general model described in Section 2.2, building on Farr et al. (2017) and Tiwari et al. (2018). Four discrete spin magnitude models are considered, the first three being special cases of Equation (4):

1. Low (L): $p(a) = 2(1-a)$, i.e., $\alpha_a = 1, \beta_a = 2$.

2. Flat (F): $p(a) = 1$, i.e., $\alpha_a = 1, \beta_a = 1$.

3. High (H): $p(a) = 2a$, i.e., $\alpha_a = 2, \beta_a = 1$.

4. Very low (V): $p(a) \propto e^{-(a/0.2)}$.

Such magnitude distributions are chosen as simple representations of low, moderate, and highly spinning individual BHs. The very low (V) population is added to capture the features of an even lower spinning population—this is motivated by the features at low spin of the parametric distribution displayed in Figure 8.

For spin orientations we consider three fixed models representing extreme cases of Equation (6):

1. Isotropic (I): $p(\cos t_i) = 1/2$; $-1 < \cos t_i < 1$, i.e., $\zeta = 0$.

2. Aligned (A): $p(\cos t_i) = \delta(\cos t_i - 1)$, i.e., $\zeta = 1, \sigma_i = 0$.

3. Restricted (R): $p(\cos t_i) = 1$; $0 < \cos t_i < 1$. This is the same as I, except the spins are restricted to point above the orbital plane.

The isotropic distribution is motivated by dynamical or similarly disordered assembly scenarios, while the aligned one better captures a population of isolated binaries, under the simplifying assumption that the stars remain perfectly aligned throughout their evolution. In order to assess any preferences in the data for binaries with $\chi_{eff} > 0$, we introduce the restricted model: it resembles the isotropic distribution but limits tilt angles to be positive. While we have mathematically defined the R model by assuming tilted spins, the same $\chi_{eff}$ distribution can be generated with nonprecessing spins.

Here we perform our inference entirely through $\chi_{eff}$, whose 12 different distributions are illustrated in Figure 13. Since we do not have conclusive results on $\beta_q$ from Figure 3, we cannot make a single simplifying assumption on the mass model, which the $\chi_{eff}$ distribution depends on. We therefore consider three limiting cases: two of these fix the mass ratio to fiducial values, $q = 1$ and $q = 0.5$. The third corresponds to a fixed-parameter model with $\alpha = 1, m_{min} = 5, m_{max} = 50$. Figure 13 illustrates the $\chi_{eff}$ distributions implied by each of these scenarios.

Following Farr et al. (2017) and Tiwari et al. (2018), we calculate the evidence and compute the Bayes factors for each of the zero-dimensional spin models. Results are provided in Table 8, with the low and isotropic distribution (LI) as the reference.

Because of degeneracies in the GW waveform between mass ratio and $\chi_{eff}$, the choice of mass distribution impacts inferences about spins. This effect explains the significant
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Top row: $p(\chi_{\text{eff}})$ under various model assumptions. Labels in each subpanel legend correspond to the tilt and magnitude models defined in Appendix B.1. Isotropic models (left) provide support for both negative and positive $\chi_{\text{eff}}$. Aligned models (middle) assume perfect alignment for each of the four magnitude distributions. Restricted models (right) have the same shape as the Isotropic ones, with support over $\chi_{\text{eff}} > 0$ only. However, they can be generated with nonprecessing spins. Bottom row: posterior on the mixture fraction $\zeta$ between isotropic and aligned distributions. $\zeta = 0$ corresponds to a completely isotropic distribution.</figcaption></figure>

<table><caption>**Table 8**

Natural Log Bayes Factors for Various Spin Distributions</caption><thead><tr><th></th><th>Very low</th><th>Low</th><th>Flat</th><th>High</th></tr></thead><tbody><tr><td><strong>Isotropic</strong></td><td>1.29</td><td>0.0</td><td>-1.04</td><td>-2.25</td></tr><tr><td><strong>Restricted</strong></td><td>3.5</td><td>3.22</td><td>1.06</td><td>-0.2</td></tr><tr><td><strong>Aligned</strong></td><td>1.39</td><td>-4.57</td><td>-13.62</td><td>-33.13</td></tr><tr><td><strong>q = 0.5</strong></td><td>Very low</td><td>Low</td><td>Flat</td><td>High</td></tr><tr><td>Isotropic</td><td>1.32</td><td>0.0</td><td>-1.12</td><td>-2.6</td></tr><tr><td>Restricted</td><td>3.55</td><td>3.23</td><td>1.0</td><td>-0.58</td></tr><tr><td>Aligned</td><td>1.52</td><td>-4.15</td><td>-12.86</td><td>-31.6</td></tr><tr><td><strong>Fixed param.</strong></td><td>Very low</td><td>Low</td><td>Flat</td><td>High</td></tr><tr><td>Isotropic</td><td>0.64</td><td>0.0</td><td>-2.0</td><td>-3.85</td></tr><tr><td>Restricted</td><td>1.83</td><td>0.8</td><td>-2.3</td><td>-5.0</td></tr><tr><td>Aligned</td><td>-2.69</td><td>-11.98</td><td>-21.98</td><td>-44.6</td></tr></tbody></table>

**Note.** The orientation models are described in Section 2. We find modest evidence for small spins. When spins are small, we cannot make strong statements about the distribution of spin orientations.

difference in Bayes factors for the third row in the table. We find again that our result moderately favors small BH spins. The restricted models with $\chi_{\text{eff}}$ strictly positive consistently produce the highest Bayes factors. For the small-spin magnitude models we cannot make strong statements about the distribution of spin orientations. Models containing highly spinning components are significantly disfavored, with high or flat aligned spins particularly selected against (e.g., FA and HA are disfavored with Bayes factors ranging in $[10^{-11}, 10^{-6}]$ and

$[10^{-21}, 10^{-13}]$, respectively). As a bracket for our uncertainty on the mass and mass ratio distribution, we evaluated the Bayes factors for the fixed-parameter model $\alpha = 2.3$, $m_{\text{min}} = 5$, $m_{\text{max}} = 50$. They differ from the third mass model in Table 8 by a factor comparable to unity.

## B.2. Spin Mixture Models

The models considered for model selection in Table 8 all assume a fixed set of spin magnitudes and tilts. There is no reason to believe, however, that the universe produces from only one of these distributions. A natural extension is to allow for a mixing fraction describing the relative abundances of perfectly aligned and isotropically distributed BH spins.

We assume that the aligned and isotropic components follow the same spin magnitude distribution. It is possible that BHs with a different distribution of spin orientations would have a different distribution of spin magnitudes, but given our weaker constraints on spin magnitudes, we focus on spin tilts sharing the same magnitude distribution.

We compute the posterior on the fraction of aligned binaries $\zeta$ in the population as per Equation (6) in the limit ($\sigma_i \to 0$). The models here are subsets of the Mixture distribution, with purely isotropic being $\zeta = 0$ and completely aligned being $\zeta = 1$. The prior on the mixing fraction is flat.

All of the models that contain a completely aligned component favor isotropy over alignment. This ability to distinguish a mixing fraction diminishes with smaller spin magnitudes. This is because such spin magnitudes yield populations that are not distinguishable to within measurement
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Posterior distribution for the fraction of informative binaries (i.e., $|\chi_{eff}| > 0.05$), and the fraction of those informative binaries with positive $\chi_{eff}$ (i.e., $\chi_{eff} > 0.05$).</figcaption></figure>

uncertainty of $\chi_{eff}$. We do not include the most-favored restricted (R) configuration, but we expect that the results would be similar. Coupled with the model selection results in the previous section, this implies that the mixing fraction is not well determined when fixed to the models (low and very low) that are favored by the data (see Figure 13). As stated above, in this case our ability to measure the mixing fraction is negligible.

### B.3. Three-bin Analysis of $\chi_{eff}$

We illustrate here how $\chi_{eff}$ measurements can provide insights into discerning spin orientation distributions. Following Farr et al. (2018), we split the range of $\chi_{eff}$ into three bins. One encompasses the fraction of uninformative binaries with $\chi_{eff}$ consistent with zero ($|\chi_{eff}| \le 0.05$); the vertical axis of Figure 14 shows the fraction of binaries lying outside of this bin. The other two capture significantly positive ($\chi_{eff} > 0.05$) and significantly negative ($\chi_{eff} < -0.05$) binaries. The width 0.05 is chosen to be of the order of the uncertainty in a typical event posterior.

The aligned spin scenario is preferred in the posterior support on the right half of Figure 14: the small fraction of binaries that are informative tend to possess $\chi_{eff}$ greater than zero. Conversely, if the spins are isotropic, there would be no preference for positive or negative $\chi_{eff}$, and the posterior in Figure 14 would peak toward the middle. However, of the 10 observed binaries, 8 are consistent with zero $\chi_{eff}$ and only 2 are informative, thus demonstrating that our ability to distinguish between the two scenarios is weak.

## Appendix C
### Importance Resampling the Single-event Likelihood

Our hierarchical population analysis uses the individual-event likelihood for each event $n = 1, \dots, N$, $\mathcal{L}(d_n|\xi, z)$ (see Section 2, Equation (11)). Individual-event analyses report posterior samples drawn from a density that is proportional to this likelihood times a prior (Veitch et al. 2015; Abbott et al. 2018a). The prior density used is uniform in detector-frame masses and proportional to the square of the luminosity distance (Veitch et al. 2015); in terms of the source-frame

masses and redshift, the prior is

$$ \begin{aligned} p(m_1^{\text{source}}, m_2^{\text{source}}, z) &\propto d_L^2(z) \frac{\partial m_1^{\text{det}}}{\partial m_1^{\text{source}}} \frac{\partial m_2^{\text{det}}}{\partial m_2^{\text{source}}} \frac{\partial d_L}{\partial z} \\ &= d_L^2(z)(1+z)^2 \frac{\partial d_L}{\partial z}. \end{aligned} \quad (18) $$

The derivative of the luminosity distance in a spatially flat universe (Hogg 1999) is

$$ \frac{\partial d_L}{\partial z} = \frac{d_L}{1+z} + (1+z) \frac{d_H}{E(z)}, \quad (19) $$

where $d_H = c/H_0$ is the Hubble distance and

$$ E(z) \equiv \sqrt{\Omega_M(1+z)^2 + \Omega_\Lambda} \quad (20) $$

in a $\Lambda$CDM universe.

Given a set of posterior samples as described above, we can transform them to samples from the likelihood over source-frame masses and redshift by importance resampling with weights that are the inverse prior

$$ w(m_1^{\text{source}}, m_2^{\text{source}}, z) = \frac{1}{p(m_1^{\text{source}}, m_2^{\text{source}}, z)}. \quad (21) $$

The integral in Equation (11) may then be approximated as

$$ \int d\xi \, dz \, \mathcal{L}(d_n|\xi, z) \frac{dN}{d\xi dz} \propto \sum_{w}^{N_{\text{sample}=1}} (m_{1,i}^{\text{source}}, m_{2,i}^{\text{source}}, z_i) \frac{dN}{d\xi_i dz_i}, \quad (22) $$

where the sum is taken over all posterior samples, and we assume that the population distribution is expressed in terms of source-frame masses (i.e., $\xi = (m_1^{\text{source}}, m_2^{\text{source}}, \dots)$). Alternately, we can construct a random resampling of the set of existing posterior samples, with sample $i$ appearing in the resampling with probability proportional to $w(m_{1,i}^{\text{source}}, m_{2,i}^{\text{source}}, z_i)$; the integral is then proportional to the average value of $dN/d\xi dz$ over the resampled set. The (unknown) constant of proportionality is related to the Bayesian evidence for event $n$; as long as a consistent method (weighted sum or resampling) is used to compare different population models, this constant is irrelevant to computing Bayes factors between models.

## References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, PhRvL, 116, 131103

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016b, PhRvL, 116, 241103

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016c, PhRvX, 6, 041015

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016d, ApJL, 818, L22

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016e, ApJL, 833, L1

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016f, PhRvL, 116, 131102

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJL, 851, L35

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, PhRvL, 119, 161101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, ApJL, 848, L13

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017d, CQGra, 34, 104002

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017e, PhRvL, 118, 221101

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018a, PRX, submitted

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2018b, LRR, 21, 3

Acernese, F., Agathos, M., Agatsuma, K., et al. 2015, CQGra, 32, 024001

Ajith, P., Hannam, M., Husa, S., et al. 2011, PhRvL, 106, 241101

Alexander, K. D., Berger, E., Fong, W., et al. 2017, ApJL, 848, L21

Ali-Haïmoud, Y., Kovetz, E. D., & Kamionkowski, M. 2017, PhRvD, 96, 123523
Amaro-Seoane, P., & Chen, X. 2016, MNRAS, 458, 3075

Ando, K., Inomata, K., Kawasaki, M., Mukaida, K., & Yanagida, T. T. 2018, PhRvD, 97, 123512

Antonini, F., Murray, N., & Mikkola, S. 2014, ApJ, 781, 45

Antonini, F., & Perets, H. B. 2012, ApJ, 757, 27

Antonini, F., & Rasio, F. A. 2016, ApJ, 831, 187

Antonini, F., Toonen, S., & Hamers, A. S. 2017, ApJ, 841, 77

Arca Sedda, M., & Benacquista, M. 2019, MNRAS, 482, 2991

Asensio Ramos, A. 2018, Bayesian Astrophysics, Canary Islands Winter School of Astrophysics (Cambridge: Cambridge Univ. Press)

Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36

Babak, S., Taracchini, A., & Buonanno, A. 2017, PhRvD, 95, 024010

Bai, Y., Barger, V., & Lu, S. 2018, arXiv:1802.04909

Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367

Baird, E., Fairhurst, S., Hannam, M., & Murphy, P. 2013, PhRvD, 87, 024035

Banerjee, S. 2017, MNRAS, 467, 524

Barkat, Z., Rakavy, G., & Sack, N. 1967, PhRvL, 18, 379

Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, MNRAS, 477, 4685

Bartos, I., Kocsis, B., Haiman, Z., & Márka, S. 2017, ApJ, 835, 165

Belczynski, K., Bulik, T., Fryer, C. L., et al. 2010, ApJ, 714, 1217

Belczynski, K., Buonanno, A., Cantiello, M., et al. 2014, ApJ, 789, 120

Belczynski, K., Heger, A., Gladysz, W., et al. 2016a, A&A, 594, A97

Belczynski, K., Holz, D. E., Bulik, T., & O'Shaughnessy, R. 2016b, Natur, 534, 512

Belczynski, K., Kalogera, V., & Bulik, T. 2002, ApJ, 572, 407

Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223

Belczynski, K., Klencki, J., Meynet, G., et al. 2017, arXiv:1706.07053

Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504

Berti, E., Cardoso, V., Gonzalez, J. A., et al. 2007, PhRvD, 76, 064034

Bethe, H. A., & Brown, G. E. 1998, ApJ, 506, 780

Bird, S., Cholis, I., Muñoz, J. B., et al. 2016, PhRvL, 116, 201301

Bond, J. R., Arnett, W. D., & Carr, B. J. 1984, ApJ, 280, 825

Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127

Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115

Byrnes, C. T., Hindmarsh, M., Young, S., & Hawkins, M. R. S. 2018, JCAP, 2018, 041

Carr, B., Kühnel, F., & Sandstad, M. 2016, PhRvD, 94, 083504

Carr, B. J., & Hawking, S. W. 1974, MNRAS, 168, 399

Chatterjee, S., Rodriguez, C. L., Kalogera, V., & Rasio, F. A. 2017, ApJL, 836, L26

Chen, H.-Y., Fishbach, M., & Holz, D. E. 2018, Natur, 562, 545

Chen, Z.-C., & Huang, Q.-G. 2018, ApJ, 864, 61

Chornock, R., Berger, E., Kasen, D., et al. 2017, ApJL, 848, L19

Chruslinska, M., Belczynski, K., Klencki, J., & Benacquista, M. 2018, MNRAS, 474, 2937

Clausen, D., Sigurdsson, S., & Chernoff, D. F. 2013, MNRAS, 428, 3618

Clesse, S., & García-Bellido, J. 2017, PDU, 15, 142

Corsaro, E., Lee, Y.-N., García, R. A., et al. 2017, NatAs, 1, 0064

Coughlin, M., Meyers, P., Thrane, E., Luo, J., & Christensen, N. 2015, PhRvD, 91, 063004

Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Sci, 358, 1556

Cowperthwaite, P. S., Berger, E., Villar, V. A., et al. 2017, ApJL, 848, L17

Damour, T. 2001, PhRvD, 64, 124013

de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243

de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166

de Mink, S. E., & Mandel, I. 2016, MNRAS, 460, 3545

Dewi, J. D. M., Podsiadlowski, P., & Sena, A. 2006, MNRAS, 368, 1742

Dominik, M., Belczynski, K., Fryer, C., et al. 2012, ApJ, 759, 52

Dominik, M., Belczynski, K., Fryer, C., et al. 2013, ApJ, 779, 72

Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2010, MNRAS, 407, 1946

Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2011, MNRAS, 416, 133

Eldridge, J. J., & Stanway, E. R. 2016, MNRAS, 462, 3302

Ertl, T., Janka, H. T., Woosley, S. E., Sukhbold, T., & Ugliano, M. 2016, ApJ, 818, 124

Farr, B., Holz, D. E., & Farr, W. M. 2018, ApJL, 854, L9

Farr, B., Kalogera, V., & Luijten, E. 2014, PhRvD, 90, 024014

Farr, W. M., Farr, B., & Littenberg, T. 2015a, Modelling Calibration Errors In CBC Waveforms, Tech. Rep. LIGO-T1400682, LIGO Scientific Collaboration and Virgo Collaboration, https://dcc.ligo.org/LIGO-T1400682/public

Farr, W. M., Gair, J. R., Mandel, I., & Cutler, C. 2015b, PhRvD, 91, 023005

Farr, W. M., Kremer, K., Lyutikov, M., & Kalogera, V. 2011a, ApJ, 742, 81

Farr, W. M., Sravan, N., Cantrell, A., et al. 2011b, ApJ, 741, 103

Farr, W. M., Stevenson, S., Miller, M. C., et al. 2017, Natur, 548, 426

Ferrari, S., & Cribari-Neto, F. 2004, Journal of Applied Statistics, 31, 799

Finn, L. S., & Chernoff, D. F. 1993, PhRvD, 47, 2198

Fishbach, M., & Holz, D. E. 2017, ApJL, 851, L25

Fishbach, M., Holz, D. E., & Farr, B. 2017, ApJL, 840, L24

Fishbach, M., Holz, D. E., & Farr, W. M. 2018, ApJL, 863, L41

Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64

Fowler, W. A., & Hoyle, F. 1964, ApJS, 9, 201

Fragione, G., Ginsburg, I., & Kocsis, B. 2018, ApJ, 856, 92

Fregeau, J. M. 2004, PhD thesis, Massachusetts Institute of Technology, http://hdl.handle.net/1721.1/29454

Freire, P. C. C., Ransom, S. M., Bégin, S., et al. 2008, ApJ, 675, 670

Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, ApJ, 749, 91

Fuller, J., Cantiello, M., Lecoanet, D., & Quataert, E. 2015, ApJ, 810, 101

Gaebel, S. M., Veitch, J., Dent, T., & Farr, W. M. 2019, MNRAS, 484, 4008

Gelman, A., Carlin, J. B., Stern, H. S., et al. 2004, Bayesian Data Analysis (Boca Raton, Fl: Chapman and Hall/CRC), doi:10.1103/PhysRev.136. B1224

Georg, J., & Watson, S. 2017, JHEP, 2017, 138

Gerosa, D., & Berti, E. 2017, PhRvD, 95, 124046

Gerosa, D., Berti, E., O'Shaughnessy, R., et al. 2018, PhRvD, 98, 084036

Gerosa, D., Kesden, M., Sperhake, U., Berti, E., & O'Shaughnessy, R. 2015, PhRvD, 92, 064016

Giacobbo, N., & Mapelli, M. 2018, MNRAS, 480, 2011

Giacobbo, N., Mapelli, M., & Spera, M. 2018, MNRAS, 474, 2959

Goldstein, A., Veres, P., Burns, E., et al. 2017, ApJL, 848, L14

González, J. A., Sperhake, U., Brügmann, B., Hannam, M., & Husa, S. 2007, PhRvL, 98, 091101

Gräfener, G., & Hamann, W. R. 2008, A&A, 482, 945

Grandclément, P., Ihm, M., Kalogera, V., & Belczynski, K. 2004, PhRvD, 69, 102002

Grindlay, J., Portegies Zwart, S., & McMillan, S. 2006, NatPh, 2, 116

Hannam, M., Schmidt, P., Bohé, A., et al. 2014, PhRvL, 113, 151101

Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288

Heger, A., & Woosley, S. E. 2002, ApJ, 567, 532

Hilbe, J. M., de Souza, R. S., & Ishida, E. E. O. 2017, Bayesian Models for Astrophysical Data: Using R, JAGS, Python, and Stan (Cambridge: Cambridge Univ. Press)

Hinder, I., Kidder, L. E., & Pfeiffer, H. P. 2018, PhRvD, 98, 044015

Hinder, I., Vaishnav, B., Herrmann, F., Shoemaker, D. M., & Laguna, P. 2008, PhRvD, 77, 081502

Hogg, D. W. 1999, arXiv:astro-ph/9905116

Hogg, D. W., Myers, A. D., & Bovy, J. 2010, ApJ, 725, 2166

Huerta, E. A., Kumar, P., Agarwal, B., et al. 2017, PhRvD, 95, 024038

Huerta, E. A., Kumar, P., McWilliams, S. T., O'Shaughnessy, R., & Yunes, N. 2014, PhRvD, 90, 084016

Husa, S., Khan, S., Hannam, M., et al. 2016, PhRvD, 93, 044006

Inayoshi, K., Kashiyama, K., Visbal, E., & Haiman, Z. 2016, MNRAS, 461, 2722

Inomata, K., Kawasaki, M., Mukaida, K., Tada, Y., & Yanagida, T. T. 2017, PhRvD, 95, 123510

Ivanova, N., Heinke, C. O., Rasio, F. A., Belczynski, K., & Fregeau, J. M. 2008, MNRAS, 386, 553

Janka, H.-T. 2012, ARNPS, 62, 407

Kalogera, V. 2000, ApJ, 541, 319

Kalogera, V., Belczynski, K., Kim, C., O'Shaughnessy, R., & Willems, B. 2007, PhR, 442, 75

Khan, S., Husa, S., Hannam, M., et al. 2016, PhRvD, 93, 044007

Kimpson, T. O., Spera, M., Mapelli, M., & Ziosi, B. M. 2016, MNRAS, 463, 2443

Klein, A., Boetzel, Y., Gopakumar, A., Jetzer, P., & de Vittori, L. 2018, PhRvD, 98, 104043

Kocsis, B., & Levin, J. 2012, PhRvD, 85, 123005

Kovetz, E. D., Cholis, I., Breysse, P. C., & Kamionkowski, M. 2017, PhRvD, 95, 103010

Kreidberg, J., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, ApJ, 757, 36

Kremer, K., Chatterjee, S., Ye, C. S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 871, 38

Kruckow, M. U., Tauris, T. M., Langer, N., Kramer, M., & Izzard, R. G. 2018, MNRAS, 481, 1908

Kudritzki, R.-P., & Puls, J. 2000, ARA&A, 38, 613

Kulkarni, S. R., Hut, P., & McMillan, S. 1993, Natur, 364, 421

Langer, N. 2012, ARA&A, 50, 107
LIGO Scientific Collaboration, Aasi, J., Abbott, B. P., et al. 2015, CQGra, 32, 074001

Littenberg, T. B., Farr, B., Coughlin, S., Kalogera, V., & Holz, D. E. 2015, ApJL, 807, L24

Liu, B., & Lai, D. 2018, ApJ, 863, 68

Loredo, T. J. 2004, in AIP Conf. Ser. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint (Melville, NY: AIP), 195

Lower, M. E., Thrane, E., Lasky, P. D., & Smith, R. 2018, PhRvD, 98, 083028

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maeder, A., & Meynet, G. 2003, A&A, 411, 543

Mandel, I. 2010, PhRvD, 81, 084029

Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634

Mandel, I., Farr, W. M., Colonna, A., et al. 2017, MNRAS, 465, 3254

Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086

Mandel, I., Haster, C.-J., Dominik, M., & Belczynski, K. 2015, MNRAS, 450, L85

Mandel, I., & O'Shaughnessy, R. 2010, CQGra, 27, 114007

Mandic, V., Bird, S., & Cholis, I. 2016, PhRvL, 117, 201102

Mapelli, M. 2016, MNRAS, 459, 3432

Mapelli, M., Colpi, M., & Zampieri, L. 2009, MNRAS, 395, L71

Mapelli, M., & Giacobbo, N. 2018, MNRAS, 479, 4391

Mapelli, M., Giacobbo, N., Ripamonti, E., & Spera, M. 2017, MNRAS, 472, 2422

Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50

Marchant, P., Renzo, M., Farmer, R., et al. 2018, arXiv:1810.13412

Margalit, B., & Metzger, B. D. 2017, ApJL, 850, L19

Margutti, R., Berger, E., Fong, W., et al. 2017, ApJL, 848, L20

McKernan, B., Ford, K. E. S., Lyra, W., & Perets, H. B. 2012, MNRAS, 425, 460

Mennekens, N., & Vanbeveren, D. 2014, A&A, 564, A134

Messick, C., Blackburn, K., Brady, P., et al. 2017, PhRvD, 95, 042001

Misner, C. W., Thorne, K. S., & Wheeler, J. A. 1973, Gravitation Physics Series (1st ed.; San Francisco, CA: Freeman)

Morscher, M., Umbreit, S., Farr, W. M., & Rasio, F. A. 2013, ApJL, 763, L15

Mortlock, D. J., Feeney, S. M., Peiris, H. V., Williamson, A. R., & Nissanke, S. M. 2018, arXiv:1811.11723

Nicholl, M., Berger, E., Kasen, D., et al. 2017, ApJL, 848, L18

Nitz, A. H., Dent, T., Dal Canton, T., Fairhurst, S., & Brown, D. A. 2017, ApJ, 849, L18

Ober, W. W., El Eid, M. F., & Fricke, K. J. 1983, A&A, 119, 61

O'Connor, E., & Ott, C. D. 2011, ApJ, 730, 70

O'Leary, R. M., Meiron, Y., & Kocsis, B. 2016, ApJL, 824, L12

O'Leary, R. M., Rasio, F. A., Fregeau, J. M., Ivanova, N., & O'Shaughnessy, R. 2006, ApJ, 637, 937

O'Shaughnessy, R., Gerosa, D., & Wysocki, D. 2017, PhRvL, 119, 011101

O'Shaughnessy, R., Kalogera, V., & Belczynski, K. 2010, ApJ, 716, 615

Özel, F., & Freire, P. 2016, ARA&A, 54, 401

Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918

Özel, F., Psaltis, D., Narayan, R., & Santos Villarreal, A. 2012, ApJ, 757, 55

Packet, W. 1981, A&A, 102, 17

Pan, Y., Buonanno, A., Taracchini, A., et al. 2014, PhRvD, 89, 084006

Peters, P. C. 1964, PhRv, 136, 1224

Petrovic, J., Langer, N., & van der Hucht, K. A. 2005, A&A, 435, 1013

Petrovich, C., & Antonini, F. 2017, ApJ, 846, 146

Pian, E., D'Avalzo, P., Benetti, S., et al. 2017, Natur, 551, 67

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Portegies Zwart, S. F., & McMillan, S. L. W. 2000, ApJL, 528, L17

Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431

Portegies Zwart, S. F., & Yungelson, L. R. 1998, A&A, 332, 173

Postnov, K. A., & Kuranov, A. G. 2019, MNRAS, 483, 3288

Postnov, K. A., & Yungelson, L. R. 2014, LRR, 17, 3

Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28

Qin, Y., Marchant, P., Fragos, T., Meynet, G., & Kalogera, V. 2019, ApJL, 870, L18

Quinlan, G. D., & Shapiro, S. L. 1987, ApJ, 321, 199

Racine, É. 2008, PhRvD, 78, 044021

Rakavy, G., & Shaviv, G. 1967, ApJ, 148, 803

Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., et al. 2018a, PhRvD, 98, 123005

Rodriguez, C. L., Amaro-Seoane, P., Chatterjee, S., & Rasio, F. A. 2018b, PhRvL, 120, 151101

Rodriguez, C. L., & Antonini, F. 2018, ApJ, 863, 7

Rodriguez, C. L., Chatterjee, S., & Rasio, F. A. 2016a, PhRvD, 93, 084029

Rodriguez, C. L., Haster, C.-J., Chatterjee, S., Kalogera, V., & Rasio, F. A. 2016b, ApJL, 824, L8

Rodriguez, C. L., & Loeb, A. 2018, ApJL, 866, L5

Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, PhRvL, 115, 051101

Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016c, ApJL, 832, L2

Roulet, J., & Zaldarriaga, M. 2019, MNRAS, 484, 4216

Sadowski, A., Belczynski, K., Bulik, T., et al. 2008, ApJ, 676, 1162

Samsing, J., MacLeod, M., & Ramirez-Ruiz, E. 2014, ApJ, 784, 71

Sasaki, M., Suyama, T., Tanaka, T., & Yokoyama, S. 2016, PhRvL, 117, 061101

Sathyaprakash, B., Abernathy, M., Acernese, F., et al. 2012, CQGra, 29, 124013

Savchenko, V., Ferrigno, C., Kuulkers, E., et al. 2017, ApJL, 848, L15

Scheepmaker, R. A., Haas, M. R., Gieles, M., et al. 2007, A&A, 469, 925

Schmidt, P., Ohme, F., & Hannam, M. 2015, PhRvD, 91, 024043

Shu, F. H., & Lubow, S. H. 1981, ARA&A, 19, 277

Sigurdsson, S., & Hernquist, L. 1993, Natur, 364, 423

Skilling, J. 2004, in AIP Conf. Ser. 735, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395 doi:10.1063/1.1835238

Soares-Santos, M., Holz, D. E., Annis, J., et al. 2017, ApJL, 848, L16

Spera, M., & Mapelli, M. 2017, MNRAS, 470, 4739

Spera, M., Mapelli, M., & Bressan, A. 2015, MNRAS, 451, 4086

Spruit, H. C. 2002, A&A, 381, 923

Stevenson, S., Berry, C. P. L., & Mandel, I. 2017a, MNRAS, 471, 2801

Stevenson, S., Ohme, F., & Fairhurst, S. 2015, ApJ, 810, 58

Stevenson, S., Vigna-Gómez, A., Mandel, I., et al. 2017b, NatCo, 8, 14906

Stone, N. C., Metzger, B. D., & Haiman, Z. 2017, MNRAS, 464, 946

Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38

Talbot, C., & Thrane, E. 2017, PhRvD, 96, 023012

Talbot, C., & Thrane, E. 2018, ApJ, 856, 173

Talon, S., & Charbonnel, C. 2005, A&A, 440, 981

Talon, S., & Charbonnel, C. 2008, A&A, 482, 597

Taracchini, A., Buonanno, A., Pan, Y., et al. 2014, PhRvD, 89, 061502

Tauris, T. M., Kramer, M., Freire, P. C. C., et al. 2017, ApJ, 846, 170

Thorne, K. S. 1974, ApJ, 191, 507

Thorne, K. S. 1983, in Gravitational Radiation, ed. N. Deruelle & T. Piran (Amsterdam: North-Holland Publishing Co.), 1

Tiwari, V. 2018, CQGra, 35, 145009

Tiwari, V., Fairhurst, S., & Hannam, M. 2018, ApJ, 868, 140

Troja, E., Piro, L., van Eerten, H., et al. 2017, Natur, 551, 71

Tutukov, A. V., & Yungelson, L. R. 1993, MNRAS, 260, 675

Ugliano, M., Janka, H.-T., Marek, A., & Arcones, A. 2012, ApJ, 757, 69

Valsecchi, F., GlebbEEK, E., Farr, W. M., et al. 2010, Natur, 468, 77

Vanbeveren, D. 2009, NewAR, 53, 27

Van Den Broeck, C. 2014, JPhCS, 484, 012008

Veitch, J., Del Pozzo, W., Cody, & Pitkin, M. 2017, johnweitch/cpnest: Minor Optimisation, version 0.1.4, Zenodo, doi:10.5281/zenodo.835874

Veitch, J., Raymond, V., Farr, B., et al. 2015, PhRvD, 91, 042003

Vink, J. S., & de Koter, A. 2005, A&A, 442, 587

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Vitale, S., & Farr, W. M. 2018, arXiv:1808.00901

Vitale, S., Lynch, R., Sturani, R., & Graff, P. 2017, CQGra, 34, 03LT01

Voss, R., & Tauris, T. M. 2003, MNRAS, 342, 1169

Wong, T.-W., Valsecchi, F., Fragos, T., & Kalogera, V. 2012, ApJ, 747, 111

Woosley, S. E. 2017, ApJ, 836, 244

Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Natur, 450, 390

Wysocki, D., Lange, J., & O'Shaughnessy, R. 2018, arXiv:1805.06442

Wysocki, D., & O'Shaughnessy, R. 2018, Rescaling VT Factors to Match Tabulated VT Averages, Tech. Rep. LIGO-T1800427, LIGO Scientific Collaboration and Virgo Collaboration, https://dcc.ligo.org/ LIGO-T1800427/public

Zevin, M., Pankow, C., Rodriguez, C. L., et al. 2017, ApJ, 846, 82

Ziosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703