Formation and coalescence sites of the first GW events

The formation and coalescence sites of the first gravitational wave events


We present a novel theoretical model to characterize the formation and coalescence sites of compact binaries in a cosmological context. This is based on the coupling between the binary population synthesis code SeBa with a simulation following the formation of a Milky Way-like halo in a well resolved cosmic volume of 4 cMpc, performed with the GAMESH pipeline. We have applied this technique to investigate when and where systems with properties similar to the recently observed LIGO/VIRGO events are more likely to form and where they are more likely to reside when they coalesce. We find that more than 70% of GW151226 and LVT151012-like systems form in galaxies with stellar mass in the redshift range [0.06 - 3] and [0.14 - 11.3], respectively. All GW150914-like systems form in low-metallicity dwarfs with at . Despite these initial differences, by the time they reach coalescence the observed events are most likely hosted by star forming galaxies with . Due to tidal stripping and radiative feedback, a non negligible fraction of GW150914-like candidates end-up in galaxies with properties similar to dwarf spheroidals and ultra-faint satellites.

galaxies: evolution, high-redshift, black hole physics – gravitational waves – stars:black holes

1 Introduction

The recent detection of gravitational waves from the LIGO/VIRGO collaboration has opened the new era of gravitational wave astronomy (Abbott et al., 2016a, b). During the first Advanced LIGO observing run (hereafter O1, from September 12 2015, to January 19 2016), two sources have been unambiguously detected (GW150914 and GW151226), while a third one (LVT151012) was below detection threshold, but with 87% probability of being of astrophysical origin (Abbott et al., 2016c). At the time of submission of this letter, the detection of a fourth event has been reported, GW1701043(Abbott et al., 2017).

All the four sources are powered by the inspiral and coalescence of two black holes (BHs) and the binary properties, estimated from the gravitational wave signals, have profound astrophysical implications (Abbott et al., 2016d). In particular, GW150914 originated from the merger of two BHs at , suggesting that stellar BHs with such large masses can form in nature, evolve in a binary system, and coalesce within a Hubble time. Comparison with stellar and binary evolution models suggests progenitor stars with low metallicity, (Abbott et al., 2016a; Belczynski et al., 2016; Mapelli, 2016). In fact, reduced mass loss at low metallicity favors the formation of more massive stellar remnants (Mapelli, Colpi & Zampieri, 2009; Mapelli et al., 2010; Belczynski et al., 2010; Spera, Mapelli & Bressan, 2015).

The limit on the metallicity potentially provides very interesting constraints on the birth environment of massive BH binaries. If the binary BH system has merged in a short time, it must have formed in a rare, low-metallicity dwarf galaxy. Alternatively, in the long merger time scenario ( Gyr), it could have formed at high redshift (), where low-metallicity star formation is expected to be more common (Abbott et al., 2016d).

The main goal of this paper is to investigate when and where massive BH binary systems with properties similar to GW150914, GW151226, and LVT151012 are more likely to form, and where they are more likely to reside at the time of their coalescence.

Following its discovery, several attempts have been made to determine the environment in which GW150914 formed. Belczynski et al. (2016) and Dvorkin et al. (2016) use different massive BH formation scenarios and the observationally inferred cosmic star formation rate density evolution with a metallicity-dependent correction to estimate the redshift evolution of birth and merger rate of massive BH binaries. Although these studies may provide useful constraints on the redshift dependent birth and merger histories of compact binary systems, (Schneider et al., 2001; Regimbau, 2011; Marassi et al., 2011; Dominik et al., 2013), they are not capable of discriminating individual formation or coalescence sites. Using observationally inferred galaxy scaling relations Lamberts et al. (2016) and Elbert, Bullock & Kaplinghat (2017) predict the low redshift BH binary merger rate as a function of the present-day host galaxy mass. They find that GW150914-like events with short merger timescales would be primarily localized in dwarf galaxies, and in massive galaxies otherwise. Conversely, using cosmological simulations of galaxies with different masses at , O’Shaughnessy et al. (2016) find that a present-day dwarf galaxy can host a considerably larger BH binary merger rate compared to a massive galaxy. However, some of the galaxy scaling relations adopted by Lamberts et al. (2016) are not directly observed for low-mass faint galaxies, and are therefore based on extrapolations. This may introduce considerable uncertainties in the estimated properties of the most likely sites of low-metallicity star formation, hence of binary BH formation. The analyses of O’Shaughnessy et al. (2016) and Elbert, Bullock & Kaplinghat (2017) adopt a delay time distribution function to characterize the merger time of binary BHs and do not determine the properties of individual binary systems formed in each galaxy.

