BBH in the Milky Way

Predicting the binary black hole population of the Milky Way with cosmological simulations


Binary black holes are the primary endpoint of massive stellar evolution. Their properties provide a unique opportunity to constrain binary evolution, which is still poorly understood. In this paper, we predict the inventory of binary black holes and their merger products in/around the Milky Way, and detail their main properties. We present the first combination of a high-resolution cosmological simulation of a Milky Way-mass galaxy with a binary population synthesis model. The hydrodynamic simulation, taken from the FIRE project, provides a cosmologically realistic star formation history for the galaxy and its stellar halo and satellites. We apply a metallicity-dependent evolutionary model to the star particles to produce individual binary black holes. We find that a million binary black holes have merged in the model Milky Way, and 3 million binaries are still present, with an average mass of 28 per binary. Because the black hole progenitors are biased towards low metallicity stars, half reside in the stellar halo and satellites and 40 per cent of the binaries were formed outside the main galaxy. This trend increases with the masses of the black holes. The numbers and mass distribution of the merged systems is compatible with the LIGO/Virgo detections. Observations of these black holes will be challenging, both with electromagnetic methods and LISA. We find that a cosmologically realistic star formation history, with self-consistent metal enrichment and Galactic accretion history, are key ingredients for determining binary black hole rates that can be compared with observations to constrain massive binary evolution.

stars:black holes; binaries:close, Galaxy: stellar content, abundances, gravitational waves

1 Introduction

The global properties of compact objects (CO) in the Milky Way (MW) provide crucial information on the star formation history of the Galaxy as well as on stellar evolution. Broadly speaking, stars born with mass evolve into white dwarfs (WD), those of evolve into neutron stars (NS), and those above turn into black holes (BH; e.g. Fryer 1999, but also see Sukhbold et al., 2016), though stars between may undergo a pair instability supernova that leaves no remnant (Fryer et al., 2012). Aside from systems that have undergone mergers or left the Galaxy due to BH kicks (Janka, 2013), the number of compact remnants quantifies past star formation. The localization of COs within a galaxy may also be indicative of the progenitor’s formation conditions (lookback time, local environment, and metallicity). The mass distributions, orbital properties and/or proper motion of the COs can inform us on stellar evolution and explosion mechanisms. In this paper, we provide detailed predictions for the expected binary black hole (BBH) population and their merger products in the MW, as well as some observational properties.

The current () population of BHs is particularly important as BHs evolve from the most massive stars, whose short lives mean that we can only directly observe the population that formed within the last  Myr. In particular, BHs provide unique information on the initial mass function of massive stars, which are key drivers of galactic evolution through chemical enrichment, stellar winds, ionizing radiation and their final explosions (Muratov et al., 2015; Geen et al., 2015). Observations indicate that most, if not all massive stars form in binary systems (Sana et al., 2012), with recent work suggesting that binaries are even more ubiquitous for lower metallicity stars (Badenes et al., 2017). As such, the properties of stellar black holes can also inform us about crucial phases of binary evolution such as supernova kicks and mass transfer (see Postnov & Yungelson 2014 for a recent review on compact binary formation).

Unfortunately, the inventory of stellar mass black holes in the MW is far from complete. So far, the only confirmed systems are found in X-ray binaries, where the BHs manifest themselves through accretion of material from a companion star (see Casares et al. 2017 for a recent review). About 60 of these systems have been detected around low-mass companion stars, and a handful around massive stars (Corral-Santana et al., 2016). The latter are systems formed within the last 20 million years and possible progenitors to BBHs. The astrometric mission GAIA could detect more than ten thousand BHs around stellar companions (Breivik et al., 2017; Mashian & Loeb, 2017). Based on the observed BBH merger rate from Abbott et al. (2016a), Elbert et al. (2017) estimate that there could be up to 100 million BHs in the MW. So far, however, no BH has been observed in a binary with another compact object in the MW.

Similarly, there have been no firm detections of stellar mass BHs without a stellar companion in the MW. Such black holes are not expected to emit electromagnetic radiation unless they are accreting from a dense environment (Agol & Kamionkowski, 2002; Maccarone, 2005). Future hard X-ray surveys or radio observations with the Square Kilometer Array (SKA) could lead to the first such detections, if the accretion rates and radiative efficiencies are high enough (Fender et al., 2013; Corbel et al., 2015). Ioka et al. (2017) further suggest that BHs formed out of merged BBHs should have a high spin and may thus produce gamma-ray emission in a jet. Year-long microlensing events with no visible lens have been tentatively attributed to BHs in the galactic bulge (Wyrzykowski et al., 2011, 2016).

Gravitational waves may be the most promising means of detecting BBHs in the MW. The first direct detection of GWs came from the merger of two stellar mass BHs (GW150914) via the Laser Interferometry Gravitational Wave Observatory (LIGO; Abbott et al., 2016c), with a handful of similar detections following. Several studies propose to distinguish binary evolution channels using observational properties of compact object mergers from gravitational waves (Belczynski et al., 2002; Voss & Tauris, 2003; Belczynski et al., 2008a; Dominik et al., 2012; Stevenson et al., 2015; Mapelli et al., 2017). With LIGO/Virgo, we observe the very last moments (less than a second) before the BBH merges, the merger itself, and the ringdown of the newly formed (and more massive) BH. Given the very short signal in comparison to the very long inspiral time of a BBH, the likelihood of detecting one of these events in the MW is close to zero.

The high masses of GW150914, strongly point toward a low-metallicity progenitor binary (Dominik et al., 2013; Belczynski et al., 2008b), with (Belczynski et al., 2016). Based on an analytic model, we showed in Lamberts et al. (2016) that the progenitors of GW150915 most likely formed either relatively recently in a dwarf galaxy (stellar mass ), or around the peak of cosmic star formation ( 10 Gyrs ago) in a galaxy that would now resemble the MW (stellar mass ). Mapelli et al. (2017) and Schneider et al. (2017) found similar results when combining a binary evolution model with the Illustris simulation and with a high-resolution dark matter-only simulation, respectively. While lower mass BBHs can be formed out of higher metallicity progenitors, BBHs are strongly biased towards sub-solar metallicity environments, as the amount of mass loss from stellar winds scales with metallicity (Belczynski et al., 2010a; Dominik et al., 2013). As such, we expect the BBH population of the MW to have a different spatial distribution than the overall stellar mass.

Current predictions for the binary compact object population of the MW are based on simplified models for its star formation history, metallicity and morphology. The MW is often approximated by a spherically symmetric bulge and a disk with a characteristic scale height (e.g. Nelemans et al. (2001); Ruiter et al. (2010); Liu & Zhang (2014)). The star formation rate in the disk is typically assumed to be constant over time and to occur at fixed metallicity (solar). The stellar halo is rarely included (except in Belczynski et al., 2010c, where the halo is assumed to form in a single burst 13 Gyr ago with a single metallicity). The inaccuracies resulting from these approximations are compounded by the uncertainties in the binary evolution models, such that it is extremely difficult to constrain binary evolution through comparisons between predictions from these simplified Galaxy models and observations.

These simplifications motivate the present analysis, where we instead apply a population synthesis model to the star formation history of a cosmological, hydrodynamic simulation of a MW-mass galaxy, first presented in Wetzel et al. (2016). Although the simulation does not specifically aim to reproduce the exact morphology of the MW, it includes a cosmologically realistic star formation and merger histories, satellite population, and stellar halo, and a self-consistent metal enrichment and gas exchange with the circumgalactic and extragalactic media. An exact prediction for the BBH distribution in the MW requires observational constraints on the ages and metallicities of the entire stellar population of the MW, and would be particularly sensitive to the oldest (and therefore faintest) stars. Unfortunately, this detailed information will remain out of reach for some time, especially outside the disk of the Galaxy. In the meantime, our technique should yield a binary black hole population thatp is statistically consistent with that of the MW.

