© ESO, 2016

1. Introduction

Stellar binaries that consist of two stellar remnants so close to each other that they can merge within the Hubble time have been known for a long time. Merging double white dwarf systems are thought to provide one potential channel to produce so-called Type Ia supernovae, which are the main producers of iron and through which the accelerated expansion of the Universe was discovered. Various double neutron star systems have also been found, most importantly the Hulse-Taylor system (Hulse & Taylor 1975; Weisberg et al. 2010) and more recently the double pulsar PSR J0737-3039 (Burgay et al. 2003; Kramer et al. 2006), whose orbital decay has unambiguously confirmed Einstein’s gravitational wave prediction. Double black hole binaries have not yet been detected, which is obviously difficult since they do not emit electromagnetic radiation. However, because they are the most massive of the double compact binaries, their gravitational-wave radiation would be much stronger than that of white dwarf or neutron star systems.

The evolution of a binary system from the initial stage of two orbiting main-sequence stars to the double compact binary stage is in most scenarios believed to require at least one so-called common-envelope phase. This is essential because stars tend to expand dramatically after their main-sequence evolution, and as an early merging of the binary is to be avoided, a large orbital separation is required to accommodate this. The common-envelope stage is then necessary to shrink the system to a compact final configuration. While the physics of this process is not yet well understood (Ivanova et al. 2013), some observational constraints exist for low-mass binaries (e.g. Han et al. 2003; Zorotovic et al. 2011). For massive binaries, observational evidence is much scarcer, and the common-envelope process is much less well understood theoretically (Taam & Sandquist 2000). For double neutron star systems, the increased number of observed systems is beginning to provide some constraints (Kalogera et al. 2004). For double black hole systems, however, the uncertainties are so large that the predicted rates of black hole binaries from the common-envelope channel are uncertain by several orders of magnitude (Abadie et al. 2010).

A completely different route towards double compact binaries is explored in this paper. It does not involve a common-envelope phase and may work only for very massive stars: it involves the chemically homogeneous evolution of rapidly rotating stars in tidally locked binaries.

The stabilizing effects of entropy and composition gradients usually prevent the efficient mixing of gas once composition gradients have been established. However, the more massive a star is, the less this is the case because of the increased importance of radiation pressure. It has been found that mixing induced by rapid rotation can, in principle, keep massive stars chemically homogeneous throughout core hydrogen burning (Maeder 1987; Langer 1992; Heger & Langer 2000). Detailed studies of this type of evolution through large grids of stellar models (Yoon & Langer 2005; Woosley & Heger 2006; Brott et al. 2011; Köhler et al. 2015; Szécsi et al. 2015) have shown that this only works at low metallicity where strong angular-momentum loss that is due to stellar winds can be avoided, such that the stars remain fully mixed until their core hydrogen is exhausted. Because in this case, rapidly spinning iron cores are produced at the end of the stars’ lives, such single-star models have been suggested as progenitors of long-duration gamma-ray bursts (LGRBs; Woosley & Heger 2006; Yoon et al. 2006).

Chemically homogeneously evolving stars avoid the strong post-main-sequence expansion because they do not maintain a massive hydrogen-rich envelope. de Mink et al. (2009) therefore suggested that massive close binaries, where rapid rotation and thus chemically homogeneous evolution can be enforced through the tidal interaction of both stars, could evolve towards black holes without ever encountering contact or mass transfer (see also Mandel & de Mink 2016; Song et al. 2016).

We here explore the evolution of close binaries with component masses above ~20 M_⊙ by computing large grids of detailed binary evolution models. For this purpose, we use the publicly available code MESA, which we extended to allow the inclusion of contact binaries with mass-ratios close to one (Sect. 2). Against our initial expectation, we find that contact-free evolution occurs only very rarely. Instead, when computing the evolution of massive overcontact binaries (MOBs), we find many systems that avoid merging during core hydrogen burning. We compute the final configurations of these binaries in Sect. 3, including the black hole masses, separations, and their mass ratios. In Sect. 4 we discuss the predicted black hole Kerr parameters, potential explosive mass loss and momentum kicks, and the connection of the MOB scenario with long-duration gamma-ray bursts (GRBs) and pair-instability supernovae (PISNe). We discuss event rates and potential LIGO detection rates in Sect. 5 before giving our conclusions in Sect. 6.

2. Methods

We here provide the first detailed binary stellar evolution models that are followed until the double black hole stage. To obtain these, we applied the MESA code (Paxton et al. 2015, 2013, 2011), which now includes all the physics required for such calculations, in particular, tidal interactions and differential rotation. We computed about2000 detailed binary-evolution sequences in six model grids for different initial metallicities and mass ratios, thereby achieving a complete coverage of the relevant parameter space. Our model calculations include the overcontact phase, which occurs in the closest simulated binaries. This constitutes the main channel for providing massive close black hole binaries.

2.1. Physics implemented in MESA and initial parameters

To model the evolution of our systems, we used the binary module in version r8118 of the MESA code¹. Opacities were calculated using CO-enhanced opacity tables from the OPAL project (Iglesias & Rogers 1996), computed using solar-scaled abundances based on Grevesse et al. (1996). Convection was modelled using the standard mixing-length theory (Böhm-Vitense 1958) with a mixing-length parameter α = 1.5, adopting the Ledoux criterion. Semiconvection was modelled according to Langer et al. (1983) with an efficiency parameter α_sc = 1.0. We included convective core overshooting during core hydrogen burning following Brott et al. (2011). The effect of the centrifugal force was implemented as in Heger & Langer (2000). Composition- and angular-momentum transport due to rotation includes the effects of Eddington-Sweet circulation, secular and dynamical shear instabilities, and the GSF instability, with an efficiency factor f_c = 1/30. This corresponds to the calibrations of the mixing efficiency in stellar models based on the VLT FLAMES Survey of Massive Stars (Brott et al. 2011, and references therein). We included the effect of magnetic fields on the transport of angular momentum as in Petrovic et al. (2005). Tidal effects were implemented as in Hurley et al. (2002) and Detmers et al. (2008) for the case of stars with a radiative envelope. We are not interested in following the nucleosynthesis in detail, therefore we used the simple networks provided with MESA basic.net for H and He burning, co_burn.net for C and O burning, and approx21.net for later phases.

Our implementation of stellar winds follows that of Yoon et al. (2006), with mass-loss rates for hydrogen-rich stars (with a surface helium abundance Y_s< 0.4) computed as in Vink et al. (2001), while for hydrogen-poor stars (Y_s> 0.7) we used the recipe of Hamann et al. (1995) multiplied by a factor of one tenth. In the range 0.4 <Y_s< 0.7, the rate was interpolated between the two. For both rates we used a metallicity scaling of (Z/Z_⊙)^0.85. We also included the enhancement of winds through rotation as in Heger & Langer (2000), and, when the rotation rate exceeded a given threshold Ω/Ω_crit> 0.98, we implicitly computed the mass-loss rate required for the rotation rate to remain just below this value.

Whenever one component in the system attempted to overflow its Roche lobe, we implicitly computed the mass-transfer rate necessary for it to remain just within the Roche volume (computed as in Eggleton 1983). The treatment of mass transfer in overcontact systems is described in the following section.