In this work, we attempt to overcome the above limitations, and to use individual galaxy properties to predict the number and properties of binary systems that each galaxy hosts.

Figure 1: The number of events like GW150914, GW151226, and LVT151012 (see the legenda) predicted by SeBa as a function of metallicity (left panel) and merger time (right panel).

2 Forming BH-BH systems

The binary population synthesis model:
We adopt the binary population synthesis code SeBa4 originally developed by Portegies Zwart & Verbunt (1996) and Nelemans, Yungelson & Portegies Zwart (2001). SeBa follows the evolution of binary systems by taking into account the relevant physics involved in their evolution: stellar composition, stellar winds, mass transfer and accretion, magnetic braking, common envelope phase, supernova kicks, and gravitational radiation. The original version has been recently modified by Mapelli et al. (2013) to include metallicity-dependent prescriptions for stellar evolution, stellar winds, and remnant formation. In particular, metallicity-dependent stellar evolution is based on the polynomial fitting formulas described in Hurley, Pols & Tout (2000). Stellar winds were substantially upgraded including Vink, de Koter & Lamers (2001) prescriptions for O-type massive stars, Vink & de Koter (2005) formalism for Wolf-Rayet stars, and Belczynski et al. (2010) recipes for luminous blue variable stars. The mass of the compact remnant was derived as described in Mapelli, Colpi & Zampieri (2009): stars with a pre-supernova mass M are assumed to collapse to BH directly, retaining most of their final mass. Binary evolution was modelled as in the default SeBa code. For more details on single and binary evolution, see Mapelli et al. (2013) and Fig. 1 of Mapelli et al. (2013) for the derived BH mass spectrum. We set the binary fraction so that all the stars are assumed to be in binary systems. Observations show that this assumption is approximately correct for stars with masses (Sana et al., 2012). For each simulation, we initialize binaries and follow their evolution with initial properties randomly selected from independent distribution functions. In particular, the primary stellar mass, , is distributed according to a Kroupa Initial Mass Function (IMF, Kroupa 2001) between . The mass of the secondary star, , is generated according to a flat distribution for the mass ratio with . The initial semi-major axis (sma) distribution is flat in log(sma) (see Portegies Zwart & Verbunt 1996) ranging from (Roche lobe contact) up to . The eccentricity of the binary is selected from a thermal distribution in the range (Heggie, 1975). With this choice of initial conditions, we run simulations for 11 different values of the initial metallicity in the range .

An extended statistical analysis of the simulations will be presented in a forthcoming paper. Here we illustrate the properties of simulated BH-BH systems with primary and secondary masses in the range estimated for GW150914 ( and ), GW151226 ( and ), and LVT151012 ( and ). In the left panel of Fig. 1 we show the number of candidate systems for the three events found in simulated samples with different metallicities. We find only two massive BH binaries with masses compatible with GW150914 and both require very low initial stellar metallicities, 5. GW151226-like systems are more common and their stellar progenitors have metallicities in the range . Due to the higher mass of the primary BH, we find that LVT151012-like systems are less frequent and that most of their stellar progenitors have metallicities . In the right panel of Fig. 1 we show the corresponding distribution of merger times. The two GW150914-like systems are both characterized by long coalescence times . This is expected as massive BH binaries originate from massive stellar binaries that, given their larger radii, must have large semi-major axis at birth to avoid merging before the compact binary system forms (Linden et al., 2010; Mapelli & Zampieri, 2014; Ziosi et al., 2014). Instead, GW151226 and LVT151012-like systems follow a relatively flat distribution in merger times in the range . Due to the initial mass ratio in the zero age main sequence (), all the GW151226 and LVT151012-like systems during their evolution go through a phase of common envelope, which drastically reduces the initial orbital separation. For these systems, the flat merger times distribution is a consequence of the flat initial distribution of orbital separations.