The Laser Interferometer Space Antenna (LISA), a space-based mission led by ESA and scheduled to launch in the mid 2030’s, will provide the first view of the double compact object population in the MW. With its 2.5 million km arms, the interferometer will be mostly sensitive to gravitational wave frequencies Hz. One of LISA’s main science goals is a first inventory of very short orbit double compact objects (DCO) in the MW. WD binaries with orbits shorter than an hour are expected to vastly dominate the signal (Nelemans et al., 2001; Ruiter et al., 2010), as they stem from more common low mass stars and their formation likely has limited metallicity dependence. However, analytic estimates (Seto, 2016; Christian & Loeb, 2017) and binary population synthesis models combined with a multi-component model for the Galaxy (Belczynski et al., 2010c; Liu & Zhang, 2014) predict that at most a few BBH may be detected. However, Belczynski et al. (2010c) predict roughly BBH in the Galaxy in their model A, where binaries survive the common envelope occurring during the Hertzprung gap. They predict roughly binaries in model B, where common envelope mass transfer during the Hertzprung gap results in a stellar merger. In both cases, they predict most of the BBHs will be in the disk, a quarter of the systems in the bulge, and a negligible contribution from the halo. We will show that this arises directly from their neglect of more complex chemical and star formation structure in the Galaxy. By applying a binary population synthesis model (similar to model B of Belczynski et al. 2010c) to a more realistic stellar population (both in terms of spatial distribution and age-metallicity space), we predict at most one BBH detection in the MW with LISA.

In this paper, we provide a detailed view of the black hole population resulting from binary black holes systems in a MW-mass galaxy. We only focus on systems that survived BBH formation and are either unmerged binary black holes (UBBH) or merged binary black holes (MBBH). We do not consider black holes from single star formation, black holes resulting from binary systems disrupted by supernova kicks, black holes with a stellar companion, or single black holes formed after a stellar merger. We specifically consider UBBH+MBBH black holes and do not present earlier phases of binary evolution where BHs may have stellar companions emitting electromagnetic radiation. Fig. 1 summarizes the systems considered in this work.

Figure 1: Possible endpoints of binary massive star formation. In this paper we only focus on unmerged binary black holes (UBBH) and merged binary black holes (MBBH), in red in the bottom left. We discard black holes from single stellar evolution, stellar mergers or binaries unbound by supernova kicks. We also do not consider black holes with white dwarf, neutron star or stellar companions. The above picture is a cartoon view to clarify the systems considered here and is not intended to be a precise representation of binary evolution.

We combine a high resolution cosmological simulation of a MW-mass halo (Wetzel et al., 2016) with a binary population synthesis model (§2). We show how the progenitors of the BHs compare with the global population of stars and determine the properties of the BH population (§3). We present observational properties of both merged (i.e. currently single BHs) and unmerged systems (binary BHs), both with gravitational waves and electromagnetic methods (§4). We discuss the importance of a detailed model of the MW (§5) and conclude (§6).

2 Methods

This paper emphasizes the importance of a realistic model for the star formation history for MW-mass galaxies with respect to previous studies of binary compact objects in the MW. We first present the key numerical and physical parameters of the simulation and show the physical quantities that are the most relevant to this study (§2.1). We then present the binary evolution model we use (§2.2) as well as our computation of GW emission (§2.3). We then explain how all these aspects are combined together (§2.4).

2.1 A realistic model of a Milky Way-mass galaxy

The inputs to our binary evolution model (the ages, metallicities, and positions of star particles) are primarily drawn from the m12i FIRE-2 simulation, also known as “Latte” (Wetzel et al., 2016). The simulation has an initial gas particle mass of 7070 M. The Latte simulation is part of the Feedback in Realistic Environment (FIRE; Hopkins et al., 2014) project3, specifically run using the improved “FIRE-2” version of the code from Hopkins et al. (2017, for details, see Section 2 therein). The simulations use the code GIZMO (Hopkins, 2015),4 with hydrodynamics solved using the mesh-free Lagrangian Godunov “MFM” method. For the gas, both the hydrodynamic and gravitational (force softening) resolutions are fully adaptive down to 1 pc. The simulations include cooling and heating from a meta-galactic background and local stellar sources from K. Star formation occurs in locally self-gravitating, dense, self-shielding molecular, Jeans-unstable gas. Stellar feedback from OB & AGB star mass-loss, type Ia & II supernovae, and multi-wavelength photo-heating and radiation pressure is directly based on stellar evolution models. Chemical enrichment stems from type Ia supernova (Iwamoto et al., 1999), core-collapse supernova (Nomoto et al., 2006), and O and AGB star winds (van den Hoek & Groenewegen, 1997; Marigo, 2001; Izzard et al., 2004).

The FIRE simulations reproduce the observed mean mass-metallicity relation both for stars and star forming gas, between and (Ma et al., 2016) down to a stellar mass of . The simulations here include the subgrid-scale numerical turbulent metal diffusion terms described in Hopkins et al. (2017), which have almost no dynamical effect at the galaxy mass scales considered here (Su et al., 2017), but produce better agreement with the internal metallicity distribution functions observed in MW satellite galaxies (Escala et al., 2017).

Our main analysis is based on galaxy m12i (from Wetzel et al., 2016, though we analyze a re-simulation with turbulent metal diffusion first presented in Bonaca et al., 2017), chosen to have a relatively “normal” merger history, but we also consider a lower-resolution version of m12i as well as two different galaxies m12b and m12c (Hopkins et al., 2017) at the same mass scale. m12i shows metallicity gradients (Ma et al., 2017) and abundances of -elements (Wetzel et al, in prep.) in the disk that are broadly consistent with observations of the MW. Its global star formation history is consistent with the MW (see Ma et al. (2017) for illustrations) although its present day star formation rate of 6 yr is somewhat higher than observed in the Milky Way. The satellite distribution around the main galaxy in m12i presents a similar mass and velocity distribution as observed around the Milky Way and M31, down to a stellar mass of , though the simulation does not contain an equivalent of the Large Magellanic Cloud; the most massive satellite is comparable to the Small Magellanic Cloud.

The simulation produces a catalog of star particles.5 For each particle, the quantities of interest here are its formation time , metallicity , and position at and at . To determine the accretion rates of each BH associated with a star particle, we also recover the properties of the closest surrounding gas particle and assign it to the star particle. We consider only particles within 300 kpc of the center of the main galaxy. This is slightly larger than the virial radius of the galaxy and allows us to largely sample the halo, satellites and streams while remaining unaffected by the boundaries of the high resolution region.

Figure 2: Properties of BBHs at their formation times for different stellar metallicities. We show the total binary masses () versus initial orbital period () or equivalently the gravitational wave frequency of an equivalent circular orbit (the orbits are not necessarily circular, but this is intended only to guide the reader). The colors are normalized to the total number of BBH formed in each metallicity bin while the percentage value shows the percentage of massive binaries that evolve into BBHs. The black/grey lines across the plot shows the maximal initial period for circular systems to be able to merge, assuming they formed at the beginning of the universe or at , respectively. The green line shows the minimal period for a merger within a Hubble time for BBHs with initial eccentricity , which is the case for about 1 per cent of the systems at all metallicities. Most systems left of the black line cannot have merged by the present day.

The simulations assume a CDM cosmology with = 0.728, = 0.272, = 0.0455, h = 0.702, = 0.807, and n = 0.961 (Planck Collaboration et al., 2016). All metallicities are defined with respect to the solar metallicity, set to .

2.2 Binary evolution model

We consider binary evolution in the field, neglecting possible GW sources from N-body stellar dynamics in globular clusters or other formation channels (Rodriguez et al., 2015; O’Leary et al., 2016; Mapelli, 2016), including Pop III stars (Kinugawa et al., 2014). We focus on the current standard picture of massive binary evolution and neglect alternate channels based on chemically homogeneous evolution (Mandel & de Mink, 2016; Marchant et al., 2016). Throughout the paper, quantities refer to BH properties and we use superscript ’’ to refer to properties of the progenitor stars in ambiguous cases (such as the mass).