We considered four different metallicities, Z_⊙/ 4, Z_⊙/ 10, Z_⊙/ 20, and Z_⊙/ 50, with Z_⊙ = 0.017 as in Grevesse et al. (1996). The helium abundance was set in such a way that it increased linearly from its primordial value Y = 0.2477 (Peimbert et al. 2007) at Z = 0 to Y = 0.28 at Z = Z_⊙. For all metallicities, we computed grids for mass-ratios q_i = M₂/M₁ = 1; for Z_⊙/ 50, we also computed grids at q_i = 0.9,0.8. The initial orbital periods were chosen from the range (P_i = 0.5−3.0) days, with an interval of 0.1 days, while the initial primary masses were taken from the range (log M₁/M_⊙ = 1.4−2.7) in intervals of 0.1 dex.

2.2. Computing massive overcontact systems

Very close binaries may evolve into contact where both binary components fill and even overfill their Roche volumes. The evolution during the overcontact phase differs from a classical common-envelope phase because co-rotation can, in principle, be maintained as long as material does not overflow the L2 point. This means that a spiral-in that is due to viscous drag can be avoided, resulting in a stable system evolving on a nuclear timescale.

As a simple approximation to the modelling of overcontact systems with a 1D code, we considered the photosphere of such an object to lie on a Roche equipotential Φ and divided it into two distinct volumes for each star, V₁(Φ) and V₂(Φ), separated by a plane crossing through L₁ perpendicular to the line joining both stars. We then associated a volume-equivalent radius with each of these, R₁(Φ) and R₂(Φ), with the radii corresponding to the potential at L₁ as the Roche-lobe radius R_RL of each. When both stars overfill their Roche-lobes, the amount of overflow of one component is a function of the mass ratio, and the amount of overflow of the other component is $\begin{array}{ccc}\frac{{\mathit{R}}_{\mathrm{2}}\mathrm{\left(}\mathrm{\Phi }\mathrm{\right)}\mathrm{-}{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{2}}}{{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{2}}}\mathrm{=}\mathit{F}\left(\mathit{q,}\frac{{\mathit{R}}_{\mathrm{1}}\mathrm{\left(}\mathrm{\Phi }\mathrm{\right)}\mathrm{-}{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{1}}}{{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{1}}}\right)\mathit{,}& & \end{array}$(1)where the function F(q,x) must satisfy the conditions F(q,0) = 0 and F(1,x) = x. By numerically integrating V₁(Φ) and V₂(Φ) for different q ratios and potential values between those at L₁ and L₂, we found the fit F(q,x) = q^-0.52x, with the Roche-lobe radius computed as in Eggleton (1983), that approximates F(q,x) with a few percent error in the range 0.1 ≤ q ≤ 1. During an overcontact phase, mass transfer is then adjusted in such a way that the amount of overflow from each component satisfies this relationship.

After both stars overflow past the outer Lagrangian point, we expect the system to merge rapidly, either as a result of mass loss from L₂ carrying a high specific angular momentum, or in consequence of a spiral-in due to the loss of co-rotation. We found that, given the volume equivalent radii for the potential at the L2 point R_L2,i, the amount of overflow of the least massive star in a system that reaches L₂ approximately satisfies the relationship $\begin{array}{ccc}\frac{{\mathit{R}}_{\mathrm{L}\mathrm{2}\mathit{,}\mathrm{2}}\mathrm{-}{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{2}}}{{\mathit{R}}_{\mathrm{RL}\mathit{,}\mathrm{2}}}\mathrm{=}\mathrm{0.299}{\tan }^{-1}\mathrm{(}\mathrm{1.84}{\mathit{q}}^{\mathrm{0.397}}\mathrm{)}& & \end{array}$(2)with an error smaller than 2% in the range 0.02 ≤ q ≤ 1. At q = 1 this means that a star needs to expand by up to 1.32 times its Roche-lobe radius before reaching L2, which leaves a significant amount of space for a binary to survive through an overcontact phase. We note that we ignored the effects of energy and element transfer through the shared envelope for the moment. However, since the systems we model in this work undergo contact as rather unevolved stars with mass ratios not far from one, we expect these effects to be of minor importance for our present study.

It is worth mentioning that VFTS352 is a massive (~30 M_⊙ + 30 M_⊙), short-period (P_orb = 1.12 d) overcontact binary that evolves on the nuclear timescale. It has a mass ratio of q = 1.008 and is thought to undergo chemically homogeneous evolution (Almeida et al. 2015). This system therefore provides direct support for the MOB scenario and for our treatment of this phase.

3. Results

Before examining the wider parameter space, we provide an example of a typical MOB evolution in the following section.

3.1. Exemplary MOB evolution

We show in Fig. 1 the evolutionary tracks of the two components of a 79 + 64 M_⊙ binary at Z = Z_⊙/ 50 with an initial period of 1.1 d from the zero-age main sequence until core helium exhaustion. This system enters an overcontact phase early during core hydrogen burning, during which it swaps mass back and forth several times. Each time, the mass-ratio becomes closer to one, so that eventually contact is avoided and the system evolves as a detached binary with two stars of about 71 M_⊙ from this time on.

Core hydrogen burning ends at an effective temperature of log T_eff ≃ 4.9, after which both stars contract to even smaller radii. As a result of the different initial evolution, the stars do not evolve exactly synchronously during the late evolutionary stages.

Figure 2 illustrates the different evolutionary stages of a typical binary in our grid from the zero-age main sequence (ZAMS) to the black-hole merger stage. While chemically homogeneous evolution is maintained throughout core hydrogen burning, mass is transferred back and forth between the two binary components in a succession of contact stages, which eventually leads to a mass ratio very close to one. During the main-sequence phase and the post-core-hydrogen-burning phase, where both stars are compact detached helium stars, stellar-wind mass loss leads to a widening of the orbit. Because the winds are metallicity dependent (Mokiem et al. 2007; Vink & de Koter 2005), metal-rich systems are found to widen strongly, which limits the tidal interaction so that the homogeneous evolution may end. Furthermore, if the orbital period is too long, any BH+BH binary that may be produced will not merge in a Hubble time. The depicted case corresponds to a case with a metallicity of Z = Z_⊙/ 20, which provides a black hole merger after 2.6 Gyr.

Fig. 1 Evolutionary tracks of both stars in a MOB in the HR diagram. The initial masses are 79 M⊙ and 64 M⊙ and the initial orbital period is 1.1 d. Both stars evolve towards nearly equal masses, such that their evolutionary tracks after the overcontact phase become almost identical.

Fig. 2 Illustration of the binary stellar evolution leading to a BH+BH merger with a high chirp mass. The initial metallicity is Z⊙/ 50, the masses of the stars in solar masses are indicated with red numbers, and the orbital periods in days are given as black numbers. A phase of contact near the ZAMS causes mass exchange. Acronyms used in the figure: ZAMS: zero-age main sequence; TAMS: termination of hydrogen burning; He-star: helium star; SN: supernova; GRB: gamma-ray burst; BH: black hole.

3.2. Example grid

An example of a grid of binary systems is shown in Fig. 3 for Z = Z_⊙/ 50 and q_i = 1. Each rectangle in the plot corresponds to one detailed binary evolution model.

As Fig. 3 shows, progenitors of massive double helium stars require initial primary masses above about30 M_⊙, and the range of periods for which they are formed widens with increasing primary mass. This broadening is the consequence of the larger convective cores and stronger winds for the more massive stars; this has a similar effect as the rotational mixing in exposing helium-rich material at the surface (Köhler et al. 2015; Szécsi et al. 2015).