The cosmological simulation:
We couple the binary information provided by SeBa with the star formation and chemical enrichment histories of all the galaxies contained in a comoving volume , experiencing a forming Milky Way (MW)-like galaxy at its center. These are obtained processing a N-Body cosmological simulation with the GAMESH pipeline (Graziani et al., 2015). The simulation has been presented in Graziani et al. (2017) and here we briefly summarize the properties relevant to the present analysis.

The N-Body simulation is performed by the code GCD+ (Kawata et al., 2013) and adopts a Planck 2013 cosmology (Planck Collaboration et al. 2014, ,, and ) and initial conditions (ICs) suitable to reproduce a MW-like halo in the present Universe, starting from a cosmological volume of about (83.53 cMpc). Once a MW-sized halo is identified at the largest scale, we re-create ICs for a multi-resolution, zoom-in simulation resolving the same halo with mass at the center of a volume of (4 cMpc) with a dark matter particle resolution mass of . Hereafter we will refer to this volume as the Local Group (LG) of the MW-like halo. The simulation runs from to and we store the outputs with a time resolution of  Myr at , and of  Myr at . In each snapshot, halos are identified by a standard FoF algorithm adopting linking parameter of and a minimum number of particles of 100. A particle-based merger tree (MT) has been computed to exactly establish each halo ancestor/descendant relationships, also accounting for all the dynamical processes regulating halo assembly: accretion, mergers, tidal stripping and halo disruption. The simulation is in good agreement with observations and with independent theoretical studies (Graziani et al., 2017). At all redshifts, the resulting halo mass function in the LG matches the prediction of the analytic Press-Schecter distribution for all halos with . Below , the central halo dynamically dominates the LG region (half of the MW halo mass is already assembled by ) and halos with mass are found to evolve faster than expected for an average cosmic region.

The baryonic evolution of galaxies is performed by GAMESH accounting for star formation, metal enrichment, Pop III/Pop II transition and supernova (SN)-driven feedback. Radiative feedback is modeled by adopting an instant reionization prescription which suppresses star formation in all galaxies hosted by mini-halos (with virial temperatures  K) found at . By calibrating the model free parameters to reproduce the stellar, gas and metal mass observed in the MW, we find the simulated galaxies to be in good agreement with recent observations of candidate MW progenitors at , with the galaxy main sequence, mass-metallicity relation and fundamental plane of metallicity relations in the redshift range .

It is worth pointing out that although the simulated cosmic volume is very small compared to the current instrumental range of the advanced LIGO detector ( Gpc), the adopted resolution allows us to simulate DM structures down to a minimum mass of , where star formation in low-metallicity environments (hence massive BH formation) is more likely to occur.

Figure 2: Properties of the formation sites of GW150914, GW151226 and LVT151012-like systems in the simulation. Left panel: percentage of systems as a function of formation redshift. Right panel: gas metallicity as a function of the stellar mass of the galaxies where the simulated LIGO-like events are formed.
Figure 3: Properties of the coalescence sites of GW150914, GW151226 and LVT151012-like systems in the simulation. Histograms show the percentage of systems hosted in galaxies binned by their star formation rate at the time of coalescence. The percentage of BH binaries hosted in galaxies with SFRs within the broader bins identified by the vertical lines are also shown, together with the average host stellar mass and metallicity. The numbers in the leftmost bin refer to BH binaries hosted in galaxies with negligible or null star formation activity, as indicated by the arrow (see text).

BH binary formation in the simulated galaxies:
In order to determine the properties of compact binaries formed in each galaxy, we assume that the stellar progenitors have the metallicity of the gas in which they form and we randomly extract from the SeBa output with the closest metallicity a number of binary systems until we reach the total mass of newly formed stars predicted by the cosmological simulation. As a result of this procedure, we can predict the metallicity-dependent formation and merger rates of compact binaries and the properties of compact binary systems hosted in each galaxy of the LG at any given time from .

3 Results