As in Lamberts et al. (2016), we use the binary stellar evolution code (BSE, Hurley et al., 2002), with modifications for massive binaries. Rapid binary population synthesis codes like BSE are based on prescriptions for the stellar evolution and interactions based on the initial masses and metaliicty of the stars. Such codes allow for a computationally efficient exploration of a wide range of parameters, but can miss subtle effects in the stellar structure and mass transfer (Eldridge & Stanway, 2016). Currently, the main uncertainties in the evolution of massive binary stars are their mass-loss rates (especially for low-metallicity stars), the outcome of the common envelope interactions and the effect of SN kicks. We use the metallicity-dependent prescription for the mass loss rates from Belczynski et al. (2010b) and the remnant mass from Belczynski et al. (2008b). We use the model from Dominik et al. (2013) for the BH initial kicks, which are drawn from a Maxwellian distribution that peaks at 265 km s. The kicks are then reduced according to the amount of material that falls back after core collapse (i.e. by a factor , such that more massive BHs experience smaller kicks). For example, a BH (roughly the median mass of the BHs in our binaries) is given a typical kick of km s. For the common envelope mass transfer, we use the so-called -formalism (Webbink, 1984) using the common envelope efficiency and the envelope binding energy is determined according to the evolutionary stages of the stars. When common envelope occurs during the Hertzprung gap (between core hydrogen burning and shell hydrogen burning), we assume that a stellar merger occurs, as the boundary between the core and envelope is too smooth to stop the inspiral (Ivanova & Taam, 2004, referred to as model B in Belczynski et al., 2007.). Globally, we assume that half of the mass lost by the donor during Roche-lobe overflow is accreted by the secondary.

We create 13 different samples from the BSE model with metallicities logarithmically spaced between and , which are the limits currently allowed in BSE. For each metallicity, we build a statistical sample of binaries, with primary masses between 20 and 120 following a Kroupa IMF (Kroupa, 2001). The secondary masses are drawn assuming that the mass ratios are uniformly distributed, and the initial period distribution is taken to be flat in log-space (Sana et al., 2012). As (single) stars above 120 are likely subject to pair-instability supernova, which does not leave any remnant, we choose not to expand our upper limit. We begin with a thermal distribution for the eccentricities (), which favors high-eccentricity systems. While the initial conditions of the binary evolution (masses and orbital parameters) are somewhat uncertain, de Mink & Belczynski (2015) showed that these uncertainties only affect the final results by at most a factor 2. We set the binary fraction to unity; all numbers can be directly rescaled to lower values. The number of binaries simulated with BSE depends on the metallicity, and is adjusted in order to ensure that our star particles are matched to a smooth distribution of BBHs. In practice, we require at least 3500 BBHs per metallicity bin, which requires an initial sample of approximately 30 000 binaries at but roughly at the highest metallicity.

At most 12 per cent of the initial binaries result in BBHs (for the lowest metallicity). About two thirds of the binaries have companions that are too low mass (initial ) to allow the formation of a second black hole. Among the massive enough binaries, about a third will undergo stellar mergers (see Fig. 1). This is consistent with the estimates of de Mink et al. (2014) that about 10 of massive stars on the main sequence are merger products. Roughly another third of the binaries that would be massive enough to yield a BBH will instead be disrupted during the one of the supernova explosions; the remaining third will result in a BBH. At metallicity beyond , BBH creation is available to less than one per cent of stellar binaries, due to mass loss from stellar winds.

Fig. 2 shows the total mass mass and orbital period at the formation of the BBH according to our binary evolution model. Assuming circular orbits, we indicate the corresponding gravitational wave frequency (twice the value of the orbital frequency). These plots shows that BBH formation is rare at solar metallicity, and limited to , in agreement with (Belczynski et al., 2016; Eldridge & Stanway, 2016). At , BBH formation occurs for about of more massive binaries, but only the lowest metallicity model routinely produces binaries with , as stars at higher metallicities lose most of their mass in strong winds due to higher opacities in their atmosphere.

More specifically, for the lower metallicity models, several channels to BBH formation appear. The systems with the shortest initial orbit (or highest frequency) result from initial stars with , a stellar mass ratio close to unity and a wide enough initial orbit to avoid stellar mergers during the common envelope phase. The narrow strip of binaries with orbits of order of a day typically have and result from more massive stars with initial mass fractions close to unity. These systems have undergone important mass transfer from the primary to the secondary. The other binaries come from initially closer orbits, with stars with and . Systems with periods of order a month or larger typically come from equal mass binaries at higher masses.

The main output of the BPS is a list of 3500 “sample” binaries for each metallicity. For each binary, we record its initial stellar masses and orbital properties, the formation time of the BBH with respect to the formation of the progenitors and its masses and orbital properties. In the next section, we explain the gravitational wave properties of a given binary.

2.3 Gravitational wave emission

After the BBH is formed, the binary evolves only via gravitational wave radiation, gradually shortening the orbit. To assess the BBH population of the galaxy at any given point in time , we first need to determine whether a binary with given properties has already merged. The time to coalescence is given by e.g. Maggiore (2008)


where is the reduced mass of the BBH, the initial eccentricity of the binary and a function depending on the orbital evolution of the system ,which is equal to unity for circular binaries. (see Eq. 4.137 in Maggiore (2008).).

If the binary has not yet merged, its semi-major axis and eccentricity , which determine its GW emission, evolve according to (Peters & Mathews, 1963)


For eccentric sources, GWs are emitted over a range of harmonics of the orbital frequency, while only the second harmonic emits for circular orbits. The total energy loss can be significantly higher than for circular orbits, resulting in faster inspirals. Even though many binaries are eccentric at birth, we find that the orbits are close to circular () by the present day. As such, we include eccentricity for the orbital evolution, but compute the present day GW emission assuming circular orbits. For circular orbits, the frequency and characteristic strain of the GW at a given distance of the source are given by


where the chirp mass is given by .

2.4 From a MW model to GW emission

The core of this paper is the combination of a binary population synthesis model (§2.2) with a cosmological model for a MW-mass galaxy (§2.1) and its star formation history as a function of localization and metallicity to determine the GW emission (§2.3). Here we describe how these three aspects are combined.

We first bin the star particles from the m12i simulation into the same 13 metallicity bins as our BPS model. Stars with metallicity below/above the range covered in the BPS model are assigned to the first/last bin. Each star is then randomly assigned a binary from the BPS model at the corresponding metallicity. This effectively associates black hole masses , , an orbital period and eccentricity of the BBH at formation with each star. We also keep track of the formation time of the BBH with respect to the formation of the progenitor stars .

To obtain the present-day distribution of BBH, we determine the time over which to evolve the binary in order to reach , the present-day age of the Universe. We evolve all the BBHs forward during according to Eq. 2 and their initial orbits and masses. For each BBH, we then have the present-day orbital parameters if it has not yet merged (from GW evolution); and spatial localization (from the simulation) and can determine its gravitational wave properties and possibilities for electromagnetic detection.

In this method, we have assumed that each star particle can be uniquely associated with a single BBH. In practice, the star particles have a mass and represent an IMF-averaged group of stars. From the Kroupa IMF, we find that there are about 12 binaries with in such a star particle. For each metallicity, we can determine the expected number of BBHs according to the fraction of massive stars that effectively turn into a BBH (see percentages in Fig. 2). Effectively, we find that star particles with typically form one BBH, stars with lower metallicity create between 1 and 2, and stars with higher metallicity rarely make BBHs. As such, we weigh all of our mock statistical BBHs with the expectation value of the number of BBH associated with an IMF-averaged stellar population of the same mass, age, and metallicity as the simulation star particle. We perform a polynomial fit of the fraction of BBHs as a function of metallicity and use this to extrapolate the probability of producing a BBH for systems beyond our initial metallicity range.6 This weighting is accounted for in all the results presented here.

Although we include BH natal kicks to determine the initial orbits of the BBH in our BPS model, we do not assign these same kicks to our particles and the location of the BH is set by the location of the star particle they stem from. As we discuss in §5, we estimate that this approximation does not impact the main results of our study.

3 Black holes in a MW-mass galaxy