In particular for binaries whose final masses are low enough to avoid the pair-instability regime (i.e. roughly for M_f< 60 M_⊙; cf. Sects. 3.8 and 3.9), the parameter space for progenitors is completely dominated by overcontact systems, many of which come from systems whose periods are so low that they overflow their Roche radii at the ZAMS. Only for initial masses well above 100 M_⊙ do we find systems that avoid contact throughout their evolution. But even in this mass regime, most systems undergo at least one contact phase.

Fig. 3 Example of a grid of binary systems (initial orbital period versus initial primary mass) with Z = Z⊙/ 50 and qi = 1. Models that reached a point at which the difference between the surface and central helium abundance in one of the stars exceeds 0.2 are considered not to be evolving chemically homogeneously and their calculation is stopped (pink). The region in which the initial orbital period is small enough as to have L2 overflow at the ZAMS is marked in black, while those systems that reach L2 overflow during the main sequence are marked in green. Systems marked in blue successfully form double helium stars. Single hatching marks systems that experience contact during the main sequence, while doubly hatched systems are in an overcontact phase at the ZAMS.

We stopped all but three (see Sect. 3.8) of our binary-evolution models at a time when the stars ended core helium burning since their fate is settled at that time, and the binary orbit will essentially not change any more until the first stellar collapse occurs (third stage in Fig. 2).

3.3. Final binary configurations

Figure 4 summarizes the distribution of the final total system masses as a function of their final orbital period for those models in our grid that succeeded in producing close pairs of helium stars. Since the initial binary periods have to be very short to enforce the rapid rotation required for homogeneous evolution, the final properties lie in a narrow strip for each metallicity, but these are distinctly different for different metallicities. For the highest considered masses, this is mainly due to the metallicity dependence of the stellar wind mass loss, which has the effect of widening the systems and reducing the mass of the stars, thus producing systems with longer orbital periods and lower masses at higher metallicity.

Figure 4 also indicates the merger times for these systems, assuming that the masses and periods do not change in the black hole formation process (cf. Sect. 3.9). All models with Z_⊙/ 4 and all but the lowest mass models with Z_⊙/ 10 produce binaries that are too wide to lead to black hole mergers within a Hubble time. The more metal-poor models, on the other hand, produce very tight He-star binaries below as well as above the mass regime where pair-instability supernovae (PISNe) are expected to lead to the complete disruption of the stars and not to the formation of black holes (Heger & Woosley 2002; Chatzopoulos & Wheeler 2012).

Fig. 4 Total masses and orbital periods at core helium depletion for systems with qi = 1 at four different metallicities. Dashed lines show constant merger times assuming direct collapse into a black hole, and the shaded region indicates the mass range at which PISNe would occur, resulting in the total disruption of the stars instead of black hole formation. The coloured bands represent the relative number of objects formed for each metallicity.

The trend of shorter merger times for lower metallicities is expected to continue towards the lowest metallicities found in the Universe. As stellar wind mass loss becomes increasingly negligible, the initial stellar radii determine the shortest possible orbital periods. As an example, stars of 60 M_⊙ have ZAMS radii of 12 R_⊙, 10.5 R_⊙, 10 R_⊙, and 3.5 R_⊙, at Z = Z_⊙, Z_⊙/ 10, Z_⊙/ 50, and Z = 0, respectively. This implies that the merger times for the lowest metallicities, in particular for Population III stars, become extremely short. While the expected number of such objects is small, this opens the exciting possibility of eventually observing primordial black hole mergers at high redshift.

3.4. Mass distribution and mass ratios

Figure 5 shows the predicted intrinsic chirp mass distribution for BH+BH mergers for our different metallicity grids, again assuming no mass loss in the BH formation process. The most prominent feature is the prediction of a clear gap in this distribution, which occurs because systems that would otherwise populate this gap do not appear since the stars explode as pair-instability supernovae without leaving a stellar remnant. The BH progenitors in the systems above the gap also become pair unstable, but the explosive burning cannot reverse the collapse, which leads straight to the formation of a black hole (Heger & Woosley 2002; Langer 2012).

There is a strong general trend towards higher chirp masses with decreasing metallicity. At the lowest metallicity (Z = Z_⊙/ 50) we also produce BHs above the PISN gap. While obviously there are fewer of them than BH systems below the gap, they may still be significant because the amplitude of the gravitational-wave signal is a strong function of the chirp mass (cf. Sect. 4).

As indicated in Fig. 5, the vast majority of merging systems have passed through a contact phase. Since both stars are relatively unevolved when they undergo contact, these contact phases result in mass transfer back and forth until a mass-ratio q ≃ 1 is achieved. This is depicted in Fig. 6, where final mass ratios are shown for systems with q_i = 0.9and0.8 and a metallicity of Z = Z_⊙/ 50. For each mass ratio, two distinct branches are visible, corresponding to systems that undergo contact and evolve to q ≃ 1, and systems that avoid contact altogether. Owing to the strong dependence of mass-loss rates with mass, even systems that avoid contact altogether evolve towards q = 1 at high masses.

Mandel & de Mink (2016) modelled this channel without contact systems and found that many binaries form double BHs from progenitors below the PISN gap, with final mass ratios in the range of 0.6 to 1, reflecting just a small shift from the initial mass ratio distribution as a result of mass loss. However, Mandel & de Mink (2016) did not perform detailed stellar evolution calculations. They checked whether their binary components underfilled their Roche radii at the ZAMS and then assumed that this will remain so in the course of the quasi-homogeneous evolution of both stars. When considered in detail, however, in particular the more massive and more metal-rich stars undergo some expansion during core hydrogen burning, even on the quasi-homogeneous path (Brott et al. 2011; Köhler et al. 2015; Szécsi et al. 2015). This is most likely due to the increase of their luminosity-to-mass ratio and the related approach to the Eddington limit (Sanyal et al. 2015). As a result, the vast majority of the binaries considered by Mandel & de Mink (2016) enter contact when computed in detail. Therefore, our final mass ratio distribution is much more strongly biased towards q = 1.

Fig. 5 Stacked distribution of chirp masses of BH+BH systems formed at different metallicities, so that they merge in less than 13.8 Gyr. The contribution from each metallicity is scaled assuming a flat distribution in Z. At very short periods, systems are in contact at the ZAMS.

Fig. 6 Mass-ratios of BH+BH systems resulting from our modelled systems for qi = 0.9 and qi = 0.8 and a metallicity of Z = Z⊙/ 50 under the assumption that no mass is lost during collapse. The shaded region indicates the limits for the occurrence of PISNe.

3.5. Merger delay times

As Fig. 4 indicates, the merger delay time, which is the time between the formation of the BH+BH binary and the eventual merger, is a strong function of metallicity, where the merger delay times (at a given BH mass) are systematically shorter for lower metallicity. Figure 7 shows the distribution of the merger delay times for the different metallicities in our grids (assuming a uniform metallicity distribution). While the typical delay time is several Gyr, which helps detecting these events at lower redshift (see Sect. 4), the delay time can be as short as 0.4 Gyr for BH+BH mergers below the PISN gap at the lowest metallicity.

