The Most Ordinary Formation of the Most Unusual Double Black Hole Merger

The Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo Collaboration (LVC) has reported the discovery of a surprisingly heavy double black hole (BH–BH) merger with component masses ${m}_{1}={85}_{-14}^{+21}\,{M}_{\odot }$ and ${m}_{2}={66}_{-18}^{+17}\,{M}_{\odot }$ and an effective spin parameter ${\chi }_{\mathrm{eff}}={0.08}_{-0.36}^{+0.27}$ at redshift z=0.82 (GW190521; Abbott et al. 2020). The corresponding merger rate density of events similar to GW190521 was estimated to be ${0.13}_{-0.11}^{+0.30}\,{\mathrm{Gpc}}^{-3}\,{\mathrm{yr}}^{-1}$.

Stars are not expected to form BHs of such masses. In particular, the pair-instability pulsation supernovae (PPSNe; Heger & Woosley 2002; Woosley et al. 2007) are associated with severe mass loss that limits BH mass and pair-instability supernovae (PSNe; (Bond et al. 1984; Fryer et al. 2001; Chatzopoulos & Wheeler 2012) are expected to completely disrupt massive stars with no resulting BH formation. These processes were believed to create the so-called "upper mass gap" in the BH mass spectrum, i.e., the lack of stellar-origin BHs in the mass range M_BH ∼ 50–135M_⊙ (Marchant et al. 2016; Mandel & de Mink 2016; Belczynski et al. 2016a; Spera & Mapelli 2017). It appeared that results of the first and second observing runs (O1/O2) of advanced LIGO/Virgo were consistent with the existence of this mass gap  (Fishbach et al. 2020). Yet, the latest LIGO/Virgo third observing run (O3) observations revealed GW190521.

This has naturally promoted proposals in which BHs in GW190521 are not products of standard stellar evolution. These proposals include dynamical formation scenarios of repeated BH mergers in dense clusters (Fragione et al. 2020; Gayathri et al. 2020; Rizzuto et al. 2020), repeated stellar mergers in dense clusters (di Carlo et al. 2019, 2020; Renzo et al. 2020), BH captures in galactic nuclei (Gondán & Kocsis 2020), and primordial black holes  (de Luca et al. 2020). Some more exotic scenarios are also being put forward such as head-on collisions of boson stars  (Calderón Bustillo et al. 2020). Alternatively, it is claimed that the LVC analysis is not the only solution to the GW190521 waveform, and the actual BH masses may be outside the upper mass gap and are consistent with standard stellar evolution  (Fishbach & Holz 2020; Moffat 2020; Nitz & Capano 2020).

In the last few years the understanding of the upper mass gap begun to change. First, it was proposed that the first population of metal-free (Population III) stars may form BHs up to  ∼70M_⊙ without violating the pair-instability physics (Woosley 2017). This was extended to  ∼85M_⊙ by recent detailed stellar evolution (Farrell et al. 2020; Tanikawa et al. 2020) and population synthesis calculations (Kinugawa et al. 2020). Second, it was proposed that for the intermediate-metallicity stars (Population II) BHs can form with masses up to 80M_⊙ (Limongi & Chieffi 2018). Third, for high-metallicity stars (Population I) the limit was increased to 70M_⊙ (Belczynski et al. 2020a). These updates on position of lower edge of the upper mass gap were the result of detailed considerations of stellar evolution processes (e.g., rotation, mixing, convection) that allow some stars to avoid the PPSN/PSN. Finally, it was shown that for low-metallicity stars (Z < 10⁻⁵–10⁻⁴), the uncertainties in the reaction rate of carbon burning (along with uncertainties on mixing/dredge-up) can potentially shift the onset of the BH upper mass gap up to 90M_⊙ (Costa et al. 2020; Farmer et al. 2020). This reaction rate concerns one of the most uncertain reactions used in stellar evolution, and yet it plays a vital important role in astrophysics (deBoer et al. 2017; Takahashi 2018; Holt et al. 2019; Sukhbold & Adams 2020).

Here, we adopt the latest results on the lower bound of the upper mass gap to test whether it is possible to (i) form BH–BH mergers with masses as reported by LVC for GW190521, and (ii) whether it is possible to form enough of them to match the LVC reported merger rate of such events. We perform our analysis in the framework of the most ordinary BH–BH merger formation scenario: the classical isolated binary evolution of Population I/II stars.

We use the population synthesis code StarTrack (Belczynski et al. 2008). We assume standard wind losses for massive stars: O/B star winds (Vink et al. 2001) and LBV winds (specific prescriptions for these winds are listed in Section 2.2 of Belczynski et al. 2010). We treat the accretion onto compact objects during the Roche lobe overflow (RLOF) and from stellar winds using the analytic approximations presented by King et al. (2001) and by Mondal et al. (2020), and limit accretion during the common envelope (CE) phase to 5% of the Bondi rate  (MacLeod et al. 2017). We employ the delayed core-collapse supernova (SN) engine in NS/BH mass calculation (Fryer et al. 2012) that allows for populating the lower mass gap between NSs and BHs (Belczynski et al. 2012; Zevin et al. 2020). The most updated description of StarTrack is given by Belczynski et al. (2020b) and the model M30 in this study describes our standard choices of input physics. In our study we employ the fallback decreased NS/BH natal kicks with σ = 265kms⁻¹, we do not allow CE survival for Hertzsprung gap donors (submodels B in our past calculations), and we assume a 100% binary fraction and a solar metallicity of Z_⊙ = 0.02.

We extend the initial mass function (IMF) to 200M_⊙ and we keep the power-law slope for massive stars α = −2.3 (in the past we have limited IMF to 150M_⊙). This is motivated by observations of massive stars; notably, three stars in the Large Magellanic Cloud (LMC; R136a, R136b,R136c: Bestenlehner et al. 2020) and two stars in the Milky Way (WR 102ka and η Car; Hillier et al. 2001; Barniske et al. 2008) are estimated to have initial masses close to or exceeding 200M_⊙. We have also adopted favorable (in terms of forming massive BHs from stars) model from Farmer et al. (2020) and from Costa et al. (2020) that avoids PPSN mass loss for helium core masses: M_He < 90M_⊙, but allows for disruption of stars above this mass threshold. Such a model requires that carbon burning rate is decreased by 2 standard deviations and that there is an episode of dredge-up during the core-helium burning phase.

Original stellar evolution formulae that we employ in StarTack are based on stellar models only up to 50M_⊙ (Hurley et al. 2000). Our extrapolation to higher masses was checked to give reasonable results in terms of He/CO core masses and/or evolutionary tracks on Hertzsprung-Russell diagram in comparison with results obtained with detailed evolutionary calculations (e.g., with Geneva or MESA codes; Belczynski et al. 2014, 2017, 2020b). However, we note that there is no current consensus on evolution of massive stars and their calculated radii, He/CO core masses, and luminosities differ from one detailed calculation to the other (e.g., compare Tables 5 and 6 in Belczynski et al. 2020b).

The results of our model are shown in Figure 1 in which we present the dependence of the final helium core mass on the initial star mass for various metallicities. The final helium core mass is a good approximation of the BH mass for most massive stars in close binaries. The most massive stars are expected to directly collapse to BHs (Fryer 1999; Basinger et al. 2020), and stars in close binaries are typically stripped of their H-rich envelopes (BH–BH merger progenitors in particular; Belczynski et al. 2016b). Down to metallicity of Z ∼ 0.001 BH masses do not exceed M_BH ∼ 50M_⊙, which is exactly what we were obtaining with our previously employed weak mass loss from PPSN based on calculations of Leung et al. (2019). Only stars with lower metallicity (Z ∼ 0.001–0.0001) are affected by our modifications and are allowed to form BHs with very high masses M_BH ∼ 50–90M_⊙. One notes the emergence of the upper mass (at adopted M_BH = 90M_⊙) for the model with Z = 0.0001 in which BHs do not form for initial star mass above M_ZAMS > 185M_⊙.

Zoom In Zoom Out Reset image size

We follow the evolution of Population I and II (Z = 0.03–0.0001) stars with the input physics described above until the formation of BH–BH mergers. We estimate the cosmological BH–BH merger rate density using redshift-dependent star formation history and metallicity evolution across cosmic time with the standard Planck-based cosmology (Belczynski et al. 2020b). Note that we may be underestimating the amount of low-metallicity stars (Chruslinska & Nelemans 2019; Chruślińska et al. 2020) and therefore our merger rates of most massive BH–BH mergers may also be underestimated.

In Figure 2 we show an example of evolution: the formation of BH–BH merger similar to GW190521 with BH masses m₁ = 84.9M_⊙ and m₂ = 64.6M_⊙. The evolution starts with a massive primary (M_ZAMS,A = 187.1M_⊙) and a lighter secondary (M_ZAMS,B = 143.2M_⊙) at very low metallicity Z = 0.0001 on a wide (semimajor axis of a = 1247R_⊙) and virtually circular orbit (e = 0.0005). The primary star evolves off the main sequence and becomes a Hertzsprung gap star expanding and initiating a stable RLOF that increases the orbital separation (a = 2055R_⊙) and strips the primary of its H-rich envelope. The primary becomes a massive Wolf–Rayet star (M_A = 84.5M_⊙) that soon collapses directly into a BH with mass M_BH,1 = 83.6M_⊙ (no natal kick, 0.9M_⊙ mass loss in neutrinos), while the secondary is still a main sequence star. When the secondary becomes a core-helium burning star it expands over its Roche lobe and initiates a CE episode. After the CE phase the orbital separation is greatly reduced (a = 12.08R_⊙), the primary BH increases its mass through accretion in the CE (M_BH,1 = 84.9M_⊙), and the secondary loses its H-rich envelope and becomes a massive Wolf–Rayet star (M_B = 65.3M_⊙). Then the secondary star undergoes a core collapse and forms directly a second massive BH (M_BH,2 = 64.6M_⊙, no natal kick). Neutrino emission induces very a small eccentricity on the BH–BH binary (e = 0.004) and slightly expands the orbit (a = 12.08R_⊙). This BH–BH system has formed after 3.6 Myr of stellar evolution and it takes another 3.9 Myr for the two BHs to merge due to emission of gravitational radiation and associated orbital angular momentum loss. Due to the very short evolutionary and gravitational-wave emission timescale, this system would form and merge near the redshift where it has been detected (z = 0.83).

Zoom In Zoom Out Reset image size

We use the Tayler–Spruit magnetic dynamo angular momentum transport (Spruit 2002) to calculate the natal spin of the primary BH (a_spin,1 = 0.052: see Equation (4) of Belczynski et al. 2020b). This spin increases due to accretion of 1.3M_⊙ in the CE (a_spin,1 = 0.105: (MacLeod et al. 2017). The spin of the secondary BH (a_spin,2 = 0.523: see Equation (15) of Belczynski et al. 2020b) is set by the tidal spin-up of the secondary star when it is a compact Wolf–Rayet star in very close binary (a = 12R_⊙ and orbital period of P_orb = 10h). Because the BHs were formed without natal kicks, we assume that their spins are aligned with the binary angular momentum vector. This allows us to assess the effective spin parameter of this system: χ_eff = 0.29, which is within 90% credible limits of the LVC estimate (χ_eff = [−0.28: 0.35]). If the tidal spin-up were not at work as envisioned, we would calculate the secondary BH natal spin from our stellar models: a_spin,2 = 0.070 and that would have resulted in χ_eff = 0.090.

In our adopted model, massive BHs form through direct collapse of the entire progenitor star into a BH. As there is no mass loss we assume no natal kick and the system not only survives the BH formation, but also remains aligned (i.e., BH spins are aligned with binary angular momentum vector). This leads to an effective precession spin parameter equal zero (χ_p = 0) as precession requires some level of misalignment. This is apparently inconsistent with the LIGO/Virgo estimate (χ_p = [0.31: 0.93]), but this estimate is very weak (Abbott et al. 2020). Misalignment may be possibly obtained by natal kicks associated with asymmetric neutrino emission (Socrates et al. 2005; Fryer & Kusenko 2006), even if there is no baryonic mass ejection at the BH formation.

We subdivide the population of BH–BH mergers into three categories: light mergers with both BHs having mass M_BH < 50M_⊙, mixed-mass mergers with one BH with mass M_BH < 50M_⊙ and another with mass M_BH > 50M_⊙, and heavy mergers with both BHs having mass M_BH > 50M_⊙. The M_BH ≈ 50M_⊙ represents the believed (old/outdated) limit for stellar-origin BH formation set by PPSN/PSN. In Table 1 we present the merger rates of BH–BH subpopulations for a volume corresponding to redshift cuts: z = 0.1 (approximate LIGO/Virgo neutron star–neutron star (NS–NS) detection horizon), z = 0.4 (BH–NS horizon), z = 0.7 (light BH–BH horizon), z = 1.0 (mixed-mass BH–BH horizon), z = 1.5 (heavy BH–BH horizon).

Our merger-rate estimates are consistent with the 90% LVC (Abbott et al. 2019) empirical estimates: for NS–NS we find 132Gpc⁻³yr⁻¹ (LVC O1/O2: 110–3840Gpc⁻³yr⁻¹), BH–NS 11.8Gpc⁻³yr⁻¹ (LVC O1/O2:  < 610Gpc⁻³yr⁻¹), light BH–BH 63.2Gpc⁻³yr⁻¹ (LVC O1/O2: 9.7–101Gpc⁻³yr⁻¹). For heavy BH–BH mergers we find a rate of  ∼ 0.04Gpc⁻³yr⁻¹ (LVC O3: 0.02–0.43Gpc⁻³yr⁻¹ rate based on the single detection of GW190521). This may seem to be a marginal match, but note that the LVC estimates are only 90% credible limits. Merger rates are subject to change with various assumptions about input physics (natal kicks, CE, cosmic evolution of metallicity; Belczynski et al. 2020b) and they will be re-evaluated once the LVC provides more restrictive estimates.

In Figure 3 we show the intrinsic (not redshifted) distribution of the total BH–BH binary mass for mergers found in the redshift range z < 1. By construction, the light BH–BH mergers are found with M_tot = 5–100M_⊙, where the lowest masses are reached for  ∼2.5 + 2.5 mergers with both BHs originating from our delayed SN engine  (Fryer et al. 2012; Belczynski et al. 2012) and thus allowed in the lower "mass gap," while the heaviest  ∼50 + 50 mergers form with PPSN mass loss (Leung et al. 2019). The heavy mergers have total mass in the range M_tot = 100–180M_⊙, although the number of BH–BH mergers rapidly declines with increasing mass. This comes from the assumption that the IMF is steep (power law with exponent −2.3) for massive stars. In fact, the overall population of BH–BH mergers show a rapid decline of number of mergers with mass from light systems to mixed (intermediate-mass) systems to heavy systems. Note that the total BH–BH binary mass declines like an exponential (evolutionary processes affecting IMF) and not like a power law that is commonly assumed in literature. GW190521 with a total mass of ${M}_{\mathrm{tot}}={150}_{-17}^{+29}\,{M}_{\odot }$ (Abbott et al. 2020) is found in the tail of the mass distribution of our heavy BH–BH mergers. If future observations will show a flatter BH–BH mass spectrum, it would be an indication that some evolutionary process must be at work. For example, in our model the natal kicks operate only for the lightest BHs (M_BH ≲ 10–15M_⊙) and are decreasing with BH mass creating a peak in total BH mass at M_BH ∼ 20M_⊙. Had we allowed natal kicks to be applied differently it would be possible to flatten the BH mass spectrum in a desired mass range and possibly place some constraints on the core-collapse asymmetries.

Zoom In Zoom Out Reset image size

In Figure 4 we show the intrinsic mass ratio (q = M_BH,2/M_BH,1 with M_BH,1 ≥ M_BH,2) distribution of BH–BH mergers found in redshift range z < 1. The light BH–BH mergers show rather flat mass ratio distribution in a broad range q = 0.2–1 and tail reaching down to q = 0.05, with two small peaks: one at q ∼ 0.25 and another at q ∼ 0.95. The latter peak is a standard result of isolated binary evolution when a rapid SN engine (that does not produce BHs in the lower mass gap: M_BH < 5M_⊙) is applied to calculate BH mass and BH–BH mergers with similar mass BHs dominate the population (e.g., Belczynski et al. 2016b). However, note that BH–BH mergers can still reach mass ratios as small as q ∼ 0.2 (Olejak et al. 2020). The former peak, and the extent of mass ratio to very small values, is the result of our application of the delayed SN engine to calculate BH masses and our assumption that the NS/BH mass limit is at 2.5M_⊙. The population of relatively abundant (IMF) low-mass BHs (e.g., these in the lower mass gap: M_BH ∼ 2.5–5M_⊙) forms in binaries with more massive BHs creating the low-q BH–BH mergers. The lowest mass ratio arises from extreme systems with 2.5 + 50M_⊙ BH–BH mergers. Even more extreme mass ratio systems are found in BH–NS merger populations  (Drozda et al. 2020).

Zoom In Zoom Out Reset image size

The heavy BH–BH mergers are limited to q ≳ 0.6 as the lowest mass BH in this subpopulation is 50M_⊙ and the heaviest 90M_⊙. As this subpopulation does not include low-mass BHs it tends to produce similar component mass BH–BH mergers with typical mass ratio of q ∼ 0.9–1. This is consistent with LVC estimate of GW190521 mass ratio $q={0.79}_{-0.29}^{+0.19}$ Abbott et al. (2020).

We extended our evolutionary model to stars up to 200M_⊙ and we limited the action of mass loss associated with pair instabilities (Farmer et al. 2020) to test whether it is possible to form BH–BH mergers resembling GW190521 that hosts 85M_⊙ BH and 66M_⊙ BHs through classical isolated-binary evolution. Such massive BHs were/are believed not to form directly from stars.

In fact, it is possible to form massive BHs in BH–BH mergers resembling GW190521 if C-burning reaction rate uncertainties that may limit the pair-instability associated mass loss are taken into account. Once such a possibility is adopted, our standard binary evolution delivers merger rates of "normal" BHs (light BHs:  <50M_⊙) and heavy BHs ( >50M_⊙) that are consistent with LIGO/Virgo observations.

The binary evolution leading to the formation of systems resembling GW190521 is relatively simple. It requires two very massive stars (M_ZAMS ∼ 150 – 200M_⊙) at low metallicity (Z ∼ 10⁻⁴) and it involves a stable RLOF and CE episode. Our standard assumptions on BH formation involves direct BH formation through standard core collapse for both BHs with no associated PPSN mass loss and with no natal kicks.

The binary evolution leading to the formation of GW190521-like mergers may or may not involve tidal spin-up of Wolf–Rayet stars that are the immediate progenitors of massive BHs. In both cases the low predicted effective spin parameter of our proposed BH–BH merger example (χ_eff = [0.09: 0.29]) is consistent with LIGO/Virgo observations (χ_eff = [−0.28: 0.35]). In either case, the measurement of GW190521 effective spin is consistent with efficient angular momentum transport in massive stars by a magnetic dynamo.

Our model predicts that effective precession spin parameter (measuring misalignment of BH spins from binary angular momentum) for GW190521-like systems is negligible χ_p = 0. This is inconsistent with the LIGO/Virgo estimate: χ_p = [0.31: 0.93]. However, this empirical estimate was exposed as highly uncertain and a non-precessing interpretation of GW190521 cannot be excluded (Abbott et al. 2020). If precession is confirmed in such mergers it either indicates that they do not form through a classical isolated binary evolution channel, or that the second BH formation is asymmetric and leads to non-negligible BH natal kick (misalignment).

Finally, we emphasize that these new results are only valid if the carbon fusion reaction rate is highly uncertain and is allowed to be ∼2 standard deviations below the standard STARLIB rate, which is unlikely but not impossible (Farmer et al. 2020), and if during core-helium burning phase there is an episode of a dredge-up (Costa et al. 2020).

K.B. acknowledges support from the Polish National Science Center grant Maestro (2018/30/A/ST9/00050). Acknowledgment of support/advice goes to A. Olejak, P. Drozda, V. Babihav, L. Oskinova, J.-P. Lasota, R. Farmer, T. Bulik, I. Mandel, S. Woosley, M. Giersz, R. O'Shaughnessy, M. Chruslinska, and J. Klencki.