Fig. 3 shows the distribution of formation time and metallicity of all the stars within 300 kpc (top) and the subset of stars that are UMBHs at (middle) as well as systems that have already merged (MBBH, bottom). The star formation extends from very early times to the present day, and is consistent with the global peak of star formation between 10 and 4 Gyr ago (Madau & Dickinson, 2014). While the stellar metallicity globally increases over time, the scatter is important at all ages. The BPS model (§2.2) highlighted that the formation of BBHs very strongly depends on metallicity. While the bulk of recent star formation is significantly super-solar, about a percent of the recent star formation is low enough metallicity to potentially allow BBH formation. Recent star formation in m12i occurs at slightly higher metallicities than observations suggest is the case in the MW (Mackereth et al., 2017), but we emphasize that the drop-off in BBH formation with increasing metallicity is so steep that our predictions would only be marginally impacted by reducing the metallicity of late-time star formation by a factor of two. However, this difference does somewhat depress the number of BBHs in the disk of our model MW.

The middle plot shows the progenitor stars of the UMBH population. The distribution is very limited above , because BBH formation becomes extremely rare (see Fig. 2). This means that in MW-mass galaxies, the bulk of the massive stars formed recently are unable to form BBHs. Star formation over the past 500 Myrs accounts for 3 of the total stellar mass but only of the mass in black hole binaries. Most of the currently present BBHs come from stars formed 8 to 10 Gyrs ago (between and ), when low metallicity star formation was the strongest.

The distribution of progenitor stars of MBBH is shown in the bottom panel. Globally, the systems that have already merged by originate as lower metallicity stars, and thus trace older star formation. This confirms the results presented in Lamberts et al. (2016). The upper limit on the metallicity is around in the BPS model presented here, and the cutoff is even more drastic than for the unmerged systems. This is because the few BBHs that are formed at higher metallicity typically have wide orbits and low chirp masses due to differences that arise during the course of stellar evolution, implying that they will never merge within a Hubble time (see Fig. 2, and also delay time distributions in Belczynski et al., 2008b; Dominik et al., 2012; Lamberts et al., 2016). The exact cutoff is strongly dependent on the details of the BPS model, such as the mass loss through stellar winds, BH natal kicks and the outcome of common envelope evolution, and will likely be revised as our understanding of massive stellar binaries improves. However, an upper limit on the metallicity is likely to persist. Models accounting for strong mixing within the stars (Mandel & de Mink, 2016; Marchant et al., 2016) are also limited to low metallicity progenitors and are likely to stem from the same progenitor population.

Figure 3: Distribution of the formation time and metallicity of all the stars within 300 kpc of the center of the galaxy in the simulation (top), those stars that are progenitors of the currently unmerged-BBHs (middle) and the progenitors of the BBHs that have already merged (bottom). All distributions are normalized to unity.

These plots show that UBBH/MBBH in MW-mass galaxies stem from a different population than most of the stars. As such, we can expect them to also have a different spatial distribution. Fig. 4 shows maps of the distribution of stellar mass (top) and mass in (merged and unmerged) BBH for a MW-mass galaxy viewed face-on (left) and edge-on (right). The stellar mass distribution shows a zoomed-out view of a spiral galaxy with a central bulge and bright disk. The halo, which is more scarcely populated extends to about 40 kpc. Beyond, we find satellites and streams due to infalling satellites. In comparison, the BBH distribution is much less concentrated in the bulge and disk; instead, the halo is much more populated. Additionally, the satellites and streams over-produce BBHs with respect to their stellar content. Satellites are low mass galaxies (see Wetzel et al. 2016 for a full description) with low metallicity star formation. Accordingly, these are prime sites for BBH production and mergers (Lamberts et al., 2016; Elbert et al., 2017; Mapelli et al., 2017; Schneider et al., 2017).

Figure 4: Maps of the projected stellar mass density (top) and the isolated black hole mass density (merged and unmerged) (bottom) for a face-on (left) and edge-on view (right). The two rows show different absolute values (in units of Solar mass per kpc), but the dynamic range of the colors is identical, allowing for comparison between the two spatial distributions. The upper row represents the total stellar mass and is not a mock observation, which would prominently show young stars and include dust attenuation. At smaller scales, galactic components such as the bulge and thin and thick disc are more prominent. The halo, streams and satellite galaxies are overrepresented in the IBH maps owing to their lower metallicities, though we note that supernovae kicks, which we do not include in this work, may smear BBHs relative to the stellar streams and potentially eject them from the satellites. As described in Wetzel et al. (2016), the cosmological simulation produces a realistic distribution of the MW satellites down to a stellar mass of .

Fig. 5 provides a quantitative description of the spatial distribution of the BHs in a MW-mass galaxy. Globally, the BBHs are preferentially on the outskirts of the galaxy. The effect is even stronger for the merged systems. 60 per cent of the systems are located more than 10 kpc from the center of the galaxy, while 95 per cent of the stellar mass is concentrated within that distance. The right panel plots BBH counts as a function of distance above or below the plane of the disk of the galaxy. Roughly 50 per cent of the BBHs are located at least 5 kpc above or below the disk mid-plane. The sudden increases in the cumulative distribution function at radii beyond 100 kpc indicate BBHs in satellite galaxies, which effectively contribute between 5 and 10 per cent of the systems. When considering the mass distribution of systems, we find that the respective dominance of the stellar halo and satellites is even stronger for the systems with total masses above , 10 per cent of which lie beyond 50 kpc in satellites and streams, even though the latter contribute less than 1 per cent of the total stellar mass.

Figure 5: Cumulative distribution function of the number of merged black holes (MBBH, red) and unmerged black holes (UBBH, black) as a function of distance to the galactic center (left) and vertical distance away from the plane of the galaxy (right). The distribution of the stellar mass is shown as a blue dashed line for comparison. For the merged systems, we distinguish different mass ranges with different grey scales. More massive binaries require lower metalicity and are further biased towards the halo and streams and satellites.

Fig. 6 shows the distribution of merged and unmerged systems as a function of total mass of the BBH (left) and metallicity of the progenitor star (right). The total number of systems is shown with the solid lines. Out of a stellar mass of 7.3 within 300 kpc, we find merged BBH systems and unmerged systems. Out of the total current stellar mass, there is a total of in black holes from binary systems (merged or not merged). Assuming that the current stellar mass traces the total star formation we find that of the stellar mass in a MW-mass galaxy turns into a black hole binary. We remind the reader that we do not account for BHs with other companions or that have been kicked out of the initial binary. We provide an estimate of the total number of black holes in a MW-mass galaxy in §5.

Figure 6: Number of merged BBH (MBBH, red) and unmerged BBH (MBBH, black) as a function of total mass (left) and progenitor metallicity (right). The vertical dashed lines show the mean masses and metallicity of the merged and unmerged systems. The lines show the total number of systems while the shaded areas only represent systems formed outside of the main galaxy. The blue line on the right shows the stellar mass for comparison (renormalized to arbitrary units). Above , the vast majority of the stars come from in-situ formation (not shown here). The comparison between stars and BBH progenitors clearly shows that there are more total stars at high metallicity, but more efficient BBH formation at lower metallicity, producing the “peak” at .

The average (median) mass is for the MBBHs. For the UBBHs the average(median) total mass is for UBBHs. This is because most of the BBHs come from progenitors with , where systems with wider obits, which have not merged yet, have more massive black holes (see Fig. 2). Both for the merged and unmerged systems, massive binaries with total masses above make up about 10 of the systems, systems above make up less than of the systems.

The average (median) metallicity of the progenitor stars is for the MBBHs and for the UBBHs. While extremely low metallicity progenitors (Z) are prime candidates for BBH formation (see §2.2), they only contribute about of the systems in MW-mass galaxies. Conversely, progenitors formed with contribute about 30 of the unmerged systems and of the merged systems. In our model, progenitors with supersolar metallicity contribute less than 1 of the formed BBH.