The decrease in delay times with lower metallicity is not found in the models of Mandel & de Mink (2016), who concluded that no high-redshift mergers are expected. The reason is that they effectively only considered one metallicity, namely the threshold metallicity for chemically homogeneous evolution by Yoon et al. (2006) of Z = 0.004. However, the components of lower metallicity binaries are more compact, allowing tighter binaries at zero age, and they have weaker winds, which produces tighter double black hole binaries. Therefore, while Mandel & de Mink (2016) predicted delay times to be longer than 3.5 Gyr, we found up to ten times lower values at our lowest metallicity (Fig. 7). Since the shortest delay time depends on the metallicity-dependent stellar radii and stellar radii of massive metal-free stars are smaller than half compared to those at Z_⊙/ 50 (Yoon et al. 2012; Szécsi et al. 2015), even merger times shorter by orders of magnitude can be expected. Therefore, even though much rarer, we argue that massive BH mergers could occur up to the redshift of Population III stars. If such mergers were detected, it would allow us to probe the evolution of massive stars in the very early Universe.

We also note that if the black holes receive kicks at birth, even higher metallicity systems may merge very rapidly if the kick reduces the pericenter distance (see Appendix A).

Fig. 7 Stacked distribution of merger delay times for different metallicities (as indicated). The meaning of the hatching is the same as in Fig. 5.

Fig. 8 Angular-momentum profiles at core helium depletion for the primary stars of binaries from our grid that result in double helium star binaries. We show stars of three different initial masses in binaries with similar initial orbital periods at metallicities of Z = Z⊙/ 50,Z⊙/ 20, and Z⊙/ 10. The curves for the specific angular momentum of the last stable orbit for a non-rotating (Schwarzschild) and critically rotating (Kerr) black hole are also included.

3.6. Observational counterparts

Our choice of including models of up to 500M_⊙ in our grids is supported by the evidence for stars of such high masses in the LMC (Crowther et al. 2010).

Various evolutionary stages of the MOB scenario are observationally confirmed by massive binary systems in nearby galaxies. As mentioned in Sect. 2, the MOB VFTS352 (Almeida et al. 2015) supports the idea of homogeneous evolution of overcontact binaries, even though it is not expected to lead to a black hole merger because of the rather high metallicity of the LMC (cf. Fig. 4). It corresponds to the first stage of our cartoon in Fig. 2.

The SMC binary HD 5980 corresponds well to the second stage of Fig. 2. It consists of two stars with masses of about 60 M_⊙ that are both very hydrogen poor, that orbit each other in about19 d. Koenigsberger et al. (2014) concluded that this system most likely emerged from homogeneous evolution. This system is well recovered in our grid at Z_⊙/ 10.

Finally, IC10 X-1 and NGC 300 X-1 are binaries that may correspond well to stage 3 of Fig. 2. Both have a short orbital period (P_orb ≃ 1.5 d for both) and contain very massive black hole primaries (>23 M_⊙ and 20 M_⊙) and similar-mass, hydrogen-free companions (~35 M_⊙ and 26 M_⊙; Barnard et al. 2008; Bulik et al. 2011). Both systems have close-matching counterparts in our Z_⊙/ 20 binary-evolution grids, with life times of up to several 10⁴ yr.

3.7. Spins

To test the possibility of producing LGRBs according to the “collapsar” scenario (Woosley 1993) from our MOB models, we compare in Fig. 8 angular-momentum profiles at the point of core helium depletion for a few systems that fall below the PISN gap. A significant amount of mass ejected during an LGRB event could modify the final orbital periods of double BHs, although we find that this probably does not play a determining role in our rate estimates (as discussed further in Sect. 4).

As Fig. 8 shows, models at a metallicity of Z_⊙/ 10 experience significant braking due to winds, and thus they are unlikely to produce LGRBs. In contrast, several systems at Z_⊙/ 50 that result in helium stars below the PISN gap have specific angular-momentum profiles above the values for the last stable orbit, assuming all mass collapses into a critically rotating black hole. The results at Z_⊙/ 20 are more ambiguous, and it is not clear whether the stars would produce an LGRB or not. For systems forming black holes above the PISN gap, wind braking is strong enough even for low metallicity to avoid the formation of LGRBs. This is confirmed when considering the Kerr parameters of our models in Fig. 9 in the different mass and metallicity regimes.

Fig. 9 Kerr parameter as function of the final system mass for our models at Z = Z⊙/ 50,Z⊙/ 20, and Z⊙/ 10, assuming a complete collapse of our helium stars to black holes. Binaries indicated by symbols with a red frame have merger times that exceed the Hubble time.

3.8. Models up to core collapse

To depict the effect of the PISN gap, we took three models with masses 200 M_⊙, 90 M_⊙ and 35 M_⊙ and metallicity Z_⊙/ 50 after helium depletion and evolved these through the late evolutionary phases. Figure 10 shows the evolution of central density and temperature of each star, together with the region in which pair production results in an adiabatic index of Γ < 4/3.

The least massive of the three stars avoids the pair-unstable region altogether and experiences core collapse after silicon depletion. In the 90 M_⊙ model, the core collapses during oxygen burning, resulting in explosive burning that injects enough energy to halt the collapse and drive an explosion. At the highest mass, the oxygen core also collapses, but explosive burning is not sufficient to stop it, and in the end burning proceeds very fast up to silicon depletion, resulting in an iron core with an infall velocity >1000 kms^-1.

Fig. 10 Evolution in the Tc−ρc-diagram for the three stellar models at Z = Z⊙/ 50 (with the masses at helium depletion as indicated) calculated to the final evolutionary stage. The shaded region shows the region that is unstable to pair creation. Both the 35 M⊙ and the 200 M⊙ stars collapse to form black holes, while the 90 M⊙ is disrupted in a PISN.

3.9. Explosive mass loss and momentum kicks

In all models below the pair-instability regime we expect the formation of black holes. If the whole star collapses without ejecting any mass or energy, the masses and periods in Fig. 4 would also represent the masses of the final black holes and the post-collapse orbital periods. On the other hand, as our helium stars tend to be rapidly rotating, some of them may experience a collapsar phase (Woosley 1993), producing LGRBs, in which part of the collapsing star is ejected, and the binary orbit may receive a supernova kick. The effect of the mass loss would be to reduce the final black hole masses (and to reduce the strength of any eventual gravitational-wave signal) and widen the system (and increase the merger time), while the effect of a kick can be to either increase or decrease the orbital period and the merger time (see Appendix A for a more detailed discussion). While the details of the collapse phase are still very uncertain, which may have an effect on the BH+BH detection rates, our main conclusions are not dependent on them.

The final angular-momentum profiles of our models (see Sect. 3.7) suggest that only the lowest mass models (M_final ≲ 40 M_⊙) at the two lowest metallicities (Z = Z_⊙/ 20, Z_⊙/ 50) may retain enough angular momentum in the core to be good LGRB candidates. Nevertheless, because of the large amount of available angular momentum, we expect many of the BHs formed in this scenario to be rapidly rotating, with the spin parameter roughly scaling inversely with the final orbital period shown in Fig. 4 (i.e. the fastest spins are expected for the lowest mass BHs at the lowest metallicity). Finally, we note that a regime of pulsational PISNe (Chatzopoulos & Wheeler 2012) lies below the disruptive PISN regime, where substantial mass loss is expected, but a BH is ultimately formed (Woosley et al. 2007).

4. Merger rates

