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 MSUB SUB and a powerlaw index of   1.3SUB1.7SUBSUP1.4SUP 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.2SUB28.2SUBSUP55.8SUP GpcSUP3SUP yrSUP1SUP 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.
['black holes', 'gravitational waves', 'statistical', 'Astrophysics - High Energy Astrophysical Phenomena']
Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo
{'No Header': "Binary Black Hole Population Properties Inferred from the First and Second Observing Runs of Advanced LIGO and Advanced Virgo  \n-  \nThis is the Accepted Manuscript version of an article accepted for publication in The Astrophysical Journal Letters and published as ApJL, 882, L24 (2019). IOP Publishing Ltd is not responsible for any errors or omissions in this version of the manuscript or any version derived from it. The Version of Record is available online at https://doi.org/10.3847/2041-8213/ab3800.  \nB. P. Abbott, 1 R. Abbott, 1 T. D. Abbott, 2 S. Abraham, 3 F. Acernese, 4, 5 K. Ackley, 6 C. Adams, 7 R. X. Adhikari, 1 V. B. Adya, 8, 9 C. Affeldt, 8, 9 M. Agathos, 10 K. Agatsuma, 11 N. Aggarwal, 12 O. D. Aguiar, 13 L. Aiello, 14, 15 A. Ain, 3 P. Ajith, 16 G. Allen, 17 A. Allocca, 18, 19 M. A. Aloy, 20 P. A. Altin, 21 A. Amato, 22 A. Ananyeva, 1 S. B. Anderson, 1 W. G. Anderson, 23 S. V. Angelova, 24 S. Antier, 25 S. Appert, 1 K. Arai, 1 M. C. Araya, 1 J. S. Areeda, 26 M. Ar'ene, 27 N. Arnaud, 25, 28 K. G. Arun, 29 S. Ascenzi, 30, 31 G. Ashton, 6 S. M. Aston, 7 P. Astone, 32 F. Aubin, 33 P. Aufmuth, 9 K. AultONeal, 34 C. Austin, 2 V. Avendano, 35 A. Avila-Alvarez, 26 S. Babak, 36, 27 P. Bacon, 27 F. Badaracco, 14, 15 M. K. M. Bader, 37 S. Bae, 38 P. T. Baker, 39 F. Baldaccini, 40, 41 G. Ballardin, 28 S. W. Ballmer, 42 S. Banagiri, 43 J. C. Barayoga, 1 S. E. Barclay, 44 B. C. Barish, 1 D. Barker, 45 K. Barkett, 46 S. Barnum, 12 F. Barone, 4, 5 B. Barr, 44 L. Barsotti, 12 M. Barsuglia, 27 D. Barta, 47 J. Bartlett, 45 I. Bartos, 48 R. Bassiri, 49 A. Basti, 18, 19 M. Bawaj, 50, 41 J. C. Bayley, 44 M. Bazzan, 51, 52 B. B'ecsy, 53 M. Bejger, 27, 54 I. Belahcene, 25 A. S. Bell, 44 D. Beniwal, 55 B. K. Berger, 49 G. Bergmann, 8, 9 S. Bernuzzi, 56, 57 J. J. Bero, 58 C. P. L. Berry, 59 D. Bersanetti, 60 A. Bertolini, 37 J. Betzwieser, 7 R. Bhandare, 61 J. Bidler, 26 I. A. Bilenko, 62 S. A. Bilgili, 39 G. Billingsley, 1 J. Birch, 7 R. Birney, 24 O. Birnholtz, 58 S. Biscans, 1, 12 S. Biscoveanu, 6 A. Bisht, 9 M. Bitossi, 28, 19 M. A. Bizouard, 25 J. K. Blackburn, 1 C. D. Blair, 7 D. G. Blair, 63 R. M. Blair, 45 S. Bloemen, 64 N. Bode, 8, 9 M. Boer, 65 Y. Boetzel, 66 G. Bogaert, 65 F. Bondu, 67 E. Bonilla, 49 R. Bonnand, 33 P. Booker, 8, 9 B. A. Boom, 37 C. D. Booth, 68 R. Bork, 1 V. Boschi, 28 S. Bose, 69, 3 K. Bossie, 7 V. Bossilkov, 63 J. Bosveld, 63 Y. Bouffanais, 27 A. Bozzi, 28 C. Bradaschia, 19 P. R. Brady, 23 A. Bramley, 7 M. Branchesi, 14, 15 J. E. Brau, 70 T. Briant, 71 J. H. Briggs, 44 F. Brighenti, 72, 73 A. Brillet, 65 M. Brinkmann, 8, 9 V. Brisson, 25, ∗ P. Brockill, 23 A. F. Brooks, 1 D. D. Brown, 55 S. Brunett, 1 A. Buikema, 12 T. Bulik, 74 H. J. Bulten, 75, 37 A. Buonanno, 36, 76 R. Buscicchio, 11 D. Buskulic, 33 C. Buy, 27 R. L. Byer, 49 M. Cabero, 8, 9 L. Cadonati, 77 G. Cagnoli, 22, 78 C. Cahillane, 1 J. Calder'on Bustillo, 6 T. A. Callister, 1 E. Calloni, 79, 5 J. B. Camp, 80 W. A. Campbell, 6 M. Canepa, 81, 60 K. C. Cannon, 82 H. Cao, 55 J. Cao, 83 E. Capocasa, 27 F. Carbognani, 28 S. Caride, 84 M. F. Carney, 59 G. Carullo, 18 J. Casanueva Diaz, 19 C. Casentini, 30, 31 S. Caudill, 37 M. Cavagli'a, 85 F. Cavalier, 25 R. Cavalieri, 28 G. Cella, 19 P. Cerd'a-Dur'an, 20 G. Cerretani, 18, 19 E. Cesarini, 86, 31 O. Chaibi, 65 K. Chakravarti, 3 S. J. Chamberlin, 87 M. Chan, 44 S. Chao, 88 P. Charlton, 89 E. A. Chase, 59 E. Chassande-Mottin, 27 D. Chatterjee, 23 M. Chaturvedi, 61 K. Chatziioannou, 90 B. D. Cheeseboro, 39 H. Y. Chen, 91 X. Chen, 63 Y. Chen, 46 H.-P. Cheng, 48 C. K. Cheong, 92 H. Y. Chia, 48 A. Chincarini, 60 A. Chiummo, 28 G. Cho, 93 H. S. Cho, 94 M. Cho, 76 N. Christensen, 65, 95 Q. Chu, 63 S. Chua, 71 K. W. Chung, 92 S. Chung, 63 G. Ciani, 51, 52 A. A. Ciobanu, 55 R. Ciolfi, 96, 97 F. Cipriano, 65 A. Cirone, 81, 60 F. Clara, 45 J. A. Clark, 77 P. Clearwater, 98 F. Cleva, 65 C. Cocchieri, 85 E. Coccia, 14, 15 P.-F. Cohadon, 71 D. Cohen, 25 R. Colgan, 99 M. Colleoni, 100 C. G. Collette, 101 C. Collins, 11 L. R. Cominsky, 102 M. Constancio Jr., 13 L. Conti, 52 S. J. Cooper, 11 P. Corban, 7 T. R. Corbitt, 2 I. Cordero-Carri'on, 103 K. R. Corley, 99 N. Cornish, 53 A. Corsi, 84 S. Cortese, 28 C. A. Costa, 13 R. Cotesta, 36 M. W. Coughlin, 1 S. B. Coughlin, 68, 59 J.-P. Coulon, 65 S. T. Countryman, 99 P. Couvares, 1 P. B. Covas, 100 E. E. Cowan, 77 D. M. Coward, 63 M. J. Cowart, 7 D. C. Coyne, 1 R. Coyne, 104 J. D. E. Creighton, 23 T. D. Creighton, 105 J. Cripe, 2 M. Croquette, 71 S. G. Crowder, 106 T. J. Cullen, 2 A. Cumming, 44 L. Cunningham, 44 E. Cuoco, 28 T. Dal Canton, 80 G. D'alya, 107 S. L. Danilishin, 8, 9 S. D'Antonio, 31 K. Danzmann, 9, 8 A. Dasgupta, 108 C. F. Da Silva Costa, 48 L. E. H. Datrier, 44 V. Dattilo, 28 I. Dave, 61 M. Davier, 25 D. Davis, 42 E. J. Daw, 109 D. DeBra, 49 M. Deenadayalan, 3 J. Degallaix, 22 M. De Laurentis, 79, 5 S. Del'eglise, 71 W. Del Pozzo, 18, 19 L. M. DeMarchi, 59 N. Demos, 12 T. Dent, 8, 9, 110 R. De Pietri, 111, 57 J. Derby, 26 R. De Rosa, 79, 5 C. De Rossi, 22, 28 R. DeSalvo, 112 O. de Varona, 8, 9 S. Dhurandhar, 3 M. C. D'ıaz, 105 T. Dietrich, 37 L. Di Fiore, 5 M. Di Giovanni, 113, 97 T. Di Girolamo, 79, 5 A. Di Lieto, 18, 19 B. Ding, 101 S. Di Pace, 114, 32 I. Di Palma, 114, 32 F. Di Renzo, 18, 19 A. Dmitriev, 11 Z. Doctor, 91 F. Donovan, 12 K. L. Dooley, 68, 85 S. Doravari, 8, 9 I. Dorrington, 68 T. P. Downes, 23 M. Drago, 14, 15 J. C. Driggers, 45 Z. Du, 83 J.-G. Ducoin, 25 P. Dupej, 44 S. E. Dwyer, 45 P. J. Easter, 6 T. B. Edo, 109 M. C. Edwards, 95 A. Effler, 7 P. Ehrens, 1 J. Eichholz, 1 S. S. Eikenberry, 48 M. Eisenmann, 33 R. A. Eisenstein, 12 R. C. Essick, 91 H. Estelles, 100 D. Estevez, 33 Z. B. Etienne, 39 T. Etzel, 1 M. Evans, 12 T. M. Evans, 7 V. Fafone, 30, 31, 14 H. Fair, 42 S. Fairhurst, 68 X. Fan, 83 S. Farinon, 60 B. Farr, 70 W. M. Farr, 11 E. J. Fauchon-Jones, 68 M. Favata, 35 M. Fays, 109 M. Fazio, 115 C. Fee, 116 J. Feicht, 1 M. M. Fejer, 49 F. Feng, 27 A. Fernandez-Galiana, 12 I. Ferrante, 18, 19 E. C. Ferreira, 13  \nT. A. Ferreira, 13 F. Ferrini, 28 F. Fidecaro, 18, 19 I. Fiori, 28 D. Fiorucci, 27 M. Fishbach, 91 R. P. Fisher, 42, 117 J. M. Fishner, 12 M. Fitz-Axen, 43 R. Flaminio, 33, 118 M. Fletcher, 44 E. Flynn, 26 H. Fong, 90 J. A. Font, 20, 119 P. W. F. Forsyth, 21 J.-D. Fournier, 65 S. Frasca, 114, 32 F. Frasconi, 19 Z. Frei, 107 A. Freise, 11 R. Frey, 70 V. Frey, 25 P. Fritschel, 12 V. V. Frolov, 7 P. Fulda, 48 M. Fyffe, 7 H. A. Gabbard, 44 B. U. Gadre, 3 S. M. Gaebel, 11 J. R. Gair, 120 L. Gammaitoni, 40 M. R. Ganija, 55 S. G. Gaonkar, 3 A. Garcia, 26 C. Garc'ıa-Quir'os, 100 F. Garufi, 79, 5 B. Gateley, 45 S. Gaudio, 34 G. Gaur, 121 V. Gayathri, 122 G. Gemme, 60 E. Genin, 28 A. Gennai, 19 D. George, 17 J. George, 61 L. Gergely, 123 V. Germain, 33 S. Ghonge, 77 Abhirup Ghosh, 16 Archisman Ghosh, 37 S. Ghosh, 23 B. Giacomazzo, 113, 97 J. A. Giaime, 2, 7 K. D. Giardina, 7 A. Giazotto, 19, † K. Gill, 34 G. Giordano, 4, 5 L. Glover, 112 P. Godwin, 87 E. Goetz, 45 R. Goetz, 48 B. Goncharov, 6 G. Gonz'alez, 2 J. M. Gonzalez Castro, 18, 19 A. Gopakumar, 124 M. L. Gorodetsky, 62 S. E. Gossan, 1 M. Gosselin, 28 R. Gouaty, 33 A. Grado, 125, 5 C. Graef, 44 M. Granata, 22 A. Grant, 44 S. Gras, 12 P. Grassia, 1 C. Gray, 45 R. Gray, 44 G. Greco, 72, 73 A. C. Green, 11, 48 R. Green, 68 E. M. Gretarsson, 34 P. Groot, 64 H. Grote, 68 S. Grunewald, 36 P. Gruning, 25 G. M. Guidi, 72, 73 H. K. Gulati, 108 Y. Guo, 37 A. Gupta, 87 M. K. Gupta, 108 E. K. Gustafson, 1 R. Gustafson, 126 L. Haegel, 100 O. Halim, 15, 14 B. R. Hall, 69 E. D. Hall, 12 E. Z. Hamilton, 68 G. Hammond, 44 M. Haney, 66 M. M. Hanke, 8, 9 J. Hanks, 45 C. Hanna, 87 M. D. Hannam, 68 O. A. Hannuksela, 92 J. Hanson, 7 T. Hardwick, 2 K. Haris, 16 J. Harms, 14, 15 G. M. Harry, 127 I. W. Harry, 36 C.-J. Haster, 90 K. Haughian, 44 F. J. Hayes, 44 J. Healy, 58 A. Heidmann, 71 M. C. Heintze, 7 H. Heitmann, 65 P. Hello, 25 G. Hemming, 28 M. Hendry, 44 I. S. Heng, 44 J. Hennig, 8, 9 A. W. Heptonstall, 1 Francisco Hernandez Vivanco, 6 M. Heurs, 8, 9 S. Hild, 44 T. Hinderer, 128, 37, 129 D. Hoak, 28 S. Hochheim, 8, 9 D. Hofman, 22 A. M. Holgado, 17 N. A. Holland, 21 K. Holt, 7 D. E. Holz, 91 P. Hopkins, 68 C. Horst, 23 J. Hough, 44 E. J. Howell, 63 C. G. Hoy, 68 A. Hreibi, 65 E. A. Huerta, 17 D. Huet, 25 B. Hughey, 34 M. Hulko, 1 S. Husa, 100 S. H. Huttner, 44 T. Huynh-Dinh, 7 B. Idzkowski, 74 A. Iess, 30, 31 C. Ingram, 55 R. Inta, 84 G. Intini, 114, 32 B. Irwin, 116 H. N. Isa, 44 J.-M. Isac, 71 M. Isi, 1 B. R. Iyer, 16 K. Izumi, 45 T. Jacqmin, 71 S. J. Jadhav, 130 K. Jani, 77 N. N. Janthalur, 130 P. Jaranowski, 131 A. C. Jenkins, 132 J. Jiang, 48 D. S. Johnson, 17 A. W. Jones, 11 D. I. Jones, 133 R. Jones, 44 R. J. G. Jonker, 37 L. Ju, 63 J. Junker, 8, 9 C. V. Kalaghatgi, 68 V. Kalogera, 59 B. Kamai, 1 S. Kandhasamy, 85 G. Kang, 38 J. B. Kanner, 1 S. J. Kapadia, 23 S. Karki, 70 K. S. Karvinen, 8, 9 R. Kashyap, 16 M. Kasprzack, 1 S. Katsanevas, 28 E. Katsavounidis, 12 W. Katzman, 7 S. Kaufer, 9 K. Kawabe, 45 N. V. Keerthana, 3 F. K'ef'elian, 65 D. Keitel, 44 R. Kennedy, 109 J. S. Key, 134 F. Y. Khalili, 62 H. Khan, 26 I. Khan, 14, 31 S. Khan, 8, 9 Z. Khan, 108 E. A. Khazanov, 135 M. Khursheed, 61 N. Kijbunchoo, 21 Chunglee Kim, 136 J. C. Kim, 137 K. Kim, 92 W. Kim, 55 W. S. Kim, 138 Y.-M. Kim, 139 C. Kimball, 59 E. J. King, 55 P. J. King, 45 M. Kinley-Hanlon, 127 R. Kirchhoff, 8, 9 J. S. Kissel, 45 L. Kleybolte, 140 J. H. Klika, 23 S. Klimenko, 48 T. D. Knowles, 39 P. Koch, 8, 9 S. M. Koehlenbeck, 8, 9 G. Koekoek, 37, 141 S. Koley, 37 V. Kondrashov, 1 A. Kontos, 12 N. Koper, 8, 9 M. Korobko, 140 W. Z. Korth, 1 I. Kowalska, 74 D. B. Kozak, 1 V. Kringel, 8, 9 N. Krishnendu, 29 A. Kr'olak, 142, 143 G. Kuehn, 8, 9 A. Kumar, 130 P. Kumar, 144 R. Kumar, 108 S. Kumar, 16 L. Kuo, 88 A. Kutynia, 142 S. Kwang, 23 B. D. Lackey, 36 K. H. Lai, 92 T. L. Lam, 92 M. Landry, 45 B. B. Lane, 12 R. N. Lang, 145 J. Lange, 58 B. Lantz, 49 R. K. Lanza, 12 A. Lartaux-Vollard, 25 P. D. Lasky, 6 M. Laxen, 7 A. Lazzarini, 1 C. Lazzaro, 52 P. Leaci, 114, 32 S. Leavey, 8, 9 Y. K. Lecoeuche, 45 C. H. Lee, 94 H. K. Lee, 146 H. M. Lee, 147 H. W. Lee, 137 J. Lee, 93 K. Lee, 44 J. Lehmann, 8, 9 A. Lenon, 39 N. Leroy, 25 N. Letendre, 33 Y. Levin, 6, 99 J. Li, 83 K. J. L. Li, 92 T. G. F. Li, 92 X. Li, 46 F. Lin, 6 F. Linde, 37 S. D. Linker, 112 T. B. Littenberg, 148 J. Liu, 63 X. Liu, 23 R. K. L. Lo, 92, 1 N. A. Lockerbie, 24 L. T. London, 68 A. Longo, 149, 150 M. Lorenzini, 14, 15 V. Loriette, 151 M. Lormand, 7 G. Losurdo, 19 J. D. Lough, 8, 9 C. O. Lousto, 58 G. Lovelace, 26 M. E. Lower, 152 H. Luck, 9, 8 D. Lumaca, 30, 31 A. P. Lundgren, 153 R. Lynch, 12 Y. Ma, 46 R. Macas, 68 S. Macfoy, 24 M. MacInnis, 12 D. M. Macleod, 68 A. Macquet, 65 F. Maga˜na-Sandoval, 42 L. Maga˜na Zertuche, 85 R. M. Magee, 87 E. Majorana, 32 I. Maksimovic, 151 A. Malik, 61 N. Man, 65 V. Mandic, 43 V. Mangano, 44 G. L. Mansell, 45, 12 M. Manske, 23, 21 M. Mantovani, 28 M. Mapelli, 51, 52, 96, 154 F. Marchesoni, 50, 41 F. Marion, 33 S. M'arka, 99 Z. M'arka, 99 C. Markakis, 10, 17 A. S. Markosyan, 49 A. Markowitz, 1 E. Maros, 1 A. Marquina, 103 S. Marsat, 36 F. Martelli, 72, 73 I. W. Martin, 44 R. M. Martin, 35 D. V. Martynov, 11 K. Mason, 12 E. Massera, 109 A. Masserot, 33 T. J. Massinger, 1 M. Masso-Reid, 44 S. Mastrogiovanni, 114, 32 A. Matas, 43, 36 F. Matichard, 1, 12 L. Matone, 99 N. Mavalvala, 12 N. Mazumder, 69 J. J. McCann, 63 R. McCarthy, 45 D. E. McClelland, 21 S. McCormick, 7 L. McCuller, 12 S. C. McGuire, 155 J. McIver, 1 D. J. McManus, 21 T. McRae, 21 S. T. McWilliams, 39 D. Meacher, 87 G. D. Meadors, 6 M. Mehmet, 8, 9 A. K. Mehta, 16 J. Meidam, 37 A. Melatos, 98 G. Mendell, 45 R. A. Mercer, 23 L. Mereni, 22 E. L. Merilh, 45 M. Merzougui, 65 S. Meshkov, 1 C. Messenger, 44 C. Messick, 87 R. Metzdorff, 71 P. M. Meyers, 98 H. Miao, 11 C. Michel, 22 H. Middleton, 98 E. E. Mikhailov, 156 L. Milano, 79, 5 A. L. Miller, 48 A. Miller, 114, 32 M. Millhouse, 53 J. C. Mills, 68 M. C. Milovich-Goff, 112 O. Minazzoli, 65, 157 Y. Minenkov, 31 A. Mishkin, 48 C. Mishra, 158 T. Mistry, 109 S. Mitra, 3 V. P. Mitrofanov, 62 G. Mitselmakher, 48 R. Mittleman, 12 G. Mo, 95 D. Moffa, 116 K. Mogushi, 85 S. R. P. Mohapatra, 12 M. Montani, 72, 73 C. J. Moore, 10 D. Moraru, 45 G. Moreno, 45 S. Morisaki, 82 B. Mours, 33 C. M. Mow-Lowry, 11 Arunava Mukherjee, 8, 9 D. Mukherjee, 23 S. Mukherjee, 105 N. Mukund, 3 A. Mullavey, 7 J. Munch, 55 E. A. Mu˜niz, 42 M. Muratore, 34 P. G. Murray, 44 A. Nagar, 86, 159, 160 I. Nardecchia, 30, 31 L. Naticchioni, 114, 32 R. K. Nayak, 161 J. Neilson, 112 G. Nelemans, 64, 37 T. J. N. Nelson, 7 M. Nery, 8, 9 A. Neunzert, 126 K. Y. Ng, 12 S. Ng, 55 P. Nguyen, 70 D. Nichols, 128, 37 S. Nissanke, 128, 37 F. Nocera, 28 C. North, 68 L. K. Nuttall, 153 M. Obergaulinger, 20 J. Oberling, 45 B. D. O'Brien, 48 G. D. O'Dea, 112 G. H. Ogin, 162 J. J. Oh, 138 S. H. Oh, 138 F. Ohme, 8, 9 H. Ohta, 82 M. A. Okada, 13 M. Oliver, 100 P. Oppermann, 8, 9 Richard J. Oram, 7 B. O'Reilly, 7 R. G. Ormiston, 43 L. F. Ortega, 48 R. O'Shaughnessy, 58 S. Ossokine, 36 D. J. Ottaway, 55 H. Overmier, 7 B. J. Owen, 84 A. E. Pace, 87 G. Pagano, 18, 19  \nM. A. Page, 63 A. Pai, 122 S. A. Pai, 61 J. R. Palamos, 70 O. Palashov, 135 C. Palomba, 32 A. Pal-Singh, 140 Huang-Wei Pan, 88 B. Pang, 46 P. T. H. Pang, 92 C. Pankow, 59 F. Pannarale, 114, 32 B. C. Pant, 61 F. Paoletti, 19 A. Paoli, 28 A. Parida, 3 W. Parker, 7, 155 D. Pascucci, 44 A. Pasqualetti, 28 R. Passaquieti, 18, 19 D. Passuello, 19 M. Patil, 143 B. Patricelli, 18, 19 B. L. Pearlstone, 44 C. Pedersen, 68 M. Pedraza, 1 R. Pedurand, 22, 163 A. Pele, 7 S. Penn, 164 C. J. Perez, 45 A. Perreca, 113, 97 H. P. Pfeiffer, 36, 90 M. Phelps, 8, 9 K. S. Phukon, 3 O. J. Piccinni, 114, 32 M. Pichot, 65 F. Piergiovanni, 72, 73 G. Pillant, 28 L. Pinard, 22 M. Pirello, 45 M. Pitkin, 44 R. Poggiani, 18, 19 D. Y. T. Pong, 92 S. Ponrathnam, 3 P. Popolizio, 28 E. K. Porter, 27 J. Powell, 152 A. K. Prajapati, 108 J. Prasad, 3 K. Prasai, 49 R. Prasanna, 130 G. Pratten, 100 T. Prestegard, 23 S. Privitera, 36 G. A. Prodi, 113, 97 L. G. Prokhorov, 62 O. Puncken, 8, 9 M. Punturo, 41 P. Puppo, 32 M. Purrer, 36 H. Qi, 23 V. Quetschke, 105 P. J. Quinonez, 34 E. A. Quintero, 1 R. Quitzow-James, 70 F. J. Raab, 45 H. Radkins, 45 N. Radulescu, 65 P. Raffai, 107 S. Raja, 61 C. Rajan, 61 B. Rajbhandari, 84 M. Rakhmanov, 105 K. E. Ramirez, 105 A. Ramos-Buades, 100 Javed Rana, 3 K. Rao, 59 P. Rapagnani, 114, 32 V. Raymond, 68 M. Razzano, 18, 19 J. Read, 26 T. Regimbau, 33 L. Rei, 60 S. Reid, 24 D. H. Reitze, 1, 48 W. Ren, 17 F. Ricci, 114, 32 C. J. Richardson, 34 J. W. Richardson, 1 P. M. Ricker, 17 K. Riles, 126 M. Rizzo, 59 N. A. Robertson, 1, 44 R. Robie, 44 F. Robinet, 25 A. Rocchi, 31 L. Rolland, 33 J. G. Rollins, 1 V. J. Roma, 70 M. Romanelli, 67 R. Romano, 4, 5 C. L. Romel, 45 J. H. Romie, 7 K. Rose, 116 D. Rosi'nska, 165, 54 S. G. Rosofsky, 17 M. P. Ross, 166 S. Rowan, 44 A. Rudiger, 8, 9, ‡ P. Ruggi, 28 G. Rutins, 167 K. Ryan, 45 S. Sachdev, 1 T. Sadecki, 45 M. Sakellariadou, 132 O. S. Salafia, 168, 169, 170 L. Salconi, 28 M. Saleem, 29 A. Samajdar, 37 L. Sammut, 6 E. J. Sanchez, 1 L. E. Sanchez, 1 N. Sanchis-Gual, 20 V. Sandberg, 45 J. R. Sanders, 42 K. A. Santiago, 35 N. Sarin, 6 B. Sassolas, 22 B. S. Sathyaprakash, 87, 68 P. R. Saulson, 42 O. Sauter, 126 R. L. Savage, 45 P. Schale, 70 M. Scheel, 46 J. Scheuer, 59 P. Schmidt, 64 R. Schnabel, 140 R. M. S. Schofield, 70 A. Schonbeck, 140 E. Schreiber, 8, 9 B. W. Schulte, 8, 9 B. F. Schutz, 68 S. G. Schwalbe, 34 J. Scott, 44 S. M. Scott, 21 E. Seidel, 17 D. Sellers, 7 A. S. Sengupta, 171 N. Sennett, 36 D. Sentenac, 28 V. Sequino, 30, 31, 14 A. Sergeev, 135 Y. Setyawati, 8, 9 D. A. Shaddock, 21 T. Shaffer, 45 M. S. Shahriar, 59 M. B. Shaner, 112 L. Shao, 36 P. Sharma, 61 P. Shawhan, 76 H. Shen, 17 R. Shink, 172 D. H. Shoemaker, 12 D. M. Shoemaker, 77 S. ShyamSundar, 61 K. Siellez, 77 M. Sieniawska, 54 D. Sigg, 45 A. D. Silva, 13 L. P. Singer, 80 N. Singh, 74 A. Singhal, 14, 32 A. M. Sintes, 100 S. Sitmukhambetov, 105 V. Skliris, 68 B. J. J. Slagmolen, 21 T. J. Slaven-Blair, 63 J. R. Smith, 26 R. J. E. Smith, 6 S. Somala, 173 E. J. Son, 138 B. Sorazu, 44 F. Sorrentino, 60 T. Souradeep, 3 E. Sowell, 84 A. P. Spencer, 44 M. Spera, 51, 52, 96, 154, 59 A. K. Srivastava, 108 V. Srivastava, 42 K. Staats, 59 C. Stachie, 65 M. Standke, 8, 9 D. A. Steer, 27 M. Steinke, 8, 9 J. Steinlechner, 140, 44 S. Steinlechner, 140 D. Steinmeyer, 8, 9 S. P. Stevenson, 152 D. Stocks, 49 R. Stone, 105 D. J. Stops, 11 K. A. Strain, 44 G. Stratta, 72, 73 S. E. Strigin, 62 A. Strunk, 45 R. Sturani, 174 A. L. Stuver, 175 V. Sudhir, 12 T. Z. Summerscales, 176 L. Sun, 1 S. Sunil, 108 J. Suresh, 3 P. J. Sutton, 68 B. L. Swinkels, 37 M. J. Szczepa'nczyk, 34 M. Tacca, 37 S. C. Tait, 44 C. Talbot, 6 D. Talukder, 70 D. B. Tanner, 48 M. T'apai, 123 A. Taracchini, 36 J. D. Tasson, 95 R. Taylor, 1 F. Thies, 8, 9 M. Thomas, 7 P. Thomas, 45 S. R. Thondapu, 61 K. A. Thorne, 7 E. Thrane, 6 Shubhanshu Tiwari, 113, 97 Srishti Tiwari, 124 V. Tiwari, 68 K. Toland, 44 M. Tonelli, 18, 19 Z. Tornasi, 44 A. Torres-Forn'e, 177 C. I. Torrie, 1 D. Toyra, 11 F. Travasso, 28, 41 G. Traylor, 7 M. C. Tringali, 74 A. Trovato, 27 L. Trozzo, 178, 19 R. Trudeau, 1 K. W. Tsang, 37 M. Tse, 12 R. Tso, 46 L. Tsukada, 82 D. Tsuna, 82 D. Tuyenbayev, 105 K. Ueno, 82 D. Ugolini, 179 C. S. Unnikrishnan, 124 A. L. Urban, 2 S. A. Usman, 68 H. Vahlbruch, 9 G. Vajente, 1 G. Valdes, 2 N. van Bakel, 37 M. van Beuzekom, 37 J. F. J. van den Brand, 75, 37 C. Van Den Broeck, 37, 180 D. C. Vander-Hyde, 42 L. van der Schaaf, 37 J. V. van Heijningen, 63 A. A. van Veggel, 44 M. Vardaro, 51, 52 V. Varma, 46 S. Vass, 1 M. Vas'uth, 47 A. Vecchio, 11 G. Vedovato, 52 J. Veitch, 44 P. J. Veitch, 55 K. Venkateswara, 166 G. Venugopalan, 1 D. Verkindt, 33 F. Vetrano, 72, 73 A. Vicer'e, 72, 73 A. D. Viets, 23 D. J. Vine, 167 J.-Y. Vinet, 65 S. Vitale, 12 T. Vo, 42 H. Vocca, 40, 41 C. Vorvick, 45 S. P. Vyatchanin, 62 A. R. Wade, 1 L. E. Wade, 116 M. Wade, 116 R. Walet, 37 M. Walker, 26 L. Wallace, 1 S. Walsh, 23 G. Wang, 14, 19 H. Wang, 11 J. Z. Wang, 126 W. H. Wang, 105 Y. F. Wang, 92 R. L. Ward, 21 Z. A. Warden, 34 J. Warner, 45 M. Was, 33 J. Watchi, 101 B. Weaver, 45 L.-W. Wei, 8, 9 M. Weinert, 8, 9 A. J. Weinstein, 1 R. Weiss, 12 F. Wellmann, 8, 9 L. Wen, 63 E. K. Wessel, 17 P. Weßels, 8, 9 J. W. Westhouse, 34 K. Wette, 21 J. T. Whelan, 58 B. F. Whiting, 48 C. Whittle, 12 D. M. Wilken, 8, 9 D. Williams, 44 A. R. Williamson, 128, 37 J. L. Willis, 1 B. Willke, 8, 9 M. H. Wimmer, 8, 9 W. Winkler, 8, 9 C. C. Wipf, 1 H. Wittel, 8, 9 G. Woan, 44 J. Woehler, 8, 9 J. K. Wofford, 58 J. Worden, 45 J. L. Wright, 44 D. S. Wu, 8, 9 D. M. Wysocki, 58 L. Xiao, 1 H. Yamamoto, 1 C. C. Yancey, 76 L. Yang, 115 M. J. Yap, 21 M. Yazback, 48 D. W. Yeeles, 68 Hang Yu, 12 Haocun Yu, 12 S. H. R. Yuen, 92 M. Yvert, 33 A. K. Zadro˙zny, 105, 142 M. Zanolin, 34 T. Zelenova, 28 J.-P. Zendri, 52 M. Zevin, 59 J. Zhang, 63 L. Zhang, 1 T. Zhang, 44 C. Zhao, 63 M. Zhou, 59 Z. Zhou, 59 X. J. Zhu, 6 A. B. Zimmerman, 90 Y. Zlochower, 58 M. E. Zucker, 1, 12 and J. Zweizig 1  \nThe LIGO Scientific Collaboration and the Virgo Collaboration  \n1 LIGO, California Institute of Technology, Pasadena, CA 91125, USA  \n2 Louisiana State University, Baton Rouge, LA 70803, USA  \n3 Inter-University Centre for Astronomy and Astrophysics, Pune 411007, India  \n4 Universit'a di Salerno, Fisciano, I-84084 Salerno, Italy  \n5 INFN, Sezione di Napoli, Complesso Universitario di Monte S.Angelo, I-80126 Napoli, Italy  \n6 OzGrav, School of Physics & Astronomy, Monash University, Clayton 3800, Victoria, Australia  \n7 LIGO Livingston Observatory, Livingston, LA 70754, USA  \n8 Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-30167 Hannover, Germany  \n27  \n16  \n13  \n9 Leibniz Universitat Hannover, D-30167 Hannover, Germany  \n10 University of Cambridge, Cambridge CB2 1TN, United Kingdom  \n11 University of Birmingham, Birmingham B15 2TT, United Kingdom  \n12 LIGO, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \nInstituto Nacional de Pesquisas Espaciais, 12227-010 S˜ao Jos'e dos Campos, S˜ao Paulo, Brazil  \n14 Gran Sasso Science Institute (GSSI), I-67100 L'Aquila, Italy  \n15 INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi, Italy  \nInternational Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bengaluru 560089, India  \n17 NCSA, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA  \n18 Universit'a di Pisa, I-56127 Pisa, Italy  \n19 INFN, Sezione di Pisa, I-56127 Pisa, Italy  \n20 Departamento de Astronom'ıa y Astrof'ısica, Universitat de Val'encia, E-46100 Burjassot, Val'encia, Spain  \n21 OzGrav, Australian National University, Canberra, Australian Capital Territory 0200, Australia  \n22  \nLaboratoire des Mat'eriaux Avanc'es (LMA), CNRS/IN2P3, F-69622 Villeurbanne, France  \n23 University of Wisconsin-Milwaukee, Milwaukee, WI 53201, USA  \n24 SUPA, University of Strathclyde, Glasgow G1 1XQ, United Kingdom  \n25 LAL, Univ. Paris-Sud, CNRS/IN2P3, Universit'e Paris-Saclay, F-91898 Orsay, France  \n26 California State University Fullerton, Fullerton, CA 92831, USA  \nAPC, AstroParticule et Cosmologie, Universit'e Paris Diderot, CNRS/IN2P3, CEA/Irfu, Observatoire de Paris, Sorbonne Paris Cit'e,  \nF-75205 Paris Cedex 13, France  \n28 European Gravitational Observatory (EGO), I-56021 Cascina, Pisa, Italy  \n29 Chennai Mathematical Institute, Chennai 603103, India  \n30 Universit'a di Roma Tor Vergata, I-00133 Roma, Italy  \n31 INFN, Sezione di Roma Tor Vergata, I-00133 Roma, Italy  \n32 INFN, Sezione di Roma, I-00185 Roma, Italy  \n33 Laboratoire d'Annecy de Physique des Particules (LAPP), Univ. Grenoble Alpes, Universit'e Savoie Mont Blanc, CNRS/IN2P3,  \nF-74941 Annecy, France  \n34 Embry-Riddle Aeronautical University, Prescott, AZ 86301, USA  \n35 Montclair State University, Montclair, NJ 07043, USA  \n36 Max Planck Institute for Gravitational Physics (Albert Einstein Institute), D-14476 Potsdam-Golm, Germany  \n37 Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands  \n38 Korea Institute of Science and Technology Information, Daejeon 34141, South Korea  \n39 West Virginia University, Morgantown, WV 26506, USA  \n40 Universit'a di Perugia, I-06123 Perugia, Italy  \n41 INFN, Sezione di Perugia, I-06123 Perugia, Italy  \n42 Syracuse University, Syracuse, NY 13244, USA  \n43 University of Minnesota, Minneapolis, MN 55455, USA  \n44 SUPA, University of Glasgow, Glasgow G12 8QQ, United Kingdom  \n45 LIGO Hanford Observatory, Richland, WA 99352, USA  \n46 Caltech CaRT, Pasadena, CA 91125, USA  \n47 Wigner RCP, RMKI, H-1121 Budapest, Konkoly Thege Mikl'os 'ut 29-33, Hungary  \n48 University of Florida, Gainesville, FL 32611, USA  \n49 Stanford University, Stanford, CA 94305, USA  \n50 Universit'a di Camerino, Dipartimento di Fisica, I-62032 Camerino, Italy  \n51 Universit'a di Padova, Dipartimento di Fisica e Astronomia, I-35131 Padova, Italy  \n52 INFN, Sezione di Padova, I-35131 Padova, Italy  \n53 Montana State University, Bozeman, MT 59717, USA  \n54 Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, 00-716, Warsaw, Poland  \n55 OzGrav, University of Adelaide, Adelaide, South Australia 5005, Australia  \n56 Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universitat Jena, D-07743 Jena, Germany  \n57 INFN, Sezione di Milano Bicocca, Gruppo Collegato di Parma, I-43124 Parma, Italy  \n58  \nRochester Institute of Technology, Rochester, NY 14623, USA  \n59 Center for Interdisciplinary Exploration & Research in Astrophysics (CIERA), Northwestern University, Evanston, IL 60208, USA  \n60 INFN, Sezione di Genova, I-16146 Genova, Italy  \n61 RRCAT, Indore, Madhya Pradesh 452013, India  \n62 Faculty of Physics, Lomonosov Moscow State University, Moscow 119991, Russia  \n63 OzGrav, University of Western Australia, Crawley, Western Australia 6009, Australia  \n64 Department of Astrophysics/IMAPP, Radboud University Nijmegen, P.O. Box 9010, 6500 GL Nijmegen, The Netherlands 65 Artemis, Universit'e Cˆote d'Azur, Observatoire Cˆote d'Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France 66 Physik-Institut, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland 67 Univ Rennes, CNRS, Institut FOTON - UMR6082, F-3500 Rennes, France 68 Cardiff University, Cardiff CF24 3AA, United Kingdom 69 Washington State University, Pullman, WA 99164, USA 70 University of Oregon, Eugene, OR 97403, USA 71 Laboratoire Kastler Brossel, Sorbonne Universit'e, CNRS, ENS-Universit'e PSL, Coll'ege de France, F-75005 Paris, France 72 Universit'a degli Studi di Urbino 'Carlo Bo,' I-61029 Urbino, Italy 73 INFN, Sezione di Firenze, I-50019 Sesto Fiorentino, Firenze, Italy 74 Astronomical Observatory Warsaw University, 00-478 Warsaw, Poland 75 VU University Amsterdam, 1081 HV Amsterdam, The Netherlands 76 University of Maryland, College Park, MD 20742, USA 77 School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA 78 Universit'e Claude Bernard Lyon 1, F-69622 Villeurbanne, France 79 Universit'a di Napoli 'Federico II,' Complesso Universitario di Monte S.Angelo, I-80126 Napoli, Italy 80 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA 81 Dipartimento di Fisica, Universit'a degli Studi di Genova, I-16146 Genova, Italy 82 RESCEU, University of Tokyo, Tokyo, 113-0033, Japan. 83 Tsinghua University, Beijing 100084, China 84 Texas Tech University, Lubbock, TX 79409, USA 85 The University of Mississippi, University, MS 38677, USA 86 Museo Storico della Fisica e Centro Studi e Ricerche 'Enrico Fermi', I-00184 Roma, Italyrico Fermi, I-00184 Roma, Italy 87 The Pennsylvania State University, University Park, PA 16802, USA 88 National Tsing Hua University, Hsinchu City, 30013 Taiwan, Republic of China 89 Charles Sturt University, Wagga Wagga, New South Wales 2678, Australia 90 Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, Ontario M5S 3H8, Canada 91 University of Chicago, Chicago, IL 60637, USA 92 The Chinese University of Hong Kong, Shatin, NT, Hong Kong 93 Seoul National University, Seoul 08826, South Korea 94 Pusan National University, Busan 46241, South Korea 95 Carleton College, Northfield, MN 55057, USA 96 INAF, Osservatorio Astronomico di Padova, I-35122 Padova, Italy 97 INFN, Trento Institute for Fundamental Physics and Applications, I-38123 Povo, Trento, Italy 98 OzGrav, University of Melbourne, Parkville, Victoria 3010, Australia 99 Columbia University, New York, NY 10027, USA 100 Universitat de les Illes Balears, IAC3-IEEC, E-07122 Palma de Mallorca, Spain 101 Universit'e Libre de Bruxelles, Brussels 1050, Belgium 102 Sonoma State University, Rohnert Park, CA 94928, USA 103 Departamento de Matem'aticas, Universitat de Val'encia, E-46100 Burjassot, Val'encia, Spain 104 University of Rhode Island, Kingston, RI 02881, USA 105 The University of Texas Rio Grande Valley, Brownsville, TX 78520, USA 106 Bellevue College, Bellevue, WA 98007, USA 107 MTA-ELTE Astrophysics Research Group, Institute of Physics, Eotvos University, Budapest 1117, Hungary 108 Institute for Plasma Research, Bhat, Gandhinagar 382428, India 109 The University of Sheffield, Sheffield S10 2TN, United Kingdom 110 IGFAE, Campus Sur, Universidade de Santiago de Compostela, 15782 Spain 111 Dipartimento di Scienze Matematiche, Fisiche e Informatiche, Universit'a di Parma, I-43124 Parma, Italy 112 California State University, Los Angeles, 5151 State University Dr, Los Angeles, CA 90032, USA 113 Universit'a di Trento, Dipartimento di Fisica, I-38123 Povo, Trento, Italy 114 Universit'a 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'omic, Universitat de Val'encia, E-46980 Paterna, Val'encia, Spain  \n120 School of Mathematics, University of Edinburgh, Edinburgh EH9 3FD, United Kingdom  \n121 Institute Of Advanced Research, Gandhinagar 382426, India  \n122 Indian Institute of Technology Bombay, Powai, Mumbai 400 076, India  \n123 University of Szeged, D'om t'er 9, Szeged 6720, Hungary  \n124 Tata Institute of Fundamental Research, Mumbai 400005, India  \n125 INAF, Osservatorio Astronomico di Capodimonte, I-80131, Napoli, Italy 126 University of Michigan, Ann Arbor, MI 48109, USA 127 American University, Washington, D.C. 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 Biaglyph[suppress]lystok, 15-424 Biaglyph[suppress]lystok, Poland 132 King's College London, University of London, London WC2R 2LS, United Kingdom 133 University of Southampton, Southampton SO17 1BJ, United Kingdom 134 University of Washington Bothell, Bothell, WA 98011, USA 135 Institute of Applied Physics, Nizhny Novgorod, 603950, Russia 136 Ewha Womans University, Seoul 03760, South Korea 137 Inje University Gimhae, South Gyeongsang 50834, South Korea 138 National Institute for Mathematical Sciences, Daejeon 34047, South Korea 139 Ulsan National Institute of Science and Technology, Ulsan 44919, South Korea 140 Universitat Hamburg, D-22761 Hamburg, Germany 141 Maastricht University, P.O. Box 616, 6200 MD Maastricht, The Netherlands 142 NCBJ, 05-400 ' Swierk-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, South Korea 147 Korea Astronomy and Space Science Institute, Daejeon 34055, South Korea 148 NASA Marshall Space Flight Center, Huntsville, AL 35811, USA 149 Dipartimento di Matematica e Fisica, Universit'a 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, United Kingdom 154 Institut fur Astro- und Teilchenphysik, Universitat 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'e 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'ora, 65-265 Zielona G'ora, Poland 166 University of Washington, Seattle, WA 98195, USA 167 SUPA, University of the West of Scotland, Paisley PA1 2BE, United Kingdom 168 NAF, Osservatorio Astronomico di Brera sede di Merate, I-23807 Merate, Lecco, Italy 169 Universit'a 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'e de Montr'eal/Polytechnique, Montreal, Quebec 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  \nAndrews University, Berrien Springs, MI 49104, USA  \n177 Max Planck Institute for Gravitationalphysik (Albert Einstein Institute), D-14476 Potsdam-Golm, Germany 178 Universit'a di Siena, I-53100 Siena, Italy  \n179 Trinity University, San Antonio, TX 78212, USA  \n180 Van Swinderen Institute for Particle Physics and Gravity, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands  \n(Received 2018 December 15; Revised 2019 July 11; Accepted 2019 July 21; Published 2019 September 9)", 'ABSTRACT': 'We present results on the mass, spin, and redshift distributions with phenomenological population models using the ten binary black hole mergers detected in the first and second observing runs completed by Advanced LIGO and Advanced Virgo. We constrain properties of the binary black hole (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 black hole in such binaries is well approximated by models with no more than 1% of black holes more massive than 45 M glyph[circledot] , and a power law index of α = 1 . 3 +1 . 4 -1 . 7 (90% credibility). We also show that BBHs are unlikely to be composed of black holes with large spins aligned to the orbital angular momentum. Modelling 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 +55 . 8 -28 . 2 Gpc -3 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 black holes via supernovae, binary interactions of massive stars, stellar cluster dynamics, and the formation history of black holes across cosmic time.', '1. INTRODUCTION': "The second LIGO/Virgo observing run (O2) spanned nine months between November 2016 through August 2017, building upon the first, four-month run (O1) in 2015. The LIGO/Virgo gravitational-wave (GW) interferometer network is comprised 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, ten binary black hole (BBH) mergers have been detected to date (Abbott et al. 2018). 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 +3 . 3 -0 . 7 M glyph[circledot] . GW170729 (Abbott et al. 2018)exceptional in several ways-is likely to be the heaviest BBH to date, having total mass 85 . 2 +15 . 4 -11 . 2 M glyph[circledot] , as well as the most distant, at redshift 0 . 48 +0 . 19 -0 . 20 . Both GW151226 and GW170729 show evidence for at least one black hole with a spin greater than zero (Abbott et al. 2016b; Abbott et al. 2018).  \nBy 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  \nBayesian inference and modelling (Gelman et al. 2004; Mandel 2010; Foreman-Mackey et al. 2014; Hilbe et al. 2017; ase 2018) which, when applied to parameterized models of the population, is able to infer populationlevel 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 (Mandel 2010; Hogg et al. 2010).  \nThe structure and parameterization of BBH populations models are guided by the physical processes and evolutionary environments in which BBH 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 commonenvelope evolution (Bethe & Brown 1998; Portegies Zwart & Yungelson 1998; Belczynski et al. 2002; Voss & Tauris 2003; Dewi et al. 2006; Belczynski et al. 2007, 2008; Dominik et al. 2013; Belczynski et al. 2014; Mennekens & Vanbeveren 2014; Spera et al. 2015; Tauris et al. 2017; Eldridge & Stanway 2016; Stevenson et al. 2017b; Chruslinska et al. 2018; Mapelli et al. 2017; Giacobbo et al. 2018; Mapelli & Giacobbo 2018; Kruckow et al. 2018; Giacobbo & Mapelli 2018) or via chemically homogeneous evolution (Marchant et al. 2016; de Mink &Mandel 2016; Mandel & de Mink 2016). Alternatively, BBHs might form via dynamical processes in stellar clus-  \nters (Portegies Zwart & McMillan 2000; Kulkarni et al. 1993; Sigurdsson & Hernquist 1993; Grindlay et al. 2006; O'Leary et al. 2006; Sadowski et al. 2008; Ivanova 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; Kimpson et al. 2016; Antonini et al. 2017; Liu & Lai 2018), gas drag and stellar scattering in accretion disks surrounding super-massive black holes (McKernan et al. 2012; Bartos et al. 2017; Stone et al. 2017). Finally, BBHs might originate as part of a primordial black hole population in the early Universe (Carr & Hawking 1974; Carr et al. 2016; Sasaki et al. 2016; Inomata et al. 2017; Inayoshi et al. 2016; Bird et al. 2016; AliHa¨ımoud et al. 2017; Clesse & Garc'ıa-Bellido 2017; Chen &Huang 2018; Ando et al. 2018), where their mass spectrum is typically proposed as having power law behavior, but spanning a much wider range of masses than stellar mass BH. Each channel contributes differently to the distributions of the mass, spin, distance, and orbital characteristics of BBHs.  \nThere are several processes common to most pathways through stellar evolution which affect the properties of the resultant BBH system. Examples include mass loss (Vink et al. 2001; Vink & de Koter 2005; Grafener & 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 black holes 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 which are not disrupted.  \nTheoretical and phenomenological models of BBH formation are explored by population synthesis. This requires modelling not only of 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  \neffects 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.  \nElectromagnetic observations and modeling of systems containing black holes have led to speculation about the existence of potential gaps in the black hole mass spectrum. Both gaps may be probed using data from current ground-based gravitational-wave interferometers, and as such, have been the target of parametric studies. At low masses, observations of X-ray binaries (XRB) combined via Bayesian population modeling (Bailyn et al. 1998; Ozel et al. 2010; Farr et al. 2011b) suggest a minimum black hole mass well above the largest neutron star masses. While the existence and nature of this gap is still uncertain (Kreidberg et al. 2012), it is proposed to exist between the most massive neutron stars ( Ozel & Freire 2016; Freire et al. 2008; Margalit & Metzger 2017) (2 . 1 -2 . 5 M glyph[circledot] ) and the lightest black holes ∼ 5 M glyph[circledot] . It is possible to constrain the existence of this lower mass gap with GW observations (Littenberg et al. 2015; Mandel et al. 2015; Kovetz et al. 2017; Mandel et al. 2017). In Section 3, we find our current GW observations do not inform the upper edge of this gap, inferring a minimum mass on the primary black hole at m min glyph[lessorsimilar] 9 M glyph[circledot] . Our volumetric sensitivity to BBH systems with masses less than 5 M glyph[circledot] 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.  \nRecently, there have been claims of an upper cutoff in the BBH mass spectrum based on the first few LIGO detections (Fishbach & Holz 2017; Talbot & Thrane 2018; Wysocki et al. 2018; Bai et al. 2018; Roulet & Zaldarriaga 2019). This might be expected as a consequence of a different supernova type, called the (pulsational) pairinstability supernova (Heger & Woosley 2002; Belczynski et al. 2016b; Woosley 2017; Spera & Mapelli 2017; Marchant et al. 2018). Evolved stars with a Helium core mass glyph[greaterorsimilar] 30 M glyph[circledot] are expected to become unstable, because efficient pair production softens their equation of state. For Helium core mass ∼ 30 - 64 M glyph[circledot] , 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  \nexpected to affect the final collapse of the star, leading to smaller black hole masses. The fate of a star with He core mass ∼ 64 - 135 M glyph[circledot] 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 black hole remnants between ∼ 50 - 150 M glyph[circledot] 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; Gerosa & Berti 2017; Fishbach et al. 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 black holes above ∼ 45 M glyph[circledot] .  \nObservational constraints on the BBH merger rate (Abbott et al. 2016c; Abbott et al. 2018) generally assume a rate density which 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. 2016a). For instance, theoretical models of isolated evolution through common envelope lead to a distribution of times to merger p ( t GW ) ∝ t -1 GW (Dominik et al. 2012; Belczynski et al. 2016a). 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 ∼ 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.  \nBlack hole spin measurements also provide a powerful tool to discriminate between different channels of BBHformation (Mandel & O'Shaughnessy 2010; Abbott et al. 2016d; Rodriguez et al. 2016c; Vitale et al. 2017; Gerosa & Berti 2017; Farr et al. 2017, 2018; Gerosa et al. 2018). For example, BBHs formed in a dynamic environ-  \nment 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 black holes and the orbital angular momentum of the system (Tutukov & Yungelson 1993; Kalogera 2000; Grandcl'ement et al. 2004; Belczynski et al. 2016a; Rodriguez et al. 2016c; Mandel & de Mink 2016; Marchant et al. 2016; Stevenson et al. 2017b; O'Shaughnessy et al. 2017; 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 which decline as the spin magnitude increases from zero, but our ability to distinguish between assumed distributions of spin orientation is very limited.  \nGW170817, 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 August of 2017. A subsequent post-merger transient (AT 2017gfo) was observed across the electromagnetic spectrum, from radio (Alexander et al. 2017), NIR/optical (Coulter et al. 2017; Soares-Santos et al. 2017; Chornock et al. 2017; Cowperthwaite et al. 2017; Nicholl et al. 2017; Pian et al. 2017), to X-ray (Troja et al. 2017; Margutti et al. 2017) and γ -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 gravitational-wave measurement is mostly compatible with the observed Galactic population ( Ozel et al. 2012). However, if GW170817 did form a black hole, it would also occupy the lower mass gap described previously.  \nWe structure the paper as follows. First, notation and models are established in Section 2. Section 3 describes our modeling of the black hole mass distribution, followed by rate distributions and evolution in Section 4. The black hole spin magnitude and orientation distributions are discussed in Section 5. We conclude in Section 6. Studies of various systematics are presented  \nin 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 1 .", '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 run (O1 and O2) (Abbott et al. 2018). We do not include marginal detections, but these likely have a minimal impact 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.  \nThe individual properties of those 10 sources were inferred using a Bayesian framework, with results summarized in Abbott et al. (2018). 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 ; a discussion of potential systematic biases in our inference are discussed in Appendix A. We also refer to Appendix B in (Abbott et al. 2018) for more details on comparisons between those two waveform families.  \nTo 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 black hole 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. 2018). 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  \nof 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: right ascension, declination, 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; Klein et al. 2018; Hinder 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.  \nIn this work, we define the mass ratio as q = m 2 /m 1 where m 1 ≥ m 2 . The frequency of gravitational wave emission is directly related to the component masses. However, due to the expansion of spacetime as the gravitational wave 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 sourceframe 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).  \nWe characterize black hole spins using the dimensionless spin parameter χ i = S i /m 2 i . Of particular interest are the magnitude of the dimensionless spin, a i = | χ i | , and the tilt angle with respect to the orbital angular momentum, ˆ L , given by cos t i = ˆ L · ˆ χ i . We also define an overall effective spin, χ eff (Damour 2001; Racine 2008; Ajith et al. 2011), which is a combination of the individual spin components along to orbital angular mo-  \nχ eff = ( χ 1 + q χ 2 ) · ˆ L 1 + q . (1)  \nχ eff is approximately proportional to the lowest order contribution to the GW waveform phase that contains spin for systems with similar masses. Additionally, χ eff is conserved throughout 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) which 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); 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.  \nThe general model family has 8 parameters to characterize the mass model; 3 to characterize each black hole's spin distribution; one parameter describing the  \nlocal merger rate, R 0 ; and one parameter characterizing redshift dependence. We refer to the set of these population parameters as θ . All of the population parameters introduced in this section are summarised 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 black hole mass. With our sample of ten binaries, we extend this analysis by considering three increasingly complex models for the distribution of black hole masses. The first extension, Model A (derived from Fishbach & Holz (2017); Wysocki et al. (2018)), allows the maximum black hole mass m max and the power-law index α to vary. In Model B (derived from Kovetz et al. (2017); Fishbach & Holz (2017); Talbot & Thrane (2018)) the minimum black hole mass m min and the mass ratio power-law index β q are also free parameters. However, the priors on Model B and C enforce a minimum of 5 M glyph[circledot] on m min - see Table 2. Explicitly, the mass distribution in Model A and Model B takes the form  \np ( m 1 , m 2 | m min , m max , α, β q ) ∝    C ( m 1 ) m -α 1 q β q if m min ≤ m 2 ≤ m 1 ≤ m max 0 otherwise , (2)  \nwhere C ( m 1 ) is chosen so that the marginal distribution is a power law in m 1 : p ( m 1 | m min , m max , α, β q ) = m -α 1 .  \nModel A fixes m min = 5 M glyph[circledot] and β 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 ) ∝ q β q . When β q = 0, C ( m 1 ) ∝ 1 / ( m 1 -m min ), as assumed in Abbott et al. (2016c, 2017e).  \nModel C (from Talbot & Thrane (2018)) further builds upon the mass distribution in Equation 2 by allowing for a second, Gaussian component at high mass, as well as introducing smoothing scales δm , which taper the hard edges of the low- and high-mass cutoffs of the primary  \nand secondary mass power-law. The second Gaussian component is designed to capture a possible build-up of high-mass black holes 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, µ m , and standard deviation, σ m , of the Gaussian component, λ m , the fraction of primary black holes in this Gaussian component, and δm the smoothing scale at the low mass end of the distribution.  \nThe full form of this distribution is  \np ( m 1 | θ ) = [ (1 -λ m ) A ( θ ) m -α 1 Θ( m max -m 1 ) + λ m B ( θ ) exp ( -( m 1 -µ m ) 2 2 σ 2 m )] S ( m 1 , m min , δm ) , p ( q | m 1 , θ ) = C ( m 1 , θ ) q β q S ( m 2 , m min , δm ) . (3)  \nThe factors A , B , and C ensure each of the power-law component, Gaussian component, and mass ratio distri-  \nbutions are correctly normalized. S is a smoothing func-  \nTable 1. Parameters describing the binary black hole population. See the text for a more thorough discussion and the functional forms of the models.  \ntion which rises from zero at m min to one at m min + δm as defined in Talbot & Thrane (2018). Θ 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 black hole 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 both black hole spin magnitudes in a binary, a i , are drawn from a beta distribution (Wysocki et al. 2018):  \np ( a i | α a , β a ) = a α a -1 i (1 -a i ) β a -1 B( α a , β a ) . (4)  \nThis 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 upper right hand panel of Figure 1. We choose to model the moments of the beta distribution using the mean ( E [ a ]) and variance (Var[ a ]), given by  \nE [ a ] = α a α a + β a ; Var[ a ] = α a β a ( α a + β a ) 2 ( α a + β a +1) . (5)  \nWe adopt a prior on the spin magnitude model parameters which are uniform over the values of E [ a ] and Var[ a ] which satisfy α a , β a ≥ 1. This choice of which values to sample avoids numerically-challenging singular spin distributions.  \nTo describe the spin orientation, we assume that the tilt angles between each black hole spin and the orbital  \nangular 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)  \np (cos t 1 , cos t 2 | σ 1 , σ 2 , ζ ) = (1 -ζ ) 4 + 2 ζ π ∏ i ∈{ 1 , 2 } exp ( -(1 -cos t i ) 2 / (2 σ 2 i ) ) σ i erf( √ 2 /σ i ) . (6)  \nWe 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 lower right panel of Figure 1.  \nThe parameter ζ denotes the fraction of binaries which are preferentially aligned with the orbital angular momentum; ζ = 1 implies all black hole spins are preferentially aligned and ζ = 0 is an isotropic distribution of spin orientations. The typical degree of spin misalignment is represented by the σ i . For spin orientations we explore two parameterized families of models:  \n- · Gaussian (G): ζ = 1.\n- · Mixture (M): 0 ≤ ζ ≤ 1.  \nThe 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.  \nFigure 1. The probability distributions for models encoded by Equations 3, 4, and 6, are shown in the left panels, upper right panel, and lower 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 glyph[circledot] rather than marginalizing for simplicity.  \n<!-- image -->', '2.5. Redshift Evolution Models': "The previous two subsections described the probability distributions of intrinsic parameters p ( ξ ) (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 glyph[lessorsimilar] 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  \nR ( z | λ ) = R 0 (1 + z ) λ , (7)  \nwhere R 0 is the rate density at z = 0. In this model, λ = 0 corresponds to a merger rate density that is uniform in comoving volume and source-frame time, while λ ∼ 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 ( λ < 0) to rate densities that peak earlier than the star-formation rate ( λ glyph[greaterorsimilar] 3).  \nThese 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.  \nThe general model family, including the distributions of masses, spins and merger redshift, is therefore given by the distribution  \nd N d ξ d z ( θ ) = R ( z ) [ d V c d z ( z ) ] T obs 1 + z p ( ξ | θ ) , (8)  \nwhere 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 obs , as measured in the detector-frame, θ is the collection of all hyper-parameters 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  \nR ( z ) = d N d V c d t ( z ) , (9)  \nwhere t is the time in the source-frame, so that Eq. 8 can be written equivalently in terms of the merger rate density:  \nd R d ξ ( z | θ ) = R 0 p ( ξ | θ )(1 + z ) λ . (10)  \n2.6. Hierarchical Population Model  \nWe perform a hierarchical Bayesian analysis, accounting for measurement uncertainty and selection effects (Loredo 2004; Abbott et al. 2016c; Wysocki et al. 2018; Fishbach et al. 2018; Mandel et al. 2018; 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 2 . The likelihood of the observed GW data given the population hyperparameters θ that describe the general astrophysical distribution, dN/dξdz , is given by the inhomogeneous Poisson likelihood:  \nL ( { d n }| θ ) ∝ e -µ ( θ ) N obs ∏ n =1 ∫ L ( d n | ξ, z ) d N d ξ d z ( θ ) d ξ d z, (11)  \nwhere µ ( θ ) is the rate constant describing the mean number of events as a function of the population hyperparameters, N obs is the number of detections, and L ( d n | ξ, z ) is the individual-event likelihood for the n th detection having parameters ξ, z .  \nIn order to calculate the expected number of detections µ ( θ ), 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 V T ( ξ ) of a network with a given sensitivity to be  \nV T ( ξ ) = T obs ∫ ∞ 0 f ( z | ξ ) d V c d z 1 1 + z d z, (12)  \nwhere the sensitivity is assumed to be constant over the observing time, T obs , as measured in the detector-frame and f ( z | ξ ) is the detection probability of a BBH with the given parameter set ξ 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 hyper-parameters θ , we can calculate the total observed spacetime volume  \n〈 V T 〉 θ = ∫ ξ p ( ξ | θ ) V T ( ξ )d ξ, (13)  \nwhere p ( ξ | θ ) 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 glyph[circledot] . We then use a semi-analytic prescription, calibrated to the simulation results, to derive the 〈 V T 〉 θ for specific hyper-parameters.  \nAllowing the merger rate to evolve with redshift, the expected number of detections is given by  \nµ ( θ ) = T obs ∫ ξ ∫ ∞ 0 p ( ξ | θ ) f ( z | ξ ) R ( z ) d V c d z 1 1 + z d z d ξ. (14)  \nIf the merger rate does not evolve with redshift, i.e., R ( z ) = R 0 , this reduces to µ ( θ ) = R 0 〈 V T 〉 θ .  \nWe note that the hyperparameter likelihood given by Eq. 11 reduces to the likelihood used in the O1 mass distribution analysis (Eq. 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 ) ∝ 1 /R 0 (Fishbach et al. 2018; Mandel et al. 2018). For consistency  \nwith previous analyses, we adopt a flat-in-log prior on the rate parameter throughout this work.", '2.7. Statistical Framework and Prior Choices': "In practice, we sample the likelihood L ( d n | ξ, 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 Eq. 11 (Hogg et al. 2010; Mandel 2010) (see Appendix C).  \nWhere 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  \nadopted in Abbott et al. (2018), we note that similar results are obtained on the rates (see Sec. 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 4 effective samples, after thinning by their autocorrelation time.  \nThe normalization factor of the posterior density in Bayes' 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.  \nWe 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., ξ , 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:  \np ( ξ new | ξ observed ) = ∫ p ( ξ new | θ ) p ( θ | ξ observed ) dθ (15)  \nIt is a predictor for future merger values ξ new given observed data ξ 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. (2018) 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 - 500 M glyph[circledot] , 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 glyph[circledot] ) black holes, a binary which is bounded away from equal mass in its posterior, or a binary with a component mass confidently below 5 M glyph[circledot] . In our survey, we also find a preponderance of observations at higher masses: six with significant posterior support above 30 M glyph[circledot] . In this section, we attempt to reconstruct the binary black hole  \nmerger 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 rate 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.  \nIf we assume the black hole masses follow a powerlaw distributed and fix the minimum black hole mass to be m min = 5 M glyph[circledot] (Model A), we find α = 0 . 4 +1 . 4 -1 . 9 and m max = 41 . 6 +9 . 6 -4 . 3 M glyph[circledot] . In Model B we infer the powerlaw index of the primary mass to be α = 1 . 3 +1 . 4 -1 . 7 with corresponding limits m min = 7 . 8 +1 . 2 -2 . 5 M glyph[circledot] and m max = 40 . 8 +11 . 8 -4 . 4 M glyph[circledot] .  \nFigure 4, shows the posterior over the population parameters present in A and B, as well as a second, Gaussian population parameterized with m max and σ m . λ m is the mixing fraction of binaries in the Gaussian population versus the power law, with λ m = 0 indicating only the power law component. The Gaussian component is centered at µ m = 29 . 8 +5 . 8 -7 . 3 M glyph[circledot] , has a width σ m = 6 . 4 +3 . 2 -4 . 2 M glyph[circledot] , 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 α = 7 . 1 +4 . 4 -4 . 8 than Models A or B, however, the posterior distribution for Model C is less informative for α glyph[greaterorsimilar] 4. This in turn means that we cannot constrain the parameter m max in Model C since the powerlaw component has negligible support above ∼ 45 M glyph[circledot] (see the upper panel of Figure 2). In the intermediate regime, ∼ 15 M glyph[circledot] -25 M glyph[circledot] , Model C infers a smaller rate than Models A or 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 black hole mass, in this model we find m min = 6 . 9 +1 . 7 -2 . 8 M glyph[circledot] . All three models produce consis-  \nTable 2. Summary of models used in Sections 3, 4, and 5, with the prior ranges for the population parameters. The fixed parameters are in bold. Each of these distributions is uniform over the stated range. All models in this Section assume rates which are uniform in the comoving volume ( λ = 0). The lower limit on m min is chosen to be consistent with Abbott et al. (2018).  \nFigure 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), and 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.  \n<!-- image -->  \ns for the marginal merger rate distribution, as is further discussed in Section 4.  \nAll 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  \n(see Equation 15). Model A obtains 44 . 0 M glyph[circledot] , Model B obtains 41 . 8 M glyph[circledot] , and Model C obtains 41 . 8 M glyph[circledot] . Therefore, all models self-consistently infer a dearth of black holes above ∼ 45 M glyph[circledot] . This is determined by the lower limit for the mass of the most massive black hole in the sample because m max can be no smaller than this value. Similarly, the models which allow m min to vary (B and C) disfavor populations with m min above glyph[similarequal] 9 M glyph[circledot] . This parameter is close to the largest allowed mass for the  \nFigure 3. One- and two-dimensional posterior distributions for the hyperparameters describing Models A and B. Large values of α correspond to a mass distribution which rapidly decays with increasing mass. Large values of β correspond to a mass-ratio distribution which prefers equal mass binaries. Also shown is the one-dimensional posterior distribution for the merger rate discussed in Abbott et al. (2018), and the stability of Model A to the removal of the GW170729 event.  \n<!-- image -->  \nleast massive black hole in the sample, for similar reasons.  \nThe lower limits we place on m min are dominated by our prior choices that constrain m min ∈ [5 , 10] M glyph[circledot] (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 glyph[circledot] . Given current sensitivities, this is to be expected (Littenberg et al. 2015; Mandel et al. 2015). In the inspiral-dominated regime, the sensitive time-  \nlume scales as V T ∼ m 15 / 6 (Finn & Chernoff 1993); extending our inferred mass distributions and merger rates into the possible lower black hole mass gap from 3-5 M glyph[circledot] ( Ozel et al. 2010; Farr et al. 2011b; Kreidberg et al. 2012) yields an expected number of detected BBH mergers glyph[lessorsimilar] 1. Thus, we are unable to place meaningful constraints on the presence or absence of a mass gap at low black hole mass.  \nFigure 4. One- and two-dimensional 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 α , β , m max , and m min describe the power-law component. The Gaussian has mean µ m and standard deviation σ m . The fraction of black holes in the Gaussian component is λ m . This model also allows for a gradual turn-on at low masses over a mass range δ m .  \n<!-- image -->  \nModels B and C also allow the distribution of mass ratios to vary according to β q . In these cases the inferred mass-ratio distribution favors comparable-mass binaries (i.e., distributions with most support near q glyph[similarequal] 1), see panel two of Figure 2. Within the context of our parameterization, we find β q = 6 . 9 +4 . 6 -5 . 7 for Model B and β q = 4 . 5 +6 . 6 -5 . 2 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 β q returns the prior for β q glyph[greaterorsimilar] 4. Thus, we cannot say much about the relative likelihood of asymmetric binaries, beyond their overall rarity.  \nThe distribution of the parameter controlling the fraction of the power law versus the Gaussian component in Model C is λ m = 0 . 3 +0 . 4 -0 . 2 , 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 amongst the three models presented in this Section, we compute the Bayes factors comparing the mass models using a nested sampler (Skilling 2004), CPNest (Veitch et al. 2017). These are shown in Table 3. Model B, which allows m min and β q to vary is preferred over Model A (ln BF A B = -1 . 42). To isolate the contributions of the Gaussian component and low mass smoothing in Model C, we compute the SavageDickey density ratio, p ( θ = 0) /p prior ( θ = 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 λ =0 C = -1 . 92); nevertheless, all models infer mass distributions 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 or sharp cutoff at low mass (ln BF δ m =0 C = 0 . 14). This is unsurprising, since we are less sensitive to structure in the mass distribution at low masses (Talbot & Thrane 2018).  \nThe analysis above includes all ten binary black hole 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 repeated these analyses after omitting the least significant detection. For  \nTable 3. The log Bayes factor comparing each of the models described in Table 2 to the most complex model, Model C. The evidence for the three mass models is computed using nested sampling, while the limits λ m = 0 and δ m = 0 of Model C are computed using the Savage-Dickey density ratio.  \nour 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 black hole mass m max reported in models A and B is decreased by about 5 M glyph[circledot] . Without GW170729, the m max distribution is 38 . 3 +7 . 3 -3 . 6 M glyph[circledot] for Model A and 37 . 3 +8 . 5 -3 . 4 M glyph[circledot] 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 the power law to be α = 2 . 3 +1 . 3 -1 . 4 at 90% confidence assuming a minimum black hole mass of 5 M glyph[circledot] and maximum total mass of 100 M glyph[circledot] . None of our models directly emulate this one, but Model A is the closest analog. When allowing m max to vary, 100 M glyph[circledot] 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.  \nIn Figure 5, we highlight the two mass gaps predicted by models of stellar evolution: the first gap between ∼ 2 and ∼ 5 M glyph[circledot] and the second between ∼ 50 M glyph[circledot] and ∼ 150 M glyph[circledot] , compared against the observed black holes. 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.  \nThe minimum mass of a black hole and the existence of a mass gap between neutron stars and black holes (lower gray shaded area, right panel of Figure 5) are currently debated. Claims ( Ozel et al. 2010; Farr et al. 2011b) of the existence of a mass gap between the heaviest neu-tr  \nFigure 5. The left-hand panel shows 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 black holes 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. In the right-hand panel, the predicted compact-object mass is shown 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 × 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-hand panel to emphasize that BNS mergers might fill the lower gap.  \n<!-- image -->  \nn stars ( ∼ 2 M glyph[circledot] ) and the lightest black holes ( ∼ 5 M glyph[circledot] ) are based on the sample of about a dozen X-ray binaries with dynamical mass measurements. However, Kreidberg et al. (2012) suggested that the dearth of observed black hole 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 modelling. We can see in Figure 5 that none of the observed binaries sit in this gap, but the sample is not sufficient to definitively confirm or refute the existence of this mass gap.  \nFrom the first six announced BBH detections, Fishbach & Holz (2017) argued that there is evidence for missing black holes with mass greater than glyph[greaterorsimilar] 40 M glyph[circledot] . The existence of this second mass gap - see the upper grey shaded area in the right panel of Figure 5 between ∼ 50 M glyph[circledot] and ∼ 150 M glyph[circledot] - has been further explored by Talbot & Thrane (2018); Wysocki et al. (2018); Bai  \net al. (2018); 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 build-up. Predictions for the maximum mass of black holes born after pulsational pair-instability supernovae are ∼ 50 M glyph[circledot] (Belczynski et al. 2016b; 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; Abbott et al. 2018; Fishbach & Holz 2017; Wysocki et al. 2018; Fishbach et al. 2018), the inferred binary black hole 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. 2018), we infer the overall BBH merger rate for two fixed-parameter populations. The first of these populations follows the powerlaw model given by Equation 2 with α = 2 . 3, β q = 0, m min = 5 M glyph[circledot] , and m max = 50 M glyph[circledot] . The second population follows a distribution in which both black hole masses are independently drawn from a flat-in-log distribution:  \np ( m 1 , m 2 ) ∝ 1 m 1 m 2 , (16)  \nsubject to the same mass cutoffs 5 M glyph[circledot] < m 2 < m 1 < 50 M glyph[circledot] as the fixed power-law population. Both the power-law and flat-in-log populations assume an isotropic and uniform-magnitude spin distribution ( α a = β a = 1). These two fixed-parameter populations are used to estimate the population-averaged sensitive volume 〈 V T 〉 with a Monte-Carlo injection campaign as described in Abbott et al. (2018), with each population corresponding to a different 〈 V T 〉 because of the strong correlation between the mass spectrum and the sensitive volume. Under the assumption of a constantin-redshift rate density, these 〈 V T 〉 estimates yield two different estimates of the rate: 57 +40 -25 Gpc -3 yr -1 for the α = 2 . 3 population, and 19 +13 -8 . 2 Gpc -3 yr -1 for the flat-in-log population (90% credibility; combining the rate posteriors from the two analysis pipelines).  \nThe 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 V T ( ξ ) for each point in parameter space with the full Monte-Carlo injection campaign described in Abbott et al. (2018), so we employ the semi-analytic methods described in Appendix A. Furthermore, while the rate calculations in Abbott et al. (2018) 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. 2015; Abbott et al. 2016e), in this work we fix a high detection threshold Abbott et al. (2018), which sets the number of detections to N obs = 10. In principle, our results  \nare 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 are enumerated in Table 4.  \nIn these calculations, we first maintain the assumption in Abbott et al. (2018) 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. Non-Evolving 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 selfconsistently determining the rate using the models for the mass and spin distribution described in the previous two sections.  \nTable 5 contains the intervals on the distribution of R 0 for all three models. For Models B and C we deduce a merger rate between 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 glyph[circledot] , whereas Models B and C favor a higher minimum mass and therefore larger populationaveraged 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. (2018). However, the fixedparameter models in Abbott et al. (2018) are disfavored by our full fit to the mass distribution, particularly with respect to the maximum mass. Our results favor maximum masses glyph[lessorsimilar] 45 M glyph[circledot] , rather than 50 M glyph[circledot] as used in Abbott et al. (2018), and power-law slopes closer to α ∼ 1. For this reason, although we infer a mass distribution slope that is similar to the flat-in-log population from Abbott et al. (2018), we infer a rate that is closer to the rate inferred for the fixed-parameter power-law model 3 . While 〈 V T 〉 gets larger (implying a smaller rate estimate) as α is decreased, decreasing m max has the opposite effect, and so the 〈 V T 〉 for the fixed-parameter power-law model is similar to the 〈 V T 〉 s for our best-fit  \na This model assumes the black holes have zero spin.  \nTable 4. Summary of models in Section 4, with prior ranges for the population parameters determining the rate models. The fixed parameter models are drawn from Abbott et al. (2018). 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.Table 5. This table lists the BBH merger rate intervals for each of the mass models tested. These rates assume no evolution in redshift, but otherwise marginalize over all other population parameters.  \nmass distributions, which favor smaller α and smaller m max .  \nWe note that while our analysis differs from the rate calculations in Abbott et al. (2018) by the choice of prior on the rate parameter (log-uniform in this work compared to a Jeffreys prior p ( R 0 ) ∝ R -0 . 5 0 in Abbott et al. (2018)), 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 +74 . 0 -33 . 3 Gpc -3 yr -1 , whereas under a Jeffreys prior this shifts by only ∼ 10% to 56 . 7 +65 . 4 -30 . 4 Gpc -3 yr -1 .', '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 including the star-formation rate, time-delay distribution, metallicity evolution, and globular cluster formation rate (Dominik et al. 2013; Belczynski et al. 2016a; Mandel & de Mink 2016; Rodriguez & Loeb 2018). Therefore, in this section, we allow the merger rate to evolve with redshift, and 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  \nthat the underlying mass distribution does not vary with redshift. We therefore fit the joint mass-redshift distribution according to the model:  \ndR dm 1 dm 2 ( z ) = R 0 p ( m 1 , m 2 | α, m max )(1 + z ) λ (17)  \nNote 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 ∼ 2, which is currently far beyond the horizon redshift for BBH detections.  \nFigure 6 shows the merger rate density as a function of redshift (blue band), compared to the rate inferred in subsection 4.1 for the non-evolving Model B (orange band). The joint posterior PDF on λ , α , and m max , marginalized over the local rate parameter R 0 , is shown  \nFigure 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. In orange, the solid line and shaded region shows the median and 90% credible interval of the rate inferred for Model B as discussed in subsection 4.1, assuming a non-evolving merger rate.  \n<!-- image -->  \nFigure 7. The posterior PDF on the redshift evolution parameter λ , mass power-law slope α , and maximum mass m max , marginalized over the local rate parameter R 0 , and assuming a flat prior on λ , α , and m 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 hyper-parameters assumed for the fixed-parameter powerlaw model. We infer a redshift evolution that is consistent with a flat in comoving volume and source-frame time merger rate ( λ = 0) with a preference for λ ≥ 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 is 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 glyph[lessorsimilar] 45 M glyph[circledot] .  \n<!-- image -->  \nin 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 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 α and m max in Section 3 are correlated by the same effect. Compared to the constraints on α and m max discussed in Section 3, which  \nassume 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 α ( m max is sufficiently well-measured that the correlation with λ is not as noticeable). Under the assumption of a constant merger rate density, Model A in Section 3 finds α = 0 . 4 +1 . 4 -1 . 9 , m max = 41 . 6 +9 . 6 -4 . 3 M glyph[circledot] , whereas allowing for redshift evolution yields α = 1 . 8 +1 . 7 -2 . 0 , m max = 41 +11 -5 M glyph[circledot] 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 α = 0 . 9 +1 . 8 -2 . 2 , m max = 38 +10 -4 M glyph[circledot] .  \nMarginalizing over the two mass distribution parameters and the redshift-evolution parameter, the merger rate density is consistent with being constant in redshift ( λ = 0), and in particular, it is consistent with the rate estimates recovered under the different mass distribution models in subsection 4.1 above. However, we find a preference for a merger rate density that increases at higher redshift ( λ ≥ 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 black holes (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 glyph[greaterorsimilar] 0 . 5, close to the O1-O2 detection horizon. Although GW170729 shifts the posterior towards larger values of λ , 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 λ = 8 . 4 +9 . 6 -9 . 5 at 90% credibility, compared to λ = 2 . 3 +9 . 9 -10 . 9 when excluding GW170729 from the analysis. With only 10 BBH detections so far, the wide range of possible values for λ 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 χ eff . For binaries  \nwhich are near equal mass, we can see from Equation 1 that only when black hole spins are high and aligned with the orbital angular momentum χ eff will be measurably greater than zero. Figure 5 in Abbott et al. (2018) illustrates the inferred χ 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 χ eff ; the rest of the posteriors cluster around χ eff = 0.  \nDespite these degeneracies, several tests have been proposed to use spins to constrain BBH formation channels (Vitale et al. 2017; Farr et al. 2017, 2018; Stevenson et al. 2017a; Talbot & Thrane 2017; Gerosa & Berti 2017; Wysocki et al. 2018; Gerosa et al. 2018). Drawing upon these methods, we now seek to estimate the black hole 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.  \nThe inferred distributions of spin magnitude are shown in Figure 8. The top panel shows the PPD as well as the median and associated uncertainties on the spin magnitude inferred from the parametric Mixture model defined in Section 2.4 and using prior distributions shown in 6. It marginalizes over all other parameters, including the mass parameters in Model C, and the spin mixture fraction. We observe that spin distributions which decline with increasing magnitude are preferred. In terms of our Beta function parameterization E [ a ] and Var[ a ], defined in Equation 5 these have mean spin E [ a ] < 1 / 2 or equivalently have β a > α a , at posterior probability 0 . 79. We find that 90% of black hole spins in BBHs are less than a ≤ 0 . 55 from the PPD, and 50% of black hole spins are less than a ≤ 0 . 27. We find similar conclusions if both black hole spins are drawn from different distributions (i.e., 90% of black hole spins on the more massive black hole are less than 0 . 6). When avoiding singular values in  \nthe 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.  \nWe also compute the posterior distribution for the magnitude of black hole spins from χ eff measurements by modeling the distribution of black hole spin magnitudes non-parametrically 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 preference for small black hole spin, inferring 90% of black holes 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 black hole spin magnitudes below 0 . 8 +0 . 15 -0 . 24 . Thus, the non-parametric 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).  \nFigure 9 shows the inferred distribution of the primary spin tilt for the more massive black hole. 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 ( ζ = 1), all black hole spin orientations are drawn from spin tilt distributions which are preferentially aligned and parameterized with σ i . In that model, the σ i distributions do not differ appreciably from the their flat priors. As such, the inferred spin tilt distribution are influenced by large σ i and the result resembles an isotropic distribution. The Mixture distribution does not return a decisive measurement of the mixture fraction, obtaining ζ = 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 ζ = 1 is ln BF = 0 . 15, indicating virtually no preference for any particular orientation distribution. While we allow both black holes to have different typical misalignment, the  \nTable 6. Summary of spin distribution models examined in Section 5.1, with prior ranges for the population parameters determining the spin models. The fixed parameters are in bold. Each of these distributions is uniform over the stated range, with boundary conditions such that the inferred parameters α a , β a must be ≥ 1. Details of the mass model listed here is described in Table 2.  \ninference 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.  \nThe mixture fraction distribution is also modelled 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.  \nIn general, we are not able to place strong constraints on the distribution of spin orientations. We elaborate in Appendix B.3 on how our black hole spin measurements are not yet informative enough to discern between isotropic and aligned orientation distribution via χ eff .", '5.2. Interpretation of Spin Distributions': "The spins of black holes are affected by a number of uncertain processes which 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 black hole should be related to the rotation of the core of its progenitor star, the amount of spin which 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; Postnov & Kuranov 2019; Arca Sedda & Benacquista 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 (Shu & Lubow 1981; Packet 1981), or 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).  \nOnce a black hole is formed, however, changing the spin magnitude is more difficult due to limitations on  \nmass accretion rates affecting how much a black hole can be spun up (Thorne 1974; Valsecchi et al. 2010; Wong et al. 2012; Qin et al. 2019). Once the binary black hole system is formed, the spin magnitudes do not change appreciably over the inspiral (Farr et al. 2014).  \nNo BBH detected to date has a component with confidently high and aligned spin magnitude. The results in the previous section imply that black holes 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 χ eff is small.  \nThe 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 which 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 binary black holes more massive than GW170729 have been observed to date, but several binaries have component masses likely between 20 -40 M glyph[circledot] . No highly asymmetric (small q ) system has been observed. Only two systems (GW151226 and GW170729) produce a χ eff distribution which is confidently different from zero; conversely, most BH binaries are consistent with χ eff near zero. These features drive our inferences about the mass and spin distribution.  \nDespite exploring a wide range of mass and spin distributions, we find the BBH merger rate density is R = 64 . 0 +73 . 5 -33 . 0 Gpc -3 yr -1 for Model A and is within R = 53 . 2 +55 . 8 -28 . 2 Gpc -3 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. 2018). We find a significant reduction in the merger rate for binary black holes with primary masses larger than ∼ 45 M glyph[circledot] . We do not have enough sensitivity to binaries with a black hole mass less than  \nFigure 8. Inferred distribution of spin magnitude for a parametric (top) and non-parametric 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 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 non-singular) results. The number of bins in the model were 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.  \n<!-- image -->  \n5 M glyph[circledot] to be able to place meaningful constraints on the minimum mass of black holes. We find mild evidence that the mass distribution of coalescing black holes may not be a pure power law, instead being slightly better fit by a model including a broad gaussian distribution  \nat high mass. We find the best-fitting models preferentially produce comparable-mass binaries (i.e., β q > 0 is preferred).  \nThe mass models in this work supersede results from an older model from O1 which inferred only the power law index (Abbott et al. 2016c, 2017e). That model found systematically larger values of α 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 glyph[circledot] , respectively. This extreme m 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 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.  \nInferring 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 which 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; Rodriguez et al. 2016b; Antonini & Rasio 2016; Mandel & de Mink 2016; Inayoshi et al. 2016; Mapelli et al. 2017; Bartos et al. 2017; Kruckow et al. 2018) produce event rates which are compatible with those from the previous observing runs (Abbott et al. 2018) 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  \nFigure 9. Inferred distribution of cosine spin tilt for the more massive black hole for two choices of prior (see Section 2.4). The dash-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.  \n<!-- image -->  \nwill allow for an exquisite probe into this dependence at large redshifts (Sathyaprakash et al. 2012; Van Den Broeck 2014; Vitale & Farr 2018).  \nWe 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 which emphasize lower spin magnitudes. Our inferences suggest 90% of coalescing black hole binaries are formed with χ eff < 0 . 3. Low spins argue against so-called second generation mergers, where at least one of the components of the binary is a black hole formed from a previous merger (Gonz'alez et al. 2007; Berti et al. 2007) and possesses spins near 0.7 (Fishbach et al. 2017).  \nGW170729 is notable in several ways: it is the most massive, largest χ 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 which 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 max . It also impacts our conclusions about redshift evolution, where its absence flattens the inferred redshift evolution.  \nRecent modelling using only the first six released events (Wysocki et al. 2018; Roulet & Zaldarriaga 2019) have 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 which 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).  \nSeveral studies have noted that population features (Mandel & O'Shaughnessy 2010; Stevenson et al. 2015; Fishbach & Holz 2017; Stevenson et al. 2017a; Zevin et al. 2017; Kovetz et al. 2017; Farr et al. 2017; Talbot & Thrane 2017; Fishbach et al. 2017; Gerosa & Berti 2017; Talbot & Thrane 2018; Farr et al. 2018; Barrett et al. 2018; Wysocki et al. 2018; Gerosa et al. 2018) and complementary physics (Abbott et al. 2016f; Zevin et al. 2017; Stevenson et al. 2017a; Chen et al. 2018) will be increasingly accessible as observations accumulate.  \nAdditional 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. 2018), 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 which were previously unobtainable.  \nThe 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'on, the Vicepresid'encia i Conselleria d'Innovaci'o, Recerca i Turisme and the Conselleria d'Educaci'o i Universitat del Govern de les Illes Balears, the Conselleria d'Educaci'o, Investigaci'o, 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 ˆ Ile-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  \nthe 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  \nCouncil 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.", 'REFERENCES': '2018, Bayesian Astrophysics, Canary Islands Winter School of Astrophysics (Cambridge University Press),  \ndoi: 10.1017/9781316182406.  \nhttps://doi.org/10.1017/9781316182406  \nAbbott, B. P., Abbott, R., Abbott, T. D., et al. 2016a, Phys. Rev. Lett., 116, 131103,  \ndoi: 10.1103/PhysRevLett.116.131103  \n-. 2016b, Phys. Rev. Lett., 116, 241103,  \ndoi: 10.1103/PhysRevLett.116.241103  \n-. 2016c, Physical Review X, 6, 041015,  \ndoi: 10.1103/PhysRevX.6.041015  \n- -. 2016d, ApJ, 818, L22,  \ndoi: 10.3847/2041-8205/818/2/L22  \n- -. 2016e, ApJ, 833, L1, doi: 10.3847/2041-8205/833/1/L1\n- -. 2016f, Phys. Rev. Lett., 116, 131102,  \ndoi: 10.1103/PhysRevLett.116.131102  \n- -. 2017a, ApJ, 851, L35, doi: 10.3847/2041-8213/aa9f0c\n- -. 2017b, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101\n- -. 2017c, ApJ, 848, L13, doi: 10.3847/2041-8213/aa920c\n- -. 2017d, Classical and Quantum Gravity, 34, 104002, doi: 10.1088/1361-6382/aa6854\n- -. 2017e, Phys. Rev. Lett., 118, 221101, doi: 10.1103/PhysRevLett.118.221101  \nAbbott, B. P., et al. 2018, To be submitted  \nAbbott, B. P., Abbott, R., Abbott, T. D., et al. 2018,  \nLiving Reviews in Relativity, 21, 3,  \ndoi: 10.1007/s41114-018-0012-9  \nAcernese, F., Agathos, M., Agatsuma, K., et al. 2015,  \nClassical and Quantum Gravity, 32, 024001,  \ndoi: 10.1088/0264-9381/32/2/024001  \nAjith, P., Hannam, M., Husa, S., et al. 2011,  \nPhys. Rev. Lett., 106, 241101,  \ndoi: 10.1103/PhysRevLett.106.241101  \nAlexander, K. D., Berger, E., Fong, W., et al. 2017, ApJ, 848, L21, doi: 10.3847/2041-8213/aa905d  \nAli-Ha¨ımoud, Y., Kovetz, E. D., & Kamionkowski, M. 2017, Phys. Rev. D, 96, 123523,  \ndoi: 10.1103/PhysRevD.96.123523  \nZevin, M., Pankow, C., Rodriguez, C. L., et al. 2017, ApJ, 846, 82, doi: 10.3847/1538-4357/aa8408  \nZiosi, B. M., Mapelli, M., Branchesi, M., & Tormen, G. 2014, MNRAS, 441, 3703, doi: 10.1093/mnras/stu824', '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 individual events, and the second is in the estimation of the detection efficiency.', 'A.1. Waveform systematics': "In Abbott et al. (2018), 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; Khan et al. 2016; Husa 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 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 glyph[circledot] 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 glyph[circledot] , 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 λ = 8 . 4 +9 . 6 -9 . 7 under the SEOBNRv3 waveform compared to λ = 8 . 4 +9 . 6 -9 . 5 under the IMRPhenomPv2 waveform. In summary, for the mass and rate part of the distributions, we find 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 which 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  \nFigure 10. Inferred distribution of cosine spin tilt for the more massive black hole for two choices of prior (see Section 2.4) with the SEOBNRv3 waveform model, with the same definitions as Figure 9.  \n<!-- image -->  \ndetector 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. (2018) to measure the total sensitive spacetime volume 〈 V T 〉 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 multi-dimensional population hyper-parameter space considered in this work. For this reason, for the parametric population studies in this work, we employ a semi-analytic method to estimate the fraction of found detections as a function of masses, spins and redshift (or equivalently, distance).  \nOur estimates of the network sensitivity are based on the semi-analytic 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 an SNR ρ ≥ ρ th in a single detector, where the threshold SNR, ρ 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 SNR, ρ opt , of the signal emitted by the BBH merger. The optimal SNR corresponds to the SNR of the signal produced by a face-on, directly overhead BBH merger with the same masses, spins and redshift. Given ρ opt , the distribution of single-detector SNRs ρ - corresponding to sources with random orientations with respect to the detector - can be calculated using the analytic distribution of angular factors Θ ≡ ρ/ρ 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 , χ 1 , χ 2 , z ), is given by the probability that ρ ≥ ρ th , or equivalently, that a randomly drawn Θ ≥ ρ th /ρ opt ( m 1 , m 2 , χ 1 , χ 2 , z ). P det referred to in this section is equivalent to the f ( z | ξ ) that appears in Equation 12 of Section 2.  \nThe semi-analytic calculation relies on two main simplifying assumptions: the detection threshold ρ 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. (2018) during O1 and the first few months of O2, and we fixed ρ th = 8. We refer to the sensitivity estimate under these assumptions as the raw semi-analytic calculation. In reality, the detector PSD fluctuates throughout the observing period. Additionally, the fixed detection threshold on SNR 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 (Nitz et al. 2017; Messick et al. 2017; Abbott et al. 2018). Consequently, the sensitive spacetime volume of a population estimated using an SNR threshold may differ from the one obtained using injections, which threshold on the pipeline-dependent detection statistic.  \nWe 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.  \nIn the first modification, which we employ throughout the mass distribution analysis (Section 3), we calibrate the raw semi-analytic method to the injection campaign in Abbott et al. (2018). The calibration takes the form of massdependent 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 〈 V T 〉 i ≡ ∫ dξp ( ξ | θ i ) V T true ( ξ ) for a set of reference hyperparameters θ i (here, mass distribution models with different exponents α and maximum masses m max ), where ξ denotes all binary parameters. To calculate 〈 V T 〉 i from injections into the PyCBC detection pipeline, we consider injections to be 'detected' if they have a detection statistic glyph[rho1] ≥ 8, where glyph[rho1] is the statistic used in the PyCBC analysis of O2 data (Nitz et al. 2017; Abbott et al. 2018). This is comparable to the detection statistic glyph[rho1] = 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. (2018), which does not fix a detection threshold, instead assigning to each trigger a probability of astrophysical origin (Farr et al. 2015). Once we have computed 〈 V T 〉 i , we correct the raw semi-analytic model V T raw described above by a factor f ( ξ ) which is a low-order polynomial in ξ : f ( ξ ) = ∑ α λ α F α ( ξ ), with F α the relevant basis polynomials. We minimize the mean-  \ncomputed by injections and ∫ dξf ( ξ ) 〈 V T 〉 raw ( ξ ) p ( ξ | θ ). If H is the precomputed matrix of weight 'moments' H k,α = ∫ dξp ( ξ | θ k ) V T ( ξ ) F α ( ξ ) then the coefficients of this least-squares expression can be computed analytically as λ = ( H T γH ) -1 H T γ 〈 V T 〉 where γ is a diagonal inverse covariance matrix characterizing the Monte Carlo integration errors of each individual 〈 V T 〉 i . This procedure yields the mass-calibrated sensitive volume 〈 V T 〉 cal .  \nThe top panel of Figure 11 shows the comparison between the raw semi-analytic 〈 V T 〉 , the calibrated 〈 V T 〉 , and the injection 〈 V T 〉 across the two-dimensional hyperparameter space of Model A for the mass distribution. We have repeated our mass distribution analysis with different choices of the 〈 V T 〉 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.  \nAs shown in Figure 11, the main effect of this calibration is to decrease 〈 V T 〉 by a factor of ∼ 1 . 6. Over the relevant part of parameter-space (i.e. the regions of the α -m 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 〈 V T 〉 calibration, although the overall rate is increased by about a factor of ∼ 1 . 6 compared to the raw semi-analytic calculation. We have verified this explicitly by repeating the analysis with and without calibrated 〈 V T 〉 .  \nFor the redshift evolution analysis (Section 4), it is not sufficient to calibrate the mass-dependence of the detection probability; we must verify that the semi-analytic calculation reproduces the proper redshift-dependence. Therefore, we pursue an alternative modification to the raw semi-analytic calculation. In this modification, we replace the single PSD of the raw semi-analytic calculation with a different PSD calculated for the Livingston detector for each fiveday 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 SNR 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.  \nWe 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 ∼ 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.  \nAnother 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 〈 V T 〉 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 det and the corresponding sensitive spacetime volume 〈 V T 〉 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 〈 V T 〉 do have a moderate influence on inferences about the spin tilt angles, presented in 5.1. When considering the effects of spin with 〈 V T 〉 , 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 does 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 (for example, choosing a different prior range of σ i as compared to 6).  \nWe 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 SNRs and detection statistics for the same sets of physical parameters. To bound the significance of this effect, we carry out the injectionbased 〈 V T 〉 estimation for both the IMRPhenomPv2 and SEOBNRv2 waveforms and find that for populations described by the two-parameter mass Model A, the waveforms produce 〈 V T 〉 estimates consistent to 10% across the relevant  \nFigure 11. Ratio between the raw semi-analytic computation of 〈 V T 〉 to the 〈 V T 〉 computed by injections into the PyCBC search pipeline (top left panel), and the same ratio for the mass-calibrated 〈 V T 〉 (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 〈 V T 〉 derived with the time-varying PSDs applied to the semi-analytic calculation V T psd to the raw and the calibrated 〈 V T 〉 . The 〈 V T 〉 for the injections is calculated for a threshold of glyph[rho1] = 8 . 0, where glyph[rho1] is the signal-noise 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 glyph[rho1] = 8 . 7 in PyCBC. We use the mass-calibrated V T cal for the parametric mass- and spin-distribution analyses in Section 3 and 5 in order to better match the injection results, while the redshift evolution analysis uses V T 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 hyper-parameters is high, indicating that systematic uncertainties in the 〈 V T 〉 estimation do not have a large impact on our results.  \n<!-- image -->  \nTable 7. Summary of intervals for each of the parameters considered in the models of Sections 3, 4 5. The reference uses the posteriors derived from the IMRPhenomPv2 waveform model and without spin effects included in 〈 V T 〉 . The second column allows for spin effects in 〈 V T 〉 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 〈 V T 〉 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 〈 V T 〉 , 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.  \nregion 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 〈 V T 〉 estimation.  \nFinally, an additional systematic uncertainty we have neglected in the 〈 V T 〉 and parametric rates calculations is the calibration uncertainty. While the event posterior samples have incorporated a marginalization over uncertainties on the calibration Farr et al. (2015) for both strain amplitude and phase, the 〈 V T 〉 estimation here does not. The amplitude calibration uncertainty results in an 18% volume uncertainty (Abbott et al. 2018), which is currently below the level of statistical uncertainty in our population-averaged merger rate estimate.", '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.  \nFigure 12. Redshift distribution of injections recovered with a false alarm rate (FAR) less than 0.1 yr -1 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 semi-analytic calculation used for the redshift evolution analysis, as described in the text. The underlying redshift distribution of the injected populations are assumed to follow a uniform in comoving volume and source-frame time distribution. The FAR threshold of 0.1 / year nearly matches the FAR of the lowest-significance GW event, GW170729, with a FAR of 0.18 / year in the GstLAL pipeline. The semi-analytic calculation closely predicts the redshift distribution of the found injections.  \n<!-- image -->', '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); Tiwari et al. (2018). Four discrete spin magnitude models are considered, the first three being special cases of Equation 4:  \n- · Low (L): p ( a ) = 2(1 -a ), i.e., α a = 1, β a = 2.\n- · Flat (F): p ( a ) = 1, i.e., α a = 1, β a = 1.\n- · High (H): p ( a ) = 2 a , i.e., α a = 2, β a = 1.\n- · Very low (V): p ( a ) ∝ e -( a/ 0 . 2)  \nSuch magnitude distributions are chosen as simple representations of low, moderate and highly spinning individual black holes. 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:  \n- · Isotropic (I): p (cos t i ) = 1 / 2; -1 < cos t i < 1, i.e., ζ = 0.\n- · Aligned (A): p (cos t i ) = δ (cos t i -1), i.e., ζ = 1, σ i = 0.\n- · 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.  \nThe isotropic distribution is motivated by dynamical or similarly disordered assembly scenarios, while the aligned one better capture 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 χ 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 χ eff distribution can be generated with nonprecessing spins.  \nHere we perform our inference entirely through χ eff , whose 12 different distributions are illustrated in Figure 13. Since we do not have conclusive results on β q from Figure 3, we cannot make a single simplifying assumption on the mass model, which the χ eff distribution depends on. We therefore consider three limiting cases: two of these  \nFigure 13. Upper row: p ( χ eff ) under various model assumptions. Labels in each subpanel legend correspond to the tilt and magnitude models defined in B.1. Isotropic models (left) provide support for both negative and positive χ eff . Aligned models (center) assume perfect alignment for each of the four magnitudes distributions. Restricted models (right) have the same shape as the Isotropic ones, with support over χ eff > 0 only. However they can be generated with nonprecessing spins. Bottom row: Posterior on the mixture fraction ζ between isotropic and aligned distributions. ζ = 0 corresponds to a completely isotropic distribution.  \n<!-- image -->  \nTable 8. Natural log Bayes factors for various spin distributions . 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.  \nfix the mass ratio to fiducial values, q = 1 and q = 0 . 5. The third corresponds to a fixed parameter model with α = 1 , m min = 5 , m max = 50. Figure 13 illustrates the χ eff distributions implied by each of these scenarios.  \nFollowing Farr et al. (2017); 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.  \nBecause of degeneracies in the GW waveform between mass ratio and χ eff , the choice of mass distribution impacts inferences about spins. This effect explains the significant difference in Bayes factors for the third row in the table. We find again our result moderately favors small black hole spins. The restricted models with χ 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 α = 2 . 3, m min = 5, m 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 black holes spins.  \nWe assume that the aligned and isotropic components follow the same spin magnitude distribution. It is possible that black holes 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.  \nWe compute the posterior on the fraction of aligned binaries ζ in the population as per Equation 6 in the limit ( σ i → 0). The models here are subsets of the Mixture distribution, with a purely isotropic being ζ = 0, and completely aligned being ζ = 1. The prior on the mixing fraction is flat.  \nAll of the models which 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 which are not distinguishable to within measurement uncertainty of χ eff . We do not include the mostfavored restricted (R) configuration, but 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) which 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 χ eff': 'We illustrate here how χ eff measurements can provide insights into discerning spin orientation distributions. Following Farr et al. (2018), we split the range of χ eff into three bins. One encompasses the fraction of uninformative binaries with χ eff consistent with zero ( | χ eff | ≤ 0 . 05); the vertical axis of Figure 14 shows the fraction of binaries lying outside of this bin. The other two capture significantly positive ( χ eff > 0 . 05), and significantly negative ( χ eff < -0 . 05) binaries. The width 0 . 05 is chosen to be of the order of the uncertainty in a typical event posterior.  \nThe aligned spin scenario is preferred in the posterior support on the right half of Figure 14: the small fraction of binaries which are informative tend to possess χ eff greater than zero. Conversely, if the spins are isotropic, there would be no preference for positive or negative χ eff , and the posterior in Figure 14 would peak towards the middle. However, of the ten observed binaries, eight are consistent with zero χ eff and only two are informative, thus demonstrating our ability to distinguish between the two scenarios is weak.', 'C. IMPORTANCE RESAMPLING THE SINGLE-EVENT LIKELIHOOD': 'Our hierarchical population analysis uses the individual-event likelihood for each event n = 1 , . . . , N , L ( d n | ξ, z ) (see Section 2, Eq. (11)). Individual-event analyses report posterior samples drawn a density that is proportional to this likelihood times a prior (Veitch et al. 2015; Abbott et al. 2018). 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  \nFigure 14. Posterior distribution for the fraction of informative binaries (i.e., | χ eff | > 0 . 05), and the fraction of those informative binaries with positive χ eff (i.e., χ eff > 0 . 05).  \n<!-- image -->  \nframe masses and redshift, the prior is  \np ( m source 1 , m source 2 , z ) ∝ d 2 L ( z ) ∂m det 1 ∂m source 1 ∂m det 2 ∂m source 2 ∂d L ∂z = d 2 L ( z ) (1 + z ) 2 ∂d L ∂z . (C1)  \nThe derivative of the luminosity distance in a spatially flat universe (Hogg 1999) is  \n∂d L ∂z = d L 1 + z +(1 + z ) d H E ( z ) , (C2)  \nwhere d H = c/H 0 is the Hubble distance and  \nE ( z ) ≡ √ Ω M (1 + z ) 2 +Ω Λ (C3)  \nin a ΛCDM universe.  \nGiven 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  \nw ( m source 1 , m source 2 , z ) = 1 p ( m source 1 , m source 2 , z ) . (C4)  \nThe integral in Eq. (11) may then be approximated as  \n∫ d ξ d z L ( d n | ξ, z ) d N d ξ d z ∝ N samp ∑ i =1 w ( m source 1 ,i , m source 2 ,i , z i ) d N d ξ i d z i , (C5)  \nwhere 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. ξ = ( m source 1 , m source 2 , . . . )). 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 source 1 ,i , m source 2 ,i , z i ) ; the integral is then proportional to the average value of d N/ d ξ d z 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.'}