The shaded area in both histograms shows the number of systems formed outside of the main galaxy. Stars are considered to have formed ex situ if their initial distance to the galactic center was more than 30 kpc (where “galactic center” refers to the center of the main progenitor of the z = 0 galaxy, at the time the stars formed). Such systems have either merged into the galaxy by now (and are present in the bulge/halo), are in the process of being merged (and are present in streams) or are still present in the dwarf galaxy where they formed (and are now in satellites). In our MW model, less than 5 per cent of the currently present stellar mass is formed ex situ (Anglés-Alcázar et al., 2017; Sanderson et al., 2017). In comparison, we find that a third of the BHs (merged and not merged) come from ex-situ formation. For total BBH masses above 60, only 40 per cent of the systems were initially formed within the main galaxy. This is because, more than 90 per cent of the stellar mass with originates outside the main galaxy, as does about half of the stellar mass with .

According to our combination of a cosmological model for a MW-mass galaxy and binary population synthesis, we find about 4 million black hole systems stemming from binary interactions, a quarter of which have merged by now. In the following section we describe the prospects of detecting these single and binary black holes with gravitational waves or electromagnetic signatures.

4 Observational perspectives

Regardless of the detection method, the spatial distribution of the BBHs is a key in determining observational survey strategies. Fig. 7 shows the distance to the the sources (MBBHs and UBBHs) with respect to the Sun, i.e. at an arbitrary point along the Solar Circle: on the disk mid-plane and 8 kpc from the center of the galaxy. About merged systems and binaries are present within a kpc. Unfortunately, the majority of both the merged and unmerged systems are located beyond 10 kpc. Most of those are in the galactic halo, off the plane of the disk (see Fig. 5). Finding these sources would require the monitoring of a large fraction of the sky with a deep survey.

Figure 7: Number of predicted merged (red) and unmerged (black) BBH systems within a given distance of a randomly chosen point on the Solar Circle (i.e. 8 kpc from the galaxy center in the plane of the disk.)

Extragalactic mergers of BBH black holes have been detected with gravitational waves (Abbott et al., 2016c) for . The local BBH merger rate inferred by LIGO is 12-213 Gpc yr (The LIGO Scientific Collaboration et al., 2017a). Assuming MW-mass galaxies per comoving Mpc (Baldry et al., 2008), this yields a galactic merger rate yr, assuming all mergers occur in MW-mass galaxies. This makes a BBH merger very unlikely in the MW (Abbott et al., 2016a). At much lower frequencies, LISA is currently the only planned GW detector that may detect BBHs before they merge. LISA will be able to detect cosmological BBHs just a few years before they merge within the LIGO band (Sesana, 2016), and it may also be able to detect BBHs in the MW or its nearby satellites thousands of years before they merge (Belczynski et al., 2010c). Given the negligible likelihood of a merger occurring within the MW itself, we focus here on the latter possibility.

Figure 8: Characteristic strain and GW frequency of the BBHs in our simulation. The black line in the upper right corner shows the current expected sensitivity curve of the LISA mission – only binaries above the line may be detectable. The color shows the number of systems in each bin. Most of the BBH population in the MW will not be detectable in these models, though some of the loudest BBHs in the disk may appear if LISA exceeds expectations on the low frequency end


Fig. 8 shows the gravitational wave frequency and characteristic strain expected from each binary in our model galaxy (Eqs. 4). It shows two distinct populations, with similar frequency distributions but different median values for the strain. The upper population is exclusively composed of binaries within the galaxy, while the lower group is composed of binaries within satellites. The “loudest” binaries ( Hz) have total masses above 40 and are all located within 30 kpc of the Sun, mostly within the disk. The black line in the upper left corner represents the expected noise curve for LISA (Amaro-Seoane et al., 2017; Klein et al., 2016). The model represented here does not predict any BBH to be detectable during the 4 years of planned operations with LISA. In our model, binary properties are drawn probabilistically (in Monte Carlo fashion) for each star particle according to its age and metallicity. We also vary the localization of the Sun along an annulus 8 kpc away from the Galactic center. The global properties of the binaries are identical for different realizations of the model, though certain realizations do predict a binary that could be detectable over a 10 year extended mission. We predict that a LISA-like telescope located on the Solar Circle would detect, at most, 1 BBH; we also find similar results for m12b and m12c.

The population represented in Fig. 8 results from a convolution between the star formation history and spatial structure of the galaxy and massive binary evolution. As shown in Figs. 2 and 6, most of the binaries result from progenitor stars with and had an initial frequency (orbit longer than a day) and quasi-circular orbits. The high frequency binaries as well as most of the low-metallicity binaries shown in Fig. 2 have merged by now. The highest frequency systems () typically have a total mass below 25 and come from progenitors with metallicity .

The detection of stellar BHs without stellar companions, whether they are single BHs or BBHs, is strongly limited as such systems do not emit any electromagnetic radiation. The current strongest indications for the presence of BHs without stellar companions in the MW come from micro-lensing events with inferred lens masses above 5 and no detectable electromagnetic counterpart (Wyrzykowski et al., 2016). Massive lenses lead to months-to-year long lensing events. Distinguishing between binaries with a short orbital period or single sources would be impossible. Without additional information on the distance to the lensed source, it is impossible to break the degeneracy between the mass of the lens and its distance  (Agol et al., 2002). As such, BH candidates can only be identified in a probabilistic sense. Surveys targeting dense regions like the bulge of the galaxy are the best suited for this type of analysis  (Wyrzykowski et al., 2016). When the distance to the source and lens are known through microlensing parallax, the mass of the lens can be determined. Surveys of the Magellanic Clouds or other nearby galaxies are the best suited for this type of search (Wyrzykowski et al., 2011; Mirhosseini & Moniez, 2017), although those events would be very rare due to the small density of the lenses. If our conclusions hold for the global distribution of BH in the MW (see Fig. 1 for all the BHs not considered here), a long-term survey of the Magellanic Clouds or M31 may be the best option to find single BHs.

Faint X-ray and radio emission may be detected if BHs are accreting surrounding gas (Maccarone, 2005). Fig. 9 provides the Bondi-Hoyle accretion rate , with respect to the Eddington accretion rate for all our BHs. The Bondi-Hoyle accretion rate is given by

Figure 9: Cumulative distribution function of the Bondi-Hoyle accretion rate of the ambient gas by merged and unmerged BBHs, normalized to the Eddington accretion rate and the Bondi-Hoyle accretion efficiency . Different colors indicate different distance cuts.

where is the velocity of the BH with respect to the ambient gas and and are the sound speed and density of the accreted gas. The parameter is an “efficiency” factor which represents our ignorance on detailed accretion physics. Observational estimates suggest based on accretion onto neutron stars (Perna et al., 2003). Based on these values of the accretion efficiency, most BH systems are expected to accrete at less than times the Eddington accretion rate. This is because only a small fraction of our BHs are located in the disk and bulge, close enough to dense molecular clouds, where accretion would be large. Moreover, the BHs we consider stem from old stars, which have high proper velocities with respect to the gas. Isolated black holes kicked out of binary systems are likely to have a similarly high velocity, and small accretion rate. The systems with the highest accretions rates are located in a star-forming satellite galaxy. The radiative properties of the accreting black holes dependent on the geometry, cooling properties, and radiative mechanisms of the accretion flow, which are highly uncertain at these low accretion rates, especially in binary BHs. But if is high in these systems, as many as 10 - 30 merged or unmerged BBH systems, mostly moving through dense molecular clouds in the inner disk or star forming satellites, could be accreting at sufficiently large rates to be detectable with all-sky X-ray surveys. Around the Milky Way, the Magellanic Clouds are the only known satellites still actively forming stars.

5 Discussion

5.1 An improved MW model

The unique combination of a high-resolution cosmological simulation of a MW-mass galaxy and a binary population synthesis model allows us to predict the number of unmerged BBHs and merged BBHs in the MW, as well as their mass distribution and localization. From there, we can determine possible detections, with electromagnetic telescopes or gravitational wave detectors.

Our MW model is based on a cosmological simulation of a MW-mass galaxy down to (Wetzel et al., 2016), with a mass resolution of . We find that the number, mass distribution and spatial distribution of the merged BBHs are nearly identical in a simulation of the same galaxy with a factor of 8 fewer particles (mass resolution of ). In contrast, the unmerged BBHs are 15 per cent more numerous in the lower resolution simulation, due to somewhat higher star formation after . This results into more systems from progenitors with , and a slightly lower mean mass for the BBH systems. The spatial distribution of the BHs is the same in both cases, aside from the satellites, which produce fewer BH in the lower-resolution simulation. Altogether, this suggests our results are not strongly sensitive to the resolution of the simulation.