We select, among all the simulated binaries, systems with properties consistent with GW150914, GW151226 and LVT151012. This is done by requiring that the system has primary and secondary BH masses in the observed range and that it coalesces at for GW150914 and GW151226, and at for LVT151012. We find , and for GW150914, GW151226, and LVT151012, corresponding to merger rates, in the LG, of , and Gpcyr, respectively. These differ from the rates obtained by Abbott et al. (2016c) based on populations with masses matching the observed events and assuming a uniform distribution in comoving volume6. Indeed, for this comparison to be meaningful a larger simulation box probing a more representative cosmological volume should be considered (see Mapelli et al. 2017 for a similar approach using the Illustris simulation).

For each class of systems, we reconstruct from the simulation output the properties of their formation and coalescence sites. Fig. 2 shows the redshift distribution (left panel) and the mass-metallicity relation (right panel) of the host galaxies at BH binary formation. The redshift range over which the systems formed reflects the corresponding merger time distribution. GW150914-like systems are characterized by long merger times and form in the range . GW151226 and LVT151012-like systems, which have a broad range of possible merger times, can form in relatively wide redshift intervals and , respectively. However, we find that while only 6% of GW151226-like systems have , among LVT151012 systems, 67% have , 48% have , and 6% formed in the pre-reionization epoch, with . In addition, all the GW150914-like systems are born in low-metallicity dwarfs with stellar mass , while 88% (70%) of GW151226 (LVT151012)-like systems form in galaxies with .

We track each GW150914, GW151226 and LVT151012-like BH system from formation to coalescence using the particle-based MT and we identify the galaxy where it resides when the merger event occurs. The results of this procedure are illustrated in Fig. 3. The histograms represent the percentage of systems hosted in galaxies characterized by their star formation rate (SFR) at the time of coalescence. The numbers shown in the plot refer to the percentage of systems in larger SFR intervals as identified by the vertical lines. For each of these intervals, we also report the average stellar mass and metallicity of the hosts. The numbers in the leftmost bin refer to systems which are not forming stars at the time of coalescence. Our results suggest that GW150914, GW151226 and LVT151012-like systems have the largest probability to be hosted in galaxies with masses , metallicities and which form stars at moderate rates /yr. These are the progenitors of the MW-like galaxy at the redshift of coalescence (Graziani et al., 2017). Smaller galaxies, with stellar masses encompassing the range of LMC, SMC and dwarfs, have a smaller probability to have hosted the first three GW events.

Finally, we find a small but not negligible percentage of GW150914-like systems ( 10%) whose coalescence could be hosted in very small galaxies, with , where star formation has been suppressed by radiative feedback. These BH binaries have formed in dwarf galaxies at which, as a consequence of stripping events, lower their (dark and baryonic) mass below the limit that renders gas infall unfavorable (virial temperatures smaller than  K) in the rising UV background accompanying reionization. Although these findings are based on a very small number statistics (1 out of 11 GW150914-like systems), they suggest that the hosts of massive BH binary mergers might be similar to the least massive among dwarf spheroidals or to the more massive among ultra-faint galaxies.

4 Conclusions

The work that we have presented is only a first step towards a more extended analysis of different classes of compact binary systems and their detectability as gravitational wave sources. Ultimately, our goal is to predict the most probable sites of formation and coalescence of different classes of compact binary systems. This information will help the search for the host galaxy properties and to identify the electromagnetic counterparts of some gravitational wave sources.

While we plan to present an extended analysis in future work, here we select, among all the simulated binaries, GW150914, GW151226 and LVT151012-like systems. We find that, despite the fact that these systems require different galaxy properties at the time of their formation, they all have a larger probability to be hosted in massive, star forming galaxies at the time of their coalescence. This is in agreement with what was recently suggested by independent studies (Lamberts et al., 2016; Elbert, Bullock & Kaplinghat, 2017). However, massive BH binaries which preferentially form in metal-poor dwarfs at , can also fail to be incorporated in more massive galaxies during hierarchical evolution. As a result, they might have evolved within galaxies which are fragile to radiative feedback effects, reaching coalescence in small non star forming dwarf satellites.