For the conventional scenario in which close double compact binaries are produced through common-envelope evolution (see Appendix B), except for a few cases (Voss & Tauris 2003; Belczynski et al. 2010; Dominik et al. 2015), the far majority of published population synthesis studies predict a much higher NS+NS merger rate per Milky Way equivalent galaxy (MWEG) than they do for the rate of BH+BH mergers. Based on a detailed comparison study of published models (Abadie et al. 2010), the NS+NS merger rate was estimated to be 100 MWEG^-1 Myr^-1, which is about 100 times higher than the rate predicted for BH+BH binaries. However, because there are more massive compact objects in BH+BH binaries than in NS+NS binaries, their emitted gravitational-wave amplitudes are significantly larger, so that the LIGO detection rates of both are approximately equal. The so-called realistic rates quoted by Abadie et al. (2010) are 40 and 20 yr^-1 for NS and mostly low-mass BH mergers, respectively, but the uncertainty in these numbers is larger than three orders of magnitude.

4.1. Chirp masses and MOB evolution

Two important assumptions need to be kept in mind with regard to these quoted rates. First, the BH+BH binaries are assumed to be composed of 10 M_⊙ BHs (even 5 M_⊙ BHs in all LIGO result papers published before 2010), corresponding to an intrinsic chirp mass, ℳ₀ of ~8.7 M_⊙ (for equal mass binaries, ℳ₀ = (1/4)^3/5M ≃ 0.435 M, where the total mass, M, is twice the BH mass, M_BH). For low metallicities, our MOB scenario predicts the formation of BHs with masses of 25−60 M_⊙ and 130−230 M_⊙ (i.e. below and above the PISN mass range, respectively), resulting in very high intrinsic chirp masses of about 20−50 M_⊙ and 115−200 M_⊙, respectively (cf. Fig. 5). Mergers like this can be seen throughout a significant part of the Universe, since the distance luminosity, ${\mathit{d}}_{\mathrm{L}}\mathrm{\propto }{\mathrm{ℳ}}_{\mathrm{0}}^{\mathrm{5}\mathit{/}\mathrm{6}}$. We note that the detected chirp masses will be redshifted to ℳ₀ (1 + z), where z is the BH+BH system redshift with respect to the detector on Earth (Finn 1996).

Secondly, all previously published rates were based on common-envelope-evolution (cf. Appendix B), which creates uncertainties in the rates by more than two orders of magnitude as a result of our poor understanding of the common-envelope physics (Dominik et al. 2012). The new BH+BH formation scenario (MOB evolution) presented in this work does not involve any common-envelope phase. Instead, it is based on much less uncertain physics (as discussed in the previous section). Equally important, it leads to the formation of much more massive BHs than in previous studies.

Assuming, as a first approximation, that the detection rate, ℛ, scales with ${\mathit{d}}_{\mathrm{L}}^{\mathrm{3}}\mathrm{\propto }{\mathrm{ℳ}}_{\mathrm{0}}^{\mathrm{5}\mathit{/}\mathrm{2}}$, we conclude that the expected LIGO detection rates for these massive BH+BH binaries could easily dominate the overall rates; they are therefore excellent candidates for the first LIGO source detection (see the more detailed discussion below). It should be noted that some previous studies (Belczynski et al. 2010; Dominik et al. 2015; Rodriguez et al. 2015) have alluded to a dominance of relatively massive BH+BH mergers in a low-metallicity environment (or, in particular, through dynamical channels in dense clusters), although without a specific detailed model for the binary case.

The expected LIGO detection rate of BH+BH binary mergers was estimated in the following manner, ℛ = r_MW × N_gal, where r_MW is the expected merger rate in an MWEG, and ${\mathit{N}}_{\mathrm{gal}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\hspace{0.17em}\mathit{\pi }\hspace{0.17em}{\left(\frac{{\mathit{d}}_{\mathrm{horizon}}}{\mathrm{Mpc}}\right)}^{\mathrm{3}}\hspace{0.17em}\mathrm{(}\mathrm{2.26}{\mathrm{)}}^{-3}\hspace{0.17em}\mathrm{(}\mathrm{0.0116}\mathrm{)}$(3)is the number of MWEGs out to a horizon distance, d_horizon (Abadie et al. 2010). Here the factor 1/2.26 is included to average over all binary orientations and sky locations, i.e. d_horizon = 2.26 d_avg (Finn & Chernoff 1993), and 1.16 × 10^-2 Mpc^-3 is the extrapolated space density of MWEGs (Kopparapu et al. 2008). For relatively low-mass BH+BH mergers, assuming M_BH = 10 M_⊙ and a corresponding average design distance luminosity of d_L ≃ 1000 Mpc for advanced LIGO (aLIGO), the estimated values are (Abadie et al. 2010) r_MW = 0.4 Myr^-1 MWEG^-1 and ℛ = 20 yr^-1 . For a massive BH+BH merger with M_BH = 60 M_⊙ (or 130 M_⊙), we obtain d_L ≃ 4.5 Gpc (or d_L ≃ 8.5 Gpc), and thus d_horizon ≃ 10 Gpc (or d_horizon ≃ 19 Gpc). The expected redshift is z = 1.4 (or z = 2.3) for standard cosmological parameters (H₀ = 69.6 km s^-1 Mpc^-1, Ω_M = 0.286, Ω_vac = 0.714).

4.2. Model assumptions for estimating aLIGO detection rates

To calculate the aLIGO detection rate of these massive BH+BH mergers, we first need to calculate the intrinsic merger rate in an MWEG. For a given metallicity, it is straightforward to calculate the rate for various events from our binary evolution grids. We used simple standard assumptions about the initial binary parameters: we assumed that (1) the orbital period distribution is flat in log P and covers the range (0.5 d−1 yr); (2) the primary mass distribution is described by a Salpeter power law (dN/d$\log \mathit{M}\mathrm{\propto }{\mathit{M}}_{\mathrm{1}}^{-1.35}$); (3) the mass-ratio distribution is flat; (4) stars more massive than 8 M_⊙ produce a core-collapse SN; and (5) there is one binary system for every three core-collapse SNe (i.e. two out of three massive stars are formed in close binaries).

Guided by the results from our grids, we required that the initial mass ratio has to be higher than 0.8 to ensure chemically homogeneous evolution for both stars. With these assumptions we calculated the fraction of systems that produce BH+BH mergers relative to the core-collapse SN rate. These are shown in the first two rows in Table 1 for various metallicities for BH+BH mergers below and above the PISN gap. These numbers imply that for Z<Z_⊙/ 10, there is typically one BH+BH merger event for about1000 core-collapse SNe. To relate these numbers to a rate for an MWEG, we need to multiply these fractions by the core-collapse rate for an MWEG. Adopting a typical rate of 0.01 yr^-1, this implies, for instance, that in an MWEG with Z = Z_⊙/ 50, the BH+BH merger rates are 6.7 Myr^-1 and 2.7 Myr^-1 below and above the PISN gap, respectively. We note that our model only predicts BH+BH mergers above the PISN gap at the lowest metallicity.

4.3. Accounting for the star-forming history and the galactic metallicity distribution throughout the Universe

The actual aLIGO detection rate depends on two main factors: (1) the detection volume within which a particular event can be detected; and (2) the cosmological distribution of the sources, which depends on the star formation history and the chemical evolution of the Universe. The detection volume can be estimated using Eq. (3), where d_horizon depends on the aLIGO sensitivity and the BH masses. However, Eq. (3) is only a good approximation for d_horizon ≲ 1 Gpc because it ignores cosmological expansion. To take this into account (in a very approximate way), we introduced a simple cut-off of N_gal = 10¹⁰ in Eq. (3). For reference, this implies that our simple model assumes that there are about three core-collapse SNe per second in the Universe. Scaling the horizon distance of the BH+BH masses produced in our grids to the design sensitivity of aLIGO, we started by calculating the aLIGO detection rate assuming that all galaxies have a particular metallicity. These rates are shown in the last two rows of Table 1, in Cols. 2–5.

One of the main problems is that our rates are a strong function of metallicity and therefore depend on the evolution of metallicity with time and the spread of the metallicity distribution at a given redshift. For example, the mean metallicity of galaxies will only be lower than Z_⊙/ 50 at redshifts higher than 5−7 (A. Fruchter, priv. comm.). On the other hand, even in the local Universe some galaxies (mostly dwarf galaxies) have extremely low metallicities. A proper calculation of this is beyond the scope of this paper. However, we can derive lower and upper limits using a local and a global approximation of the metallicity distribution.

For the former, we followed Langer & Norman (2006), who computed the formation rate of stars with a metallicity below a threshold metallicity Z^∗ in the local Universe. For their fiducial parameters, they found fractions compared to the local total star formation rate of 0.61, 0.01, 0.0025, and 0.0004 for the four metallicities used in our binary evolution models (Z_⊙/ 4, Z_⊙/ 10, Z_⊙/ 20, and Z_⊙/ 50, respectively). If the merger delay times were negligible (which they are not; see below), and the detection volume were restricted to the low-redshift Universe – which, from the point of view of the statistics mentioned above remains roughly true for redshifts up to 2 ~ 3 – the above factors would need to be applied to obtain local detection rates for the various metal-poor sources.

A global upper limit to the rates can be obtained by considering the metallicity distribution of all massive stars that have ever formed in the Universe. Using the global metallicity distribution provided by C. Kobayashi (based on the simulations presented in Taylor & Kobayashi 2015), we estimate the fraction of massive stars that have formed with metallicities of Z_⊙/ 10, Z_⊙/ 20, and Z_⊙/ 50 to be 0.086, 0.052, and 0.068, respectively (using appropriate linear binning).

4.4. Resulting predictions for aLIGO detection rates

Using these metallicity weightings, we can estimate ranges for the aLIGO detection rates (at the design sensitivity) of 19−550 yr^-1 for BH+BH mergers below the PISN gap and of 2.1−370 yr^-1 above the PISN gap, compare the last column in Table 1. Here, the first quoted number of the ranges corresponds to the local approximation, the second to the global approximation. Even for the ongoing first science run (O1) of aLIGO the prospects for detection are promising. Given that the sensitivity of aLIGO is currently about one-third of the design sensitivity, the expected detection rate is ~4% of our calculated values (see last column of Table 1).

The lower mass BH+BH mergers are more likely to sample the low-redshift Universe, therefore the lower limit may be more applicable for the mergers below the PISN gap. On the other hand, because the most massive BH+BH mergers can be detected throughout most of the visible Universe, the upper limit may be more appropriate for the mergers above the PISN gap. We note that for these, the redshift factor 1/(1 + z) has not yet been taken into account, as we did not compute the redshift distribution of events. Even our lower limits suggest, however, that aLIGO should detect BH+BH mergers from the MOB scenario. The most massive mergers, which probe a large portion of the whole Universe, might well be the dominant source for aLIGO detections (Flanagan & Hughes 1998; Abadie et al. 2010).

Another factor that helps detecting BH+BH mergers from low-metallicity populations is not taken into account in the above estimates: because of the possibly long merger delay times (see Fig. 7), even systems that were formed in the early Universe may merge at low redshift (cf. Mandel & de Mink 2016).

A remaining caveat of concern for our estimated detection rates is related to the relatively low gravitational-wave frequencies of the more massive BH+BH binaries. The emitted frequencies during in-spiral are expected to peak approximately at the innermost stable circular orbit (ISCO) before the plunge-in phase and the actual merging. However, even determining the ISCO for a merging binary system is non-trivial and depends, for example, on the spins of the BHs (Balmelli & Damour 2015). This requires numerical or sophisticated (semi-)analytical calculations within general relativity and cannot simply be estimated using a test particle in a Kerr field. For the BH+BH mergers above the PISN gap, the emitted frequencies are most likely ≤100 Hz, and with redshift corrections, the frequencies to be detected are easily lower by a factor of two or more. A frequency this low is close to the (seismic noise) edge of the detection window of aLIGO. However, the waveform amplitudes of the more massive BH+BH binaries are larger (for a given distance) and are also enlarged further by a factor of (1 + z). Finally, it may be possible that higher frequency signals from the ring down of the single rapidly spinning BH produced might be detectable, despite their expected smaller wave amplitudes.

An important question to address is whether the first generation of LIGO should have detected such massive BH+BH merger events. Given that the sensitivity of the first generation of LIGO was about ten times lower, the number of detections should have been 1000 times lower. Even for our upper limits, it is therefore not surprising that there have been no detections during the previous science runs of the first-generation LIGO detectors (Abadie et al. 2012).

5. Concluding remarks

We emphasize that, unlike other channels, the MOB channel for the formation of merging BH binaries is quite robust, relying on reasonably well understood stellar evolution physics. The main uncertainty is the treatment of mixing in rapidly rotating stars, but even here we can derive some confidence from the fact that our models are able to reproduce observed local counterparts of various stages in the MOB scenario (see Fig. 2), such as HD 5980, IC10 X-1, and NGC 300 X-1. The MOB channel predicts the formation of very massive compact BH+BH binaries with a BH mass ratio close to 1 and a bimodal BH-mass distribution from BHs formed below and above the PISN regime.

The detection of GWs from BH+BH mergers with LIGO (and potentially other current or future GW detectors) will not only start a revolution in observational astronomy and test general relativity in its highly dynamic strong-field regime, but will also have a great effect on our understanding of very massive stars throughout the Universe, including their fate as gamma-ray bursts or pair-instability supernovae.

Acknowledgments

P.M. and N.L. are grateful to Bill Paxton for his continuous help in extending the MESA code to contain all the physics required for this project over the last years. We are thankful to Selma de Mink and Norbert Wex for helpful discussions, and to Ed van den Heuvel for useful comments on an earlier version of this paper. We thank Chiaki Kobayashi for providing us with the global metallicity distribution of massive stars and Andy Fruchter for plots of the mean metallicity evolution with redshift. N.L.’s Alexander von Humboldt Professorship and PP’s Humboldt Research Award provided essential support for this research. The Geryon2 cluster housed at the Centro de Astro-Ingenieria UC was used for the calculations performed in this paper. The BASAL PFB-06 CATA, Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130008 provided funding for several improvements to the Geryon2 cluster.

References

Appendix A: Dynamical implications of black hole kicks

After a tight binary system with two BHs is formed, the continuous radiation of gravitational waves will cause loss of orbital energy and a resulting shrinkage of the orbit (Einstein 1918). The timescale until the two binary components finally merge depends on the BH masses, the orbital period, and the eccentricity of the system (Peters 1964). The initial parameters of BH+BH binaries depend not only on the stellar evolution of the progenitor stars, but also on the physics of the BH formation process itself.

The formation of BHs may be accompanied by a momentum kick, similar in nature to those imparted onto newborn NSs (Janka 2012), in particular if the core of the progenitor does not collapse directly to a BH but in a two-step process (Brandt et al. 1995). However, while the magnitude of such kicks for NSs is fairly well constrained from pulsar observations (Hobbs et al. 2005), the magnitude of kicks imparted onto a BH during its formation is rather uncertain. So far, the measured masses of observed stellar-mass BHs are all relatively low: 4−16 M_⊙ for all the Galactic sources (McClintock et al. 2014), about16 M_⊙ for M33 X-7 (Orosz et al. 2007) and 24−33 M_⊙ for IC 10 X1 (Prestwich et al. 2007).

Fig. A.1 Simulations of the dynamical effects of 10 000 asymmetric SNe on BH binary systems. The initial pre-SN system contains two 50 M⊙ stars in a binary with an orbital period of 2.0 days. During the formation of a BH, it is assumed that 10 M⊙ is lost instantaneously and a kick velocity is imparted with a magnitude between 0−300 km s-1, drawn from a flat distribution and an isotropic (random) distribution of kick directions. The upper left panel shows the resulting post-SN systems in the orbital period-eccentricity plane. The colours indicate the merger time of the post-SN system as a result of gravitational-wave radiation (see the distribution in the upper right panel). The lower panels show the distributions in orbital period and eccentricity of the post-SN binaries. The case of a purely symmetric SN (no kick) is indicated with a black dot in all panels.

Their inferred BH kicks from observations and theoretical arguments span from basically no kicks (Nelemans et al. 1999) to BH kicks of several 100 km s^-1 (Janka 2013) based on hydrodynamical kicks associated with asymmetric mass ejection and subsequent BH acceleration.

For massive stellar mass BHs, the situation is different since their progenitor stars more likely collapse directly to form a BH without a SN explosion, leading to BHs with very low kick velocities. Thus, a bimodality of the BH kick velocity distribution seems possible (Janka 2013). Another, and possibly related, problem is the amount of mass loss during BH formation (ejected baryonic mass and rest mass energy carried away by neutrinos), which also affects the orbital period and the eccentricity of the post-collapse system.

To quantify the combined effects of mass loss and kicks, we assume in the following example that a 50 M_⊙ Wolf-Rayet (WR) star collapses to form a BH with a gravitational mass of 40 M_⊙. In other words, we assume that 10 M_⊙ of mass is lost by a combination of baryonic mass loss and/or losses through neutrinos from outside the event horizon. A large amount of baryonic mass loss may be expected if the progenitor core rotates rapidly and the collapse is associated with an LGRB and an LGRB supernova.

We first consider the no-kick case for a circular pre-collapse system with two 50 M_⊙ stars and an orbital period of 2.0 days.

If there is no mass loss and no kick imparted to the newborn BH, the orbital parameters of the binary system remain unchanged and the system would merge in 550 Myr. The result of an instantaneous mass loss of 10 M_⊙ during the BH formation produces a post-collapse system with an orbital period of 2.52 days and an eccentricity of 0.11 (and stellar masses of 40 M_⊙ and 50 M_⊙), which will merge in 1180 Myr. If BH formation is accompanied by an additional momentum kick, however, the outcome can change significantly.

In Fig. A.1 we plot the dynamical consequences for a surviving binary system, using the same initial parameters as before, in which a BH is produced with an asymmetric kick velocity and instantaneous mass loss, following the recipe of Hills (1983). In a BH+BH binary, two kicks may be imparted, but here we restrict our example to just one kick to better illustrate the principal dynamical effects – the analysis can easily be generalized to two kicks without changing the main conclusions. We performed 10 000 trials, assuming a flat distribution of the BH kick magnitudes between 0−300 km s^-1 and an isotropic (random) distribution of directions. If a randomly orientated kick of a fixed magnitude of 300 km s^-1 (600 km s^-1) were applied in all cases, it would result in the disruption of about 7.3 per cent (43 per cent) of the cases. Using a flat distribution between 0−300 km s^-1, only 0.4 per cent of all systems are disrupted. Thus even for a wide range of assumed BH kick values, the survival rate of BH+BH binaries only changes by less than a factor of two.

More important for our investigation here is the merger timescale that is due to gravitational-wave radiation. Figure A.1 shows that the effect of a kick can either widen or shorten the post-collapse orbit. However, the merger timescale for a binary with given component masses is a function of both orbital period and eccentricity. Given that systems which widen more during BH formation (as a consequence of instantaneous effective mass loss and a kick) are also the systems attaining the higher eccentricities, the net effect of a kick on the resulting merger time is surprising small. For BH formation without a kick, the merger timescale of the resulting binary is 1180 Myr. Using a flat kick distribution between 0−300 km s^-1 results in roughly half (47 per cent) of the surviving systems to merge on a shorter timescale, and the other half to merge on a longer timescale (and only 2.5 per cent exceeding a Hubble time).

Applying a strong kick of 600 km s^-1 actually causes a larger fraction of surviving systems to merge (67 per cent) on a shorter timescale (<1180 Myr) than in the symmetric case without a kick. Therefore, we can safely conclude that although BH kicks may, in general, widen a number of systems, the resulting merger timescale distribution will not change one of the main findings of this paper, namely that the LIGO detection rate is likely to be dominated by quite massive BH+BH mergers.

Appendix B: Comparison to the standard BH+BH formation scenario

The standard formation scenario of BH+BH binaries involves a number of highly uncertain aspects of binary interactions (Fig. B.1, Tauris & van den Heuvel 2006). The main uncertainties include, in particular, the treatment of common-envelope evolution (Ivanova et al. 2013) and the efficiency of accretion and spin-up from mass transfer. These lead to uncertainties in the expected merger rates of several orders of magnitude (Abadie et al. 2010). In contrast, the MOB scenario presented here mostly relies on reasonably well understood physics of the evolution of massive stars, although there are still significant uncertainties, for instance, in the treatment of stellar winds (Langer 2012), rotational mixing (Maeder & Meynet 2012), and the BH formation itself (Heger et al. 2003; Ugliano et al. 2012; Pejcha & Thompson 2015).

Fig. B.1 Illustration of the binary stellar evolution leading to a BH+BH binary according to the standard scenario. Acronyms used in the figure. ZAMS: zero-age main sequence; RLO: Roche-lobe overflow (mass transfer); WR-star: Wolf-Rayet star; SN: supernova; BH: black hole; HMXB: high-mass X-ray binary; CE: common envelope.

To produce a tight BH binary that will merge within a Hubble time, we could consider a very close massive binary system with an initial orbital period of a few days. The problem with such a model is that these systems would mostly be expected to merge early during their evolution when the more massive star evolves off the main sequence, expands, and starts to transfer mass to its companion star. Although other models have been proposed that evolve without a common-envelope phase, for example to explain the formation of IC 10 X-1 (de Mink et al. 2009) and M33 X-7 (Valsecchi et al. 2010), they often require some degree of fine-tuning to work and did not follow the evolution to the end to produce a binary with two BHs.

To avoid this problem, it has been common practice to model the formation of BH+BH systems starting from relatively wide systems and let the systems evolve through a common-envelope phase following the high-mass X-ray binary (HMXB) phase after the formation of the first BH (see Fig. B.1). There are currently no self-consistent hydrodynamical simulations for modelling the spiral-in of BHs inside a massive envelope; in particular, it is unclear whether the BH will experience hypercritical accretion, and under which conditions the systems will merge completely; all of this leads to large uncertainties in the number of post-common-envelope systems and their separations (Ivanova et al. 2013).

Another problem that is often ignored is that the fate of a massive star in a binary depends on when it loses its hydrogen-rich envelope. As first pointed out by Brown (Brown et al. 1999) and confirmed in later calculations (Brown et al. 2001; see also Petermann et al., in prep.), if a massive star loses its hydrogen-rich envelope before or early during helium core burning, it ends its evolution with a much smaller iron core and is more likely to produce an NS than a BH. This means that the formation of a BH in a close binary may require that its progenitor loses its envelope only after helium core burning. To produce not just one, but two BHs, this may require extreme fine-tuning and may even be impossible in the standard scenario in Fig. B.1. The problem may be avoided in the so-called double-core scenario (Brown 1995; Dewi et al. 2006), where the BH progenitors both evolve beyond helium core burning before experiencing a common-envelope phase in which the cores of both stars spiral in, producing a close binary of two helium stars that subsequently collapse to form BHs. However, this scenario requires significant fine-tuning since the initial masses of the two stars have to be very close to each other (q_i> 0.96), and it may be impossible to produce quite massive BHs because very massive stars tend to avoid mass transfer after helium core burning (although this may be possible at sufficiently low metallicity (Z ≲ 0.1 Z_⊙)).

Given that all measured stellar-mass BHs in the Milky Way have masses in the range of 4−16 M_⊙, most population-synthesis models used for estimating LIGO detection rates were previously restricted to initial progenitor stellar masses of up to about 100 M_⊙. The discovery of very massive stars (Crowther et al. 2010; Hainich et al. 2014) in the R136 region of the Large Magellanic Cloud with masses of up to 300 M_⊙, however, suggests that BH+BH binaries may form with significantly more massive components, thus enabling LIGO to detect the merger of such massive BH+BH binaries, with chirp masses easily exceeding 30 M_⊙, out to long distances (see discussion in Sect. 4).

It has previously been argued (Belczynski et al. 2014; Rodriguez et al. 2015) that massive BH+BH binaries to be potentially detected by LIGO can only form through dynamical channels in dense stellar environments. We here demonstrated that close binaries in the Galactic disk with very massive stars undergoing chemically homogeneous evolution (and which therefore do not expand after leaving the main sequence) can form massive BH+BH binaries that merge within a Hubble time at sufficiently low metallicity. An important consequence of our scenario is that massive BHs of a given mass can be produced from stars with a lower ZAMS mass – especially at low metallicity – than in the standard BH+BH formation scenario.

So far, no stellar-mass BH+BH binaries have been discovered anywhere, and there is only one known potential progenitor system in the Milky Way, Cyg X-3. The nature of the compact component in this system is still uncertain. Belczynski et al. (2013) have argued that it contains a 2−4.5 M_⊙ BH and a 7.5−14.2 M_⊙ WR-star, which would make it a potential progenitor for a BH+BH system. However, the final destiny of this system is unclear, and it might also become a BH-NS binary or even NS+NS binary (if the first-formed compact object is a neutron star). It is also possible that the system is disrupted as a result of the explosion of the WR-star.

All Tables

All Figures

	Fig. 1 Evolutionary tracks of both stars in a MOB in the HR diagram. The initial masses are 79 M⊙ and 64 M⊙ and the initial orbital period is 1.1 d. Both stars evolve towards nearly equal masses, such that their evolutionary tracks after the overcontact phase become almost identical.
In the text

Fig. 2 Illustration of the binary stellar evolution leading to a BH+BH merger with a high chirp mass. The initial metallicity is Z⊙/ 50, the masses of the stars in solar masses are indicated with red numbers, and the orbital periods in days are given as black numbers. A phase of contact near the ZAMS causes mass exchange. Acronyms used in the figure: ZAMS: zero-age main sequence; TAMS: termination of hydrogen burning; He-star: helium star; SN: supernova; GRB: gamma-ray burst; BH: black hole.

In the text

Fig. 3 Example of a grid of binary systems (initial orbital period versus initial primary mass) with Z = Z⊙/ 50 and qi = 1. Models that reached a point at which the difference between the surface and central helium abundance in one of the stars exceeds 0.2 are considered not to be evolving chemically homogeneously and their calculation is stopped (pink). The region in which the initial orbital period is small enough as to have L2 overflow at the ZAMS is marked in black, while those systems that reach L2 overflow during the main sequence are marked in green. Systems marked in blue successfully form double helium stars. Single hatching marks systems that experience contact during the main sequence, while doubly hatched systems are in an overcontact phase at the ZAMS.

In the text

Fig. 4 Total masses and orbital periods at core helium depletion for systems with qi = 1 at four different metallicities. Dashed lines show constant merger times assuming direct collapse into a black hole, and the shaded region indicates the mass range at which PISNe would occur, resulting in the total disruption of the stars instead of black hole formation. The coloured bands represent the relative number of objects formed for each metallicity.

In the text

	Fig. 5 Stacked distribution of chirp masses of BH+BH systems formed at different metallicities, so that they merge in less than 13.8 Gyr. The contribution from each metallicity is scaled assuming a flat distribution in Z. At very short periods, systems are in contact at the ZAMS.
In the text

	Fig. 6 Mass-ratios of BH+BH systems resulting from our modelled systems for qi = 0.9 and qi = 0.8 and a metallicity of Z = Z⊙/ 50 under the assumption that no mass is lost during collapse. The shaded region indicates the limits for the occurrence of PISNe.
In the text

	Fig. 7 Stacked distribution of merger delay times for different metallicities (as indicated). The meaning of the hatching is the same as in Fig. 5.
In the text

Fig. 8 Angular-momentum profiles at core helium depletion for the primary stars of binaries from our grid that result in double helium star binaries. We show stars of three different initial masses in binaries with similar initial orbital periods at metallicities of Z = Z⊙/ 50,Z⊙/ 20, and Z⊙/ 10. The curves for the specific angular momentum of the last stable orbit for a non-rotating (Schwarzschild) and critically rotating (Kerr) black hole are also included.

In the text

	Fig. 9 Kerr parameter as function of the final system mass for our models at Z = Z⊙/ 50,Z⊙/ 20, and Z⊙/ 10, assuming a complete collapse of our helium stars to black holes. Binaries indicated by symbols with a red frame have merger times that exceed the Hubble time.
In the text

	Fig. 10 Evolution in the Tc−ρc-diagram for the three stellar models at Z = Z⊙/ 50 (with the masses at helium depletion as indicated) calculated to the final evolutionary stage. The shaded region shows the region that is unstable to pair creation. Both the 35 M⊙ and the 200 M⊙ stars collapse to form black holes, while the 90 M⊙ is disrupted in a PISN.
In the text

Fig. A.1 Simulations of the dynamical effects of 10 000 asymmetric SNe on BH binary systems. The initial pre-SN system contains two 50 M⊙ stars in a binary with an orbital period of 2.0 days. During the formation of a BH, it is assumed that 10 M⊙ is lost instantaneously and a kick velocity is imparted with a magnitude between 0−300 km s-1, drawn from a flat distribution and an isotropic (random) distribution of kick directions. The upper left panel shows the resulting post-SN systems in the orbital period-eccentricity plane. The colours indicate the merger time of the post-SN system as a result of gravitational-wave radiation (see the distribution in the upper right panel). The lower panels show the distributions in orbital period and eccentricity of the post-SN binaries. The case of a purely symmetric SN (no kick) is indicated with a black dot in all panels.

In the text

	Fig. B.1 Illustration of the binary stellar evolution leading to a BH+BH binary according to the standard scenario. Acronyms used in the figure. ZAMS: zero-age main sequence; RLO: Roche-lobe overflow (mass transfer); WR-star: Wolf-Rayet star; SN: supernova; BH: black hole; HMXB: high-mass X-ray binary; CE: common envelope.
In the text