More importantly, the simulation provides a fully cosmological model for the formation of a Milky Way-mass galaxy, self-consistently modeling its stellar disk, satellite population, and stellar halo, each with cosmologically driven formation histories and self-consistent metal enrichment. This is a significant improvement over previous estimates of the population in the MW. The left panel of Fig. 10 compares the star formation rate as a function of time and metallicity in the simulation, versus other models that are used to compute the binary population in the MW. Most models assume a constant star-formation rate at Solar metallicity (straight monochromatic lines in Fig. 10) in the disk and bulge, if the latter is included at all. When present, the halo is assumed to come from a single star burst at (Voss & Tauris, 2003; Belczynski et al., 2010c; Ruiter et al., 2010; Liu & Zhang, 2014). All these models neglect super-solar metallicity star formation. Although this does not impact BBH formation, it may impact the formation of binary neutron stars and/or white dwarfs (Nelemans et al., 2001; Ruiter et al., 2010). Only Mennekens & Vanbeveren (2014) include time-variable, metallicity-dependent star formation (colored squares), although their model globally over-predicts early star formation and underestimates the variation in the mean metallicity over time compared to observational constraints in MW-mass galaxies (e.g. Behroozi et al. (2013). None of the models in the literature include scatter in the metallicity at a given time, while the dashed lines (and also Fig. 3) and observations (Garcia Perez et al., 2017) show that there is a typical range of about a dex between the highest and lowest metallicity at a given time. This is crucial as the lowest metallicity systems most significantly contribute to the formation of BBHs.

Models in the the literature also neglect the crucial impact of galactic mergers. In Fig. 6 we have shown that about a third of the BBHs present in the MW come from progenitors formed in a satellite galaxy. At the highest masses, more than 60 per cent of the BBHs comes from ex situ formation. This is somewhat different from the findings of Chakrabarti et al. (2017) who predict that most of the progenitors of BBH mergers come from the outer disk of massive galaxies. Neglecting galactic mergers underestimates the global BBH population, specifically high mass binaries and binaries in the galactic halo.

Belczynski et al. (2010c) assume stars form at solar metallicity over the past 10 Gyrs. In contrast, the simulation shows that the majority of stars formed during the last 7 Gyrs have a supersolar metallicity, and are unavailable to BBH formation. On the other hand, they neglect most of the star formation between , where we predict the majority of the BBHs originate (see Fig. 3). This may explain why their model B (which is similar to our BPS model) predicts six times fewer binaries than our model. They also neglect the contribution from ex-situ stars and effectively miss 40 per cent of the binaries. Although our overall number of binaries is much higher, we do not predict BBHs to be detected with LISA, whereas model A in Belczynski et al. (2010c) predicts a few. This is because most of our BBHs are too far from the Sun, in the outer disk and mostly in the halo.

Figure 10: Comparison of the star formation rate as a function of metallicity and time used in the m12i simulation compared to previous models for the MW used to derive BH populations in the Galaxy (left) and other simulations of MW-type galaxies (right) from the FIRE project. In both plots, the circles show the median values in the m12i simulation, color-coded by their relative star formation rate and the black lines show the 95 scatter in metallicity for stars forming at the same time. On the left plot, monochromatic lines show models that assume a constant star formation rate (e.g. Voss & Tauris (2003); Ruiter et al. (2010); Belczynski et al. (2010c); Liu & Zhang (2014)), while colored points are shaded by the relative star formation rate at a given time (e.g. Mennekens & Vanbeveren 2014). Burst of star formation are shown with stars symbols. On the right plot, we show the m12b simulation with inverted triangles and m12c simulation with triangles.

We also show the star formation history of two other simulations of MW-mass galaxies run with the same resolution and physics (right panel of Fig. 10). The m12b simulation has a 27 per cent higher present-day stellar mass than m12i, but produces fewer BHs. This is because most of its star formation is at too high metallicity to significantly produce BBHs (see inverted triangles in Fig. 10). m12c, meanwhile, has 7 per cent less present-day stellar mass, but also produces 15 (13) per cent fewer merged (unmerged) BBHs. This is because its low metallicity star formation around is lower than in our reference simulation. The global properties (mass, progenitor metallicity and spatial distribution) of the BHs are similar in the three simulations. The comparison with different MW-mass galaxies highlights the robustness of our calculation, and indicates that our simulation is a reliable representation of the Milky Way galaxy, even though it is not an exact reproduction. The comparison also shows that although the present-day stellar mass provides a first indication of the BBH content of a galaxy, the details of its star formation history, and more specifically its metallicity can have an important effect.

Our host galaxy has a higher present-day star formation rate than the MW. As recent star formation has a limited contribution to BBH formation (see Fig. 3), we do not expect this to influence the quality of our Milky Way model. Although the global properties of the satellites are consistent with observations (Wetzel et al., 2016), the simulation does not show a massive satellite like the Large Magellanic Cloud. The latter still presents star formation with (Harris & Zaritsky, 2009) and is thus a prime candidate for BBH formation. As such, our estimate of the number of UBBH and MBBH within 300 kpc is likely a lower limit.

5.2 Pathways towards detections

The accurate spatial model in our fully cosmological simulation, with respect to simplified disk and halo models, allows us to provide some predictions of the detectability of the sources with electromagnetic or gravitational signatures. Although BH natal kicks are included in our binary evolution model, we do not model the impact of BH kicks on the location and proper motion of the BHs in the galaxy simulation. Black holes kicks are assumed to be smaller than neutron star kicks, due to material falling back. As we are specifically focusing on BBHs surviving SN kicks, we are biased toward the BBHs with small natal kicks (otherwise the binary would have been disrupted). Additionally, most of our systems are found in the stellar halo, which is dominated by random motions, where BH natal kicks are likely to have a limited impact of the statistical properties of the total BBH distribution. Systems formed in the disk may be kicked into the halo, especially if the kicks occur before most the mass of the host galaxy is accreted (). For systems formed in satellite galaxies, kicks may be powerful enough to escape the shallow potential well of the dwarf galaxy but it is unlikely that they would be able to escape the potential of the main host. As a result, the distribution of ex-situ systems may be more randomized, and BBHs may not closely trace stellar streams and satellites. Overall, we expect most of the merged and unmerged BBHs to be in the galactic halo, and we argue that this result is robust to our assumption that BBHs trace the location of their parent star particles.

The prospects for detecting black holes without stellar companions remain limited. Still, with a refined model of the star formation history in the MW, non-detections with X-ray, radio and microlensing surveys will put constraints on binary evolution models. Conversely, a much higher number of detected systems will be equally constraining, and may inform us of additional BBH formation channels (Marchant et al., 2016; Rodriguez et al., 2015). Our understanding of massive binary evolution is still very uncertain, especially the modeling of stellar mass loss through different evolutionary phases, the properties of common-envelope mass transfer, and supernova kicks. Specifically, we have assumed that common envelope interactions during the Hertzprung gap result in stellar mergers, which may underestimate the BBH production rate. Conversely, we assume the BH natal kicks are reduced with respect to neutron star natal kicks, which may overestimate the number of systems surviving both supernova explosions. As our work highlights the importance of improved galactic models, we do not explore the wide parameter space of current binary models. The different BBH merger rates predicted by Mapelli et al. (2017) show how different binary evolution models will eventually be ruled in or out by observational data, provided appropriate models are used for the star formation.

Based on our assumptions for binary evolution, we predict that LISA will detect zero or one BBH within the MW halo. While certain aspects of binary evolution are likely to be revised, we emphasize that the BBH predictions of our model are somewhat on the high end (see a discussion in Lamberts et al., 2016) and that fine-tunning it in order to produce more galactic BBHs may overproduce the observed BH merger rate (The LIGO Scientific Collaboration et al., 2017a). Unless a mechanism allows for the formation of BBHs at close-to solar metallicity, BBHs are strongly biased towards the galactic halo, which strongly reduces their GW signal on Earth relative to disk populations. Given the frequency and strain distribution of our binaries, we predict a handful of detections if LISA exceeds its expectations at its low-frequency end. Even a single detection with accurate localization will provide information on the binary evolution mechanism and the conditions of its formation.

We find that about a million binary black holes have merged in the MW by the present day. Focusing on the mergers within the last 2 Gyrs, we find a merger rate of yr. Our model, which assumes a binary fraction of unity, is on the higher end of the measured rate. We find the mean mass of the systems to be 26.5, which is in line with the announced detections (Abbott et al., 2016b, 2017a, 2017b; The LIGO Scientific Collaboration et al., 2017b). We find that sources with represent 8 per cent of the merged systems, meaning there could be roughly 80,000 such BBHs in our Galaxy. The progenitors of these systems likely formed in a satellite galaxy and are now present in the halo. Similarly to Lamberts et al. (2016), we find that mergers in MW-mass galaxies primarily stem from progenitors with .

Initial studies predicted roughly 10 BH in the Milky Way (single and binary), based on stellar evolution models (Shapiro & Teukolsky, 1983; van den Heuvel, 1992) and yields from massive stars (Samland, 1998). Based on metallicity-dependent star formation models of galaxies of all masses, Elbert et al. (2017) predict black holes in the MW, with 10 per cent of them above 30 . Our model predicts 3 million binary black holes, and about a million merged systems. Our calculation only accounts for binary black holes from field binaries, and we only track systems where both stars form a BH, where the binary does not undergo a stellar merger, and where the binary survives the natal kicks. As such, we do not account for black holes with lower mass companions (white dwarfs, neutron stars, or low-mass stars). Thousands of black holes with stellar companions will likely by detected by the astrometric GAIA satellite (Mashian & Loeb, 2017; Breivik et al., 2017) and possibly their H emission (Casares, 2017). Given that even at the lowest metallicity, only 12 per cent of the massive binaries turn unto a BH binary, there could be 10 times more black holes that have been kicked out of a binary. As such, the first electromagnetic detection of isolated black holes in the MW will most likely be a single black hole. A more accurate determination of their masses, proper motions, spatial distribution and ultimately, detectability, is left for a further study.

6 Conclusions

In this paper we provide the first estimates of the isolated black hole population resulting from binary black hole systems. It is based on the combination of a high resolution cosmological simulation of a MW-mass galaxy and a binary population synthesis model. The simulation provides a physically-motivated, metallicity-dependent star formation history as well as a complete description of the galactic morphology and merger activity over time. The simulation models both resolved and unresolved turbulent metal diffusion, which provides a realistic metal distribution, rather than just the average value (Ma et al., 2017; Escala et al., 2017). The stellar metallicity is a key parameter for massive binary evolution. Using a standard binary evolution model, we compute a metallicity-dependent library of binary black holes for 13 metallicities between and , then match them to the stars in the simulation. This provides a self-consistent distribution of black holes within 300 kpc of the center of the galaxy, including their localization, masses, orbital properties and the properties of their stellar progenitors. The main properties of these binaries are summarized below:

  • We find that a million binary black hole binaries have already merged in the MW and 3 million systems are still in binary black holes. Our Milky-Way like galaxy has turned 0.15 per cent of its stellar mass into binary black holes, including the ones that have merged already.

  • The mean progenitor metallicity of the merged (unmerged) systems is  , and only 1 per cent of the binaries come from super-solar metallicity progenitors. This means that most of the stellar mass in MW-mass galaxies is effectively unavailable for the formation of black hole binaries. The binary systems thus strongly trace star formation around while the already-merged systems mostly trace star formation from

  • The strong dependence on low metallicity star formation results in half of the binaries (merged or not) being located beyond 10 kpc of the galactic center. In comparison, 90 per cent of the stellar mass is located within 10 kpc. The galactic halo, streams and satellite galaxies are rich in BBHs.

  • We find about 50,000 merged binaries with masses comparable to the mergers detected with LIGO. Consistent with Lamberts et al. (2016), we find that these systems were typically formed outside of the MW and are now in the halo or still in their host satellite galaxy. They stem from progenitors with metallicity below .

  • The detection of merged and unmerged binary black holes without stellar companions will remain a challenge for the foreseeable future. Most of the binaries have an orbit that is too wide and/or are too distant for LISA to detect their gravitational wave emission. Depending on the radiation efficiency and accretion mechanism, the accretion of surrounding gas could lead to detections with all-sky radio or X-ray surveys like SKA or eROSITA. A few systems may be detected by microlensing surveys. In any case, observational constraints will lead to a better understanding of massive binary evolution, which is still poorly understood.

  • About a third of the binary black holes were not initially formed in the galaxy, and have been brought in by the accretion of satellite dwarf galaxies. 60 per cent of the systems with total mass above 60 were formed ex situ. This highlights the importance of accounting for galactic mergers when predicting black hole populations and merger rates in MW-mass galaxies.

  • The mean total mass is 26 for the merged systems and a total mean mass of 28 for the binary systems. Roughly 10 per cent of the systems have total masses above 50 and about one per cent of the systems have a mass above 80 .

This paper provides the first estimate of the population of binary black holes and their merger remnants in a MW-mass galaxy that incorporates a realistic star-formation history and galactic halo structure and merger history based on a hydrodynamic simulation. This allows us to reliably estimate low metallicity star formation and localize binary black holes. Based on this precise galactic model, future detections of such systems in our vicinity will then allow us to constrain key parameters of massive stellar binary evolution.


Astrid Lamberts would like to thank V. Ravi, H. Vedantham, M. Heida, C. Henderson, Y. Shvarzvald, S. Novati and S. Taylor for discussions about observational implications of this work. Numerical calculations were run on the Caltech compute cluster “Wheeler,” allocations from XSEDE TG-AST130039 and PRAC NSF.1713353 supported by the NSF, and NASA HEC SMD-16-7592.

Support for AL and PFH was provided by an Alfred P. Sloan Research Fellowship, NASA ATP Grant NNX14AH35G, and NSF Collaborative Research Grant 1715847 and CAREER grant 1455342. Support for SGK was provided by NASA through Einstein Postdoctoral Fellowship grant number PF5-160136 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. EQ was supported in part by NSF grant AST-1715070 and a Simons Investigator Award from the Simons Foundation. JSB was supported by NSF grant AST-1518291 and by NASA through HST theory grants (programs AR-13921, AR-13888, and AR-14282.001) awarded by STScI, which is op- erated by the Association of Universities for Research in Astron- omy (AURA), Inc., under NASA contract NAS5-26555. CAFG was supported by NSF through grants AST-1412836, AST-1517491, AST-1715216 and CAREER award AST-1652522 and by NASA trough grant NXX-15AB22G. AW was supported by NASA through grants HST-GO-14734 and HST-AR-15057 from STScI. DK acknowledges support from NSF grants AST-1412153 and AST-1715101 and the Cottrell Scholar Award from the Research Corporation for Science Advancement. RES was supported by an NSF Astronomy & Astrophysics Postdoctoral Fellowship under grant AST-1400989. This study was initiated during K.


  1. pubyear: 2018
  2. pagerange: Predicting the binary black hole population of the Milky Way with cosmological simulationsPredicting the binary black hole population of the Milky Way with cosmological simulations
  5. Whenever we refer to the simulation, we use the words star, particle and star particle interchangeably.
  6. Due to the strong dependence on metallicity (as shown in Fig. 2), this extrapolation produces a negligible number of additional BBHs and therefore has a minimal impact on our results.


  1. Abbott B. P., et al., 2016a, preprint, (arXiv:1602.03842)
  2. Abbott B. P., et al., 2016b, Physical Review X, 6, 041015
  3. Abbott B. P., et al., 2016c, Phys. Rev. Lett., 116, 061102
  4. Abbott B. P., et al., 2017a, Phys. Rev. Lett., 118, 221101
  5. Abbott B. P., et al., 2017b, Phys. Rev. Lett., 119, 141101
  6. Agol E., Kamionkowski M., 2002, \mnras, 334, 553
  7. Agol E., Kamionkowski M., Koopmans L. V. E., Blandford R. D., 2002, \apjl, 576, L131
  8. Amaro-Seoane P., et al., 2017, preprint, (arXiv:1702.00786)
  9. Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2017, \mnras, 470, 4698
  10. Badenes C., et al., 2017, preprint, (arXiv:1711.00660)
  11. Baldry I. K., Glazebrook K., Driver S. P., 2008, \mnras, 388, 945
  12. Behroozi P. S., Wechsler R. H., Conroy C., 2013, \apj, 770, 57
  13. Belczynski K., Kalogera V., Bulik T., 2002, \apj, 572, 407
  14. Belczynski K., Taam R. E., Kalogera V., Rasio F. A., Bulik T., 2007, \apj, 662, 504
  15. Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008a, \apjs, 174, 223
  16. Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008b, \apjs, 174, 223
  17. Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010a, \apj, 714, 1217
  18. Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010b, \apj, 714, 1217
  19. Belczynski K., Benacquista M., Bulik T., 2010c, \apj, 725, 816
  20. Belczynski K., Holz D. E., Bulik T., O’Shaughnessy R., 2016, preprint, (arXiv:1602.04531)
  21. Bonaca A., Conroy C., Wetzel A., Hopkins P. F., Kereš D., 2017, \apj, 845, 101
  22. Breivik K., Chatterjee S., Larson S. L., 2017, preprint, (arXiv:1710.04657)
  23. Casares J., 2017, preprint, (arXiv:1711.03553)
  24. Casares J., Jonker P. G., Israelian G., 2017, preprint, (arXiv:1701.07450)
  25. Chakrabarti S., Chang P., O’Shaughnessy R., Brooks A., Shen S., Bellovary J., Gladysz W., Belczynski C., 2017, preprint, (arXiv:1710.09407)
  26. Christian P., Loeb A., 2017, \mnras, 469, 930
  27. Corbel S., et al., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 53
  28. Corral-Santana J. M., Casares J., Muñoz-Darias T., Bauer F. E., Martínez-Pais I. G., Russell D. M., 2016, \aap, 587, A61
  29. Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, \apj, 759, 52
  30. Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2013, \apj, 779, 72
  31. Elbert O. D., Bullock J. S., Kaplinghat M., 2017, preprint, (arXiv:1703.02551)
  32. Eldridge J. J., Stanway E. R., 2016, \mnras, 462, 3302
  33. Escala I., et al., 2017, preprint, (arXiv:1710.06533)
  34. Fender R. P., Maccarone T. J., Heywood I., 2013, \mnras, 430, 1538
  35. Fryer C. L., 1999, \apj, 522, 413
  36. Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, \apj, 749, 91
  37. Garcia Perez A. E., et al., 2017, preprint, (arXiv:1712.01297)
  38. Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, \mnras, 448, 3248
  39. Harris J., Zaritsky D., 2009, \aj, 138, 1243
  40. Hopkins P. F., 2015, \mnras, 450, 53
  41. Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, \mnras, 445, 581
  42. Hopkins P. F., et al., 2017, preprint, (arXiv:1702.06148)
  43. Hurley J. R., Tout C. A., Pols O. R., 2002, \mnras, 329, 897
  44. Ioka K., Matsumoto T., Teraki Y., Kashiyama K., Murase K., 2017, \mnras, 470, 3332
  45. Ivanova N., Taam R. E., 2004, \apj, 601, 1058
  46. Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, \apjs, 125, 439
  47. Izzard R. G., Tout C. A., Karakas A. I., Pols O. R., 2004, \mnras, 350, 407
  48. Janka H.-T., 2013, \mnras, 434, 1355
  49. Kinugawa T., Inayoshi K., Hotokezaka K., Nakauchi D., Nakamura T., 2014, \mnras, 442, 2963
  50. Klein A., et al., 2016, \prd, 93, 024003
  51. Kroupa P., 2001, \mnras, 322, 231
  52. Lamberts A., Garrison-Kimmel S., Clausen D. R., Hopkins P. F., 2016, \mnras, 463, L31
  53. Liu J., Zhang Y., 2014, \pasp, 126, 211
  54. Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, \mnras, 456, 2140
  55. Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, \mnras, 467, 2430
  56. Maccarone T. J., 2005, \mnras, 360, L30
  57. Mackereth J. T., et al., 2017, \mnras, 471, 3057
  58. Madau P., Dickinson M., 2014, \araa, 52, 415
  59. Maggiore M., 2008, Gravitational Waves. Oxford University Press
  60. Mandel I., de Mink S. E., 2016, \mnras, 458, 2634
  61. Mapelli M., 2016, \mnras, 459, 3432
  62. Mapelli M., Giacobbo N., Ripamonti E., Spera M., 2017, \mnras, 472, 2422
  63. Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, \aap, 588, A50
  64. Marigo P., 2001, \aap, 370, 194
  65. Mashian N., Loeb A., 2017, \mnras, 470, 2611
  66. Mennekens N., Vanbeveren D., 2014, \aap, 564, A134
  67. Mirhosseini A., Moniez M., 2017, preprint, (arXiv:1711.10898)
  68. Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, \mnras, 454, 2691
  69. Nelemans G., Yungelson L. R., Portegies Zwart S. F., 2001, \aap, 375, 890
  70. Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nuclear Physics A, 777, 424
  71. O’Leary R. M., Meiron Y., Kocsis B., 2016, \apjl, 824, L12
  72. Perna R., Narayan R., Rybicki G., Stella L., Treves A., 2003, \apj, 594, 936
  73. Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  74. Planck Collaboration et al., 2016, \aap, 594, A13
  75. Postnov K. A., Yungelson L. R., 2014, Living Reviews in Relativity, 17, 3
  76. Rodriguez C. L., Morscher M., Pattabiraman B., Chatterjee S., Haster C.-J., Rasio F. A., 2015, Physical Review Letters, 115, 051101
  77. Ruiter A. J., Belczynski K., Benacquista M., Larson S. L., Williams G., 2010, \apj, 717, 1006
  78. Samland M., 1998, \apj, 496, 155
  79. Sana H., et al., 2012, Science, 337, 444
  80. Sanderson R. E., et al., 2017, preprint, (arXiv:1712.05808)
  81. Schneider R., Graziani L., Marassi S., Spera M., Mapelli M., Alparone M., de Bennassuti M., 2017, \mnras, 471, L105
  82. Sesana A., 2016, Physical Review Letters, 116, 231102
  83. Seto N., 2016, \mnras, 460, L1
  84. Shapiro S. L., Teukolsky S. A., 1983, Black holes, white dwarfs, and neutron stars: The physics of compact objects
  85. Stevenson S., Ohme F., Fairhurst S., 2015, \apj, 810, 58
  86. Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, \mnras, 471, 144
  87. Sukhbold T., Ertl T., Woosley S. E., Brown J. M., Janka H.-T., 2016, \apj, 821, 38
  88. The LIGO Scientific Collaboration et al., 2017a, preprint, (arXiv:1711.05578)
  89. The LIGO Scientific Collaboration et al., 2017b, preprint, (arXiv:1711.05578)
  90. Voss R., Tauris T. M., 2003, \mnras, 342, 1169
  91. Webbink R. F., 1984, \apj, 277, 355
  92. Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, \apjl, 827, L23
  93. Wyrzykowski L., et al., 2011, \mnras, 416, 2949
  94. Wyrzykowski Ł., et al., 2016, \mnras, 458, 3012
  95. de Mink S. E., Belczynski K., 2015, \apj, 814, 58
  96. de Mink S. E., Sana H., Langer N., Izzard R. G., Schneider F. R. N., 2014, \apj, 782, 7
  97. van den Heuvel E. P. J., 1992, Technical report, Endpoints of stellar evolution: The incidence of stellar mass black holes in the galaxy
  98. van den Hoek L. B., Groenewegen M. A. T., 1997, \aaps, 123, 305
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