The research leading to these results has received funding from the European Research Council under the Grant Agreement n. 306476. MM and MS acknowledge financial support from the Italian Ministry of Education, University and Research (MIUR) through grant FIRB 2012 RBFR12PM1F. MM acknowledges financial support from the MERAC Foundation.


  1. pagerange: The formation and coalescence sites of the first gravitational wave eventsReferences
  2. pubyear: 2017
  3. The analysis presented is limited to the first three events.
  5. A similar condition () applies to GW170104-like systems, given the masses of the primary () and secondary () black holes.
  6. Abbott et al. (2016c) find , and Gpcyr for GW150914, GW151226, and LVT151012, where we quote the median values with 90% confidence from their Table II.


  1. Abbott B. P. et al., 2016c, Phys. Rev. X, 6, 041015
  2. Abbott B. P. et al., 2016b, Physical Review Letters, 116, 241103
  3. Abbott B. P. et al., 2016a, Physical Review Letters, 116, 061102
  4. Abbott B. P. et al., 2016d, The Astrophysical Journal Letters, 818, L22
  5. Abbott B. P. et al., 2017, Physical Review Letters, 118, 22
  6. Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, The Astrophysical Journal, 714, 1217
  7. Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016, NATURE, 534, 512
  8. Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2013, Astrophysical Journal, 779, 72
  9. Dvorkin I., Vangioni E., Silk J., Uzan J.-P., Olive K. A., 2016, mnras, 461
  10. Elbert O. D., Bullock J. S., Kaplinghat M., 2017, ArXiv e-prints 1703.02551
  11. Graziani L., de Bennassuti, Schneider R., Kawata D., Salvadori S., 2017, MNRAS, in press
  12. Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, Monthly Notices of the Royal Astronomical Society, 449, 3137
  13. Heggie D. C., 1975, Monthly Notices of the Royal Astronomical Society, 173, 729
  14. Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  15. Kawata D., Okamoto T., Gibson B. K., Barnes D. J., Cen R., 2013, MNRAS, 428, 1968
  16. Kroupa P., 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
  17. Lamberts A., Garrison-Kimmel S., Clausen D. R., Hopkins P. F., 2016, Monthly Notices of the Royal Astronomical Society, 463, L31
  18. Linden T., Kalogera V., Sepinsky J. F., Prestwich A., Zezas A., Gallagher J. S., 2010, ApJ, 725, 1985
  19. Mapelli M., 2016, MNRAS, 459, 3432
  20. Mapelli M., Colpi M., Zampieri L., 2009, MNRAS, 395, L71
  21. Mapelli M., Giacobbo N., Ripamonti E., Spera M., 2017, MNRAS
  22. Mapelli M., Ripamonti E., Zampieri L., Colpi M., Bressan A., 2010, MNRAS, 408, 243
  23. Mapelli M., Zampieri L., 2014, ApJ, 794, 7
  24. Mapelli M., Zampieri L., Ripamonti E., Bressan A., 2013, Monthly Notices of the Royal Astronomical Society, 429, 2298
  25. Marassi S., Schneider R., Corvino G., Ferrari V., Portegies Zwart S., 2011, Physical Review D, 84, 124037
  26. Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001, Astronomy and Astrophysics, 375, 890
  27. O’Shaughnessy R., Bellovary J. M., Brooks A., Shen S., Governato F., Christensen C. R., 2016, Monthly Notices of the Royal Astronomical Society
  28. Planck Collaboration et al., 2014, Astronomy & Astrophysics, 571, A16
  29. Portegies Zwart S. F., Verbunt F., 1996, Astronomy and Astrophysics, 309, 179
  30. Regimbau T., 2011, Research in Astronomy and Astrophysics, 11, 369
  31. Sana H. et al., 2012, Science, 337, 444
  32. Schneider R., Ferrari V., Matarrese S., Portegies Zwart S. F., 2001, Monthly Notices of the Royal Astronomical Society, 324, 797
  33. Spera M., Mapelli M., Bressan A., 2015, Monthly Notices of the Royal Astronomical Society, 451, 4086
  34. Vink J. S., de Koter A., 2005, A&A, 442, 587
  35. Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, Astronomy and Astrophysics, 369, 574
  36. Ziosi B. M., Mapelli M., Branchesi M., Tormen G., 2014, MNRAS, 441, 3703
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description