Massive Stellar Triples Leading to Sequential Binary Black Hole Mergers in the Field

The importance of interactions between massive stars in isolated binaries has become increasingly recognized in recent decades (e.g., Podsiadlowski et al. 1992; Sana et al. 2012). Recent studies indicate that early B-type and O-type stars are almost exclusively part of higher-order configurations, such as triple and quadruple systems (Sana et al. 2014; Moe & Di Stefano 2017). If future surveys confirm this, our understanding of massive stellar evolution will have to include the increased complexity of multiple-body interactions that were previously mostly considered in dense dynamical environments such as nuclear, globular, or open clusters (e.g., Sigurdsson & Hernquist 1993; Leigh & Geller 2013).

An alternative to electromagnetic methods to study high-mass stellar multiplicity is gravitational-wave (GW) astronomy, as massive stars are believed to be the progenitors of stellar-mass black holes (BHs). Stellar-mass binary black holes (BBHs) are believed to form predominantly in field binaries (e.g., Belczynski et al. 2002; Neijssel et al. 2019) and clusters (e.g., Portegies Zwart & McMillan 2000; Leigh et al. 2014; Rodriguez et al. 2019). There have been efforts to try to understand how best to segregate these two different origins, mostly based on eccentricities and spins, yet there is no definitive consensus on the origin of current GW sources (e.g., Abbott et al. 2019).

One candidate signature for the cluster origin of a GW is a BH mass in the pair-instability supernova (PISN) mass gap.¹¹ The PISNe are initiated by an electron–positron pair instability that eventually leads to explosive oxygen burning in the core of massive stars (see Langer 2012, and references therein). The PISNe do not leave behind remnants; therefore, a gap is expected in the mass distribution of BHs for stars with helium core masses in the range of ≈64–133 M_⊙ (Heger & Woosley 2002). Consequently, isolated binary evolution theory does not predict individual BHs in that regime (Stevenson et al. 2019; van Son et al. 2020). As multiple BH mergers can populate the mass gap, the discovery of BBH systems with one component lying within the mass gap has been considered a smoking gun for cluster origin (Rodriguez et al. 2019; Samsing & Hotokezaka 2020). The recent detection of GW190521 (Abbott et al. 2020a), with at least one BH within the mass gap, adds to the conundrum.

Here we give an overview of which isolated massive stellar triples experience a sequential merger of BBHs. We investigate the potential origin of such configurations, put them in the context of GW observations, and focus on the masses and spins of sequential mergers leading to BBHs in the mass gap. We propose that GW170729 is of isolated triple origin and suggest that the mass gap event GW190521 was not formed in the field. Finally, we highlight the importance of massive stellar triples in a broader astronomical context.

Here we outline our main assumptions for the key physics of the formation of a sequential merger starting at the zero-age main sequence. We consider isolated hierarchical triple systems, composed of an inner binary with masses M₁ ≥ M₂ and an outer companion of mass M₃. In the sequential mergers that we consider here, the inner BBH merges first, and afterward, the remnant merges with the outer BH. We adopt circular coplanar prograde orbits, as supported by observations of compact triples (Tokovinin 2017). For circular coplanar prograde orbits, we do not expect Lidov–Kozai cycles and neglect other three-body dynamical effects (see, e.g., Naoz 2016, for a review).

The triple must remain dynamically stable from the zero-age main sequence until the inner BBH merger, which holds if (Mardling & Aarseth 2001)

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\geqslant {\left(\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\right)}_{\mathrm{crit}}\equiv \displaystyle \frac{2.8}{1-{e}_{\mathrm{out}}}{\left[\displaystyle \frac{(1+{q}_{\mathrm{out}})(1+{e}_{\mathrm{out}})}{\sqrt{1-{e}_{\mathrm{out}}}}\right]}^{2/5},\end{eqnarray} \tag{ 1 }$$

where a is the semimajor axis, e is the eccentricity, and q is the mass ratio (the inner and outer orbits are specified by the in and out subscripts). The outer mass ratio is q_out ≡ M₃/(M₁ + M₂).

We use the synthetic BBH population from Riley et al. (2020) based on isolated binary evolution to investigate the orbital properties of the inner binary (see Appendix A for further details). The models include mass loss, mass transfer, supernovae, and chemically homogeneous evolution (CHE). In Figure 1, we present the intrinsic mass distribution of merging BBHs. We focus on low-metallicity stars (Z ≲ 10⁻³) in order to neglect mass loss, spin-down, and orbital changes due to stellar winds.

Zoom In Zoom Out Reset image size

As each star of the triple evolves, it will eventually become a BH with mass M_BH,i (with i = 1, 2, 3), for which we assume the following. The BHs have a minimum mass of ${M}_{\mathrm{BH},\min }\,=2.5\ {M}_{\odot }$. There is a mass gap in the range 43 ≲ M_gap/M_⊙ ≲ 124 (du Buisson et al. 2020) due to PISNe. The exact lower limit of the mass gap is uncertain and might be as high as M_BH ≈ 50 M_⊙ (for an overview, see Stevenson et al. 2019, and references therein). Furthermore, we assume that stars with carbon–oxygen core masses above 11 M_⊙ experience complete fallback (Fryer et al. 2012) and negligible neutrino mass loss (Müller et al. 2016), likely suppressing BH natal kicks. All BHs are born as slow rotators (Fuller & Ma 2019).

During the inner BBH merger, mass is lost by radiation from the center of mass of the merging BBH. This leads to a GW Blaauw kick similar to that in spherically symmetric supernovae (Blaauw 1961). The fraction of radiated mass with respect to the mass of the merging BBH (f_rad) depends on the masses and spins of the system (Appendix B). This modifies the orbit of the tertiary,

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{\mathrm{out},\mathrm{post}}}{{a}_{\mathrm{out},\mathrm{pre}}}={\left[2-\displaystyle \frac{{M}_{\mathrm{BH},1}+{M}_{\mathrm{BH},2}+{M}_{\mathrm{BH},3}}{{M}_{\mathrm{BBH},\mathrm{in}}+{M}_{\mathrm{BH},3}}\right]}^{-1},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{e}_{\mathrm{out},\mathrm{post}}=\displaystyle \frac{{f}_{\mathrm{rad}}({M}_{\mathrm{BH},1}+{M}_{\mathrm{BH},2})}{{M}_{\mathrm{BBH},\mathrm{in}}+{M}_{\mathrm{BH},3}},\end{eqnarray} \tag{ 3 }$$

where M_BBH,in ≈ (1 − f_rad)(M_BH,1 + M_BH,2). For nonspinning BBHs and assuming ${e}_{\mathrm{out},\mathrm{pre}}\approx 0$, f_rad ≈ 0.05, e_out,post ≈ 0.05, and ${a}_{\mathrm{out},\mathrm{post}}\approx 1.06\times {a}_{\mathrm{out},\mathrm{pre}}$. Furthermore, conservation of momentum gives rise to a recoil kick, which has a magnitude of zero when q_BBH,in ≈ 0 or 1 and can be as high as v_recoil ≲ 175 km s⁻¹ for q_BBH,in ≈ 0.36 (González et al. 2007). As our synthetic population favors q_BBH,in ≳ 0.9 (Figures 1 and 2), the magnitude of the recoil kick should be small, and we ignore it. The radiated mass fraction is maximal (${f}_{\mathrm{rad},\max }\approx \,0.12$) when we consider a merger from an equal-mass maximally spinning BBH that is aligned with the orbital spin. This constrains the postmerger eccentricity and separation to e_out ≲ 0.14 and ${a}_{\mathrm{out},\mathrm{post}}\lesssim 1.16\times {a}_{\mathrm{out},\mathrm{pre}}$, respectively. Therefore, the systems we consider here cannot be unbound during the inner BBH merger.

Zoom In Zoom Out Reset image size

We assume the outer orbit is almost circular at all times. In this case, the GW inspiral time can be approximated with (Peters 1964)

$$\begin{eqnarray}&&{T}_{c}\approx \displaystyle \frac{5}{256}\displaystyle \frac{{c}^{5}{a}_{\mathrm{out}}^{4}}{{G}^{3}{M}_{\mathrm{BBH},\mathrm{in}}{M}_{\mathrm{BH},3}({M}_{\mathrm{BBH},\mathrm{in}}+{M}_{\mathrm{BH},3})}.\end{eqnarray} \tag{ 4 }$$

We estimate the effective spin of the sequential merger as

$$\begin{eqnarray}\begin{array}{rcl}{\chi }_{\mathrm{eff}} & = & \displaystyle \frac{({M}_{\mathrm{BBH},\mathrm{in}}{{\boldsymbol{\chi }}}_{\mathrm{BBH},\mathrm{in}}+{M}_{\mathrm{BH},3}{{\boldsymbol{\chi }}}_{\mathrm{BH},3})\cdot {\hat{L}}_{{\rm{N}}}}{{M}_{\mathrm{BBH},\mathrm{in}}+{M}_{\mathrm{BH},3}}\\ & = & \ \displaystyle \frac{| {{\boldsymbol{\chi }}}_{\mathrm{BBH},\mathrm{in}}| }{1+{M}_{\mathrm{BH},3}/{M}_{\mathrm{BBH},\mathrm{in}}},\end{array}\end{eqnarray} \tag{ 5 }$$

where χ is the dimensionless component spin of the BH and ${\hat{L}}_{{\rm{N}}}$ is the unit vector parallel to the system's orbital angular momentum. For a merger of an equal-mass inner nonrotating BBH, the remnant has a spin magnitude of ∣χ_BBH,in∣ ≈ 0.68 (Boyle et al. 2008).

3.1. BH Mass Parameter Space for Sequential Mergers

3.2. Types of Triples

For the evolution of the stars, we consider standard evolving and compact constituent stars. Standard evolving stars rotate slowly, develop a composition gradient, and can expand up to thousands of solar radii during their evolution (e.g., Podsiadlowski et al. 1992). This expansion is avoided by stars in metal-free (Population III) environments (Marigo et al. 2001), in certain mass–metallicity regimes (Shenar et al. 2020), or that rotate rapidly such that rotational mixing is induced (Maeder 1987), i.e., CHE. The orbits of CHE binaries can therefore remain more compact throughout their evolution to BBH formation as compared to binaries with standard evolving stars (de Mink & Mandel 2016; Marchant et al. 2016; Riley et al. 2020). This is favorable for the sequential merger channel in order for the triple to remain dynamically stable and compact enough to lead to two mergers within a Hubble time (Equations (1) and (4)). We therefore consider four distinct types of stellar triples (illustrated in Figure 3): (1) an inner binary with compact stars and an outer standard evolving star, (2) all compact stars, (3) at least one standard evolving star in the inner binary with an outer standard evolving companion, and (4) at least one standard evolving star in the inner binary with an outer compact star.

Zoom In Zoom Out Reset image size

3.3. Populating the (PISN) Mass Gap

From the full BH mass range (Figure 1), we extract sequential mergers with total mass within the mass gap, classify them in subregions (Figure 2), discuss their evolutionary pathways (Figure 3), and present their demographics (Figure 4 and Table 1).

Zoom In Zoom Out Reset image size

Table 1.  Summary of Expected Demographics of Different Types of Sequential Mergers

Region	A-IMRI	A	U	B	C-IMRI
Rank	4	2	1	3	5
Region	Red	Red	Blue	Blue	Yellow
Triple type	1	1 and 2	1 and 2	1 and 2	3
Mass gap MBBH,in	No	Yes	Uncertain	Yes	No
MBH,3	[124, 243]	[124, 205]	[2.5, 43]	[14, 43]	[2.5, 43]
MBBH,in	[5, 43]	[43, 86]	[37, 86]	[57, 86]	[126.5, 245.5]
MBH,3 + MBBH,in	[129, 248]	[167, 248]	[80, 100]	[100, 129]	[129, 248]
qout	[0, 0.1]	[0.1, 0.69]	[0, 1]	[0.16, 0.74]	[0, 0.1]
aout	[52, 120]	[97, 135]	[35, 69]	[57, 82]	[45, 120]
χeff	[0, 0.1]	[0.1, 0.27]	[0.4, 0.68]	[0.38, 0.58]	[0.5, 0.68]
(aout/ain)∣crit	[4.8, 13.3]	[4.0, 5.7]	[2.8, 3.8]	[3.0, 3.3]	[2.8, 3.1]

Note. This table summarizes the results from Section 3 and all figures. For each quantity of interest, we present the minimum and maximum values in square brackets. We have subjectively provided a ranking for the most (1) to least (5) likely scenario to form a sequential merger based on the inner BBH mass distribution, the condition for triple stability, and the sequential merger time. Masses and separations are in solar units.

Download table as:  ASCIITypeset image

3.3.1. Red Subregion A

The predicted properties of systems in red subregion A are ${M}_{\mathrm{BBH},\mathrm{in}}\ll {M}_{\mathrm{gap},\max }\leqslant {M}_{\mathrm{BH},3}$ and 0.1 ≲ χ_eff ≲ 0.27. The evolution of an example system is illustrated in the top panel of Figure 3. Consider a CHE inner binary (triple type 1 or 2) with M₁ ≈ M₂ ≈ 40 M_⊙, R₁ ≈ R₂ ≈ 6 R_⊙, a_in ≳ 18 R_⊙. All stars in this triple will experience complete fallback, which effectively leaves the inner and outer orbits unchanged (Section 2). Assuming the outer star collapses to a BH with mass M_BH,3 ≈ 140 M_⊙, then a_out ≳ (a_out/a_in)∣_crit × a_in ≈ 4.2 × (18 R_⊙) ≈ 76 R_⊙ in order for the triple to be stable (Equation (1) and Figure 4). If the separation after the inner BBH merges is less than ${a}_{{\rm{o}}{\rm{u}}{\rm{t}},max}\lesssim 122\,{{\rm{R}}}_{\odot }$, then the sequential merger can occur within a Hubble time (Equation (4)). The critical ratio of a_out/a_in to maintain stability and the maximum orbital separation to achieve a merger within a Hubble time depend on the mass combinations of the triple and vary within a factor of a few for the combinations of interest here (Figure 4).

If the outer star was initially a standard evolving star (triple types 1 and 3), it might initiate a mass transfer phase onto the inner binary early in the evolution of the system (step A1 in Figure 3). This could occur if the radius of a standard evolving star exceeds the Roche lobe radius, i.e., R₃ ≳ 0.43 × a_in for q_out ≈ 140/80. Assuming the stellar radii approach ∼100–1000 R_⊙ at maximum (Podsiadlowski et al. 1992; Riley et al. 2020), tertiary-driven mass transfer occurs for outer orbits up to several ∼100–1000 R_⊙.

During this mass transfer phase, where the tertiary donor is significantly heavier than the inner binary, the outer orbit likely shrinks (based on angular momentum considerations), and the hydrogen envelope is stripped off the outer star. If the inner binary avoids a merger during this stage (see, e.g., Leigh et al. 2020), at the end of the mass transfer phase, the tertiary is a stripped helium star reminiscent of the compact outer star considered in triple type 2.

In this section, so far we have considered inner binaries comprised of CHE stars that remain close to each other throughout their evolution to BBHs. However, for standard evolving inner binaries (triple types 3 and 4), the orbital separation can change drastically between the zero-age main sequence and BBH formation due to mass and angular momentum exchanges during mass transfer episodes. From our synthetic population, the orbits of standard evolving binaries can remain as small as a_in ≈ 70 R_⊙ throughout their evolution, but for most systems, their orbits expand to hundreds or thousands of solar radii (see, e.g., Leigh et al. 2020). Assuming the optimistic case of a standard evolving inner binary with a_in ≈ 70 R_⊙ maximally, the outer separation must be a_out ≥ 2.8 × (70 R_⊙) ≈ 196 R_⊙ in order for the triple to be stable throughout its full evolution. Such an orbit does not merge by GW emission alone in a Hubble time (Figure 4). We conclude that in red subregion A, only triples with inner compact binaries can lead to sequential mergers within a Hubble time.

3.3.2. Blue Subregion B

The predicted properties of systems in blue subregion B are ${M}_{\mathrm{BBH},\mathrm{in}}/2\approx {M}_{\mathrm{BH},\mathrm{out}}\leqslant {M}_{\mathrm{gap},\min }$ and 0.38 ≲ χ_eff ≲ 0.58. The evolution of an example system is illustrated in the middle panel of Figure 3. Furthermore, the synthetic population suggests M_BH,1 ≈ M_BH,2 in this region; with similar BH masses, the evolutionary timescales for all component stars are similar as well. Consider a triple with M₁ ≈ M₂ ≈ M₃ ≈ 40 M_⊙. We again assume a CHE inner binary (triple type 1 or 2) with R₁ ≈ R₂ ≈ 6 R_⊙ and a_in ≳ 18 R_⊙. This triple is dynamically stable if a_out ≳ 3.3 × (18 R_⊙) ≈ 60 R_⊙. The outer separation must be ${a}_{max}\lesssim 76\,{{\rm{R}}}_{\odot }$ for the sequential merger to occur within a Hubble time due to GW emission. Hence, the small possible range for outer separations of 60 ≲ a_out/R_⊙ ≲ 76 constrains the feasibility and frequency of this subchannel. Population III BBH progenitors with initial masses M₁ ≳ 40 M_⊙ and q ≈ 1 have initial separations a_in ≳ 20 R_⊙ (Inayoshi et al. 2017), which lead to slightly more stringent constraints than for CHE binaries.

If the tertiary was a standard evolving star (triple types 1 and 3), it could initiate a mass transfer episode onto the inner binary. Due to the mass ratio, we expect this mass transfer to proceed in a stable manner. Furthermore, we expect the outer orbit to widen, typically due to angular momentum evolution, which makes the sequential merger less likely to occur within a Hubble time (see, e.g., Belczynski et al. 2002, for an overview in the context of BBH progenitors). In summary, in blue subregion B, only triple compact binaries in a fine-tuned configuration can lead to sequential mergers within a Hubble time.

3.3.3. Yellow Subregion C

The masses of systems in yellow subregion C are ${M}_{\mathrm{BH},3}\,\ll {M}_{\mathrm{gap},\max }\leqslant {M}_{\mathrm{BBH},\mathrm{in}}$. The evolution of an example system is illustrated in the bottom panel of Figure 3. Even though our synthetic binary population does not predict any inner BBHs in this region (Appendix), BBHs above the mass gap have been suggested for both compact (Marchant et al. 2016; du Buisson et al. 2020) and standard evolving (Mangiagli et al. 2019) binaries. Regarding the former, as inner binaries experiencing CHE (triple types 1 or 2) are expected to have mass ratios q_in ≈ 1 and ${M}_{\mathrm{BBH},\mathrm{in}}\gtrsim 2\times {M}_{\mathrm{gap},\max }\approx 248\ {M}_{\odot }$, these systems can only bring about sequential mergers with a total mass above the mass gap. Regarding the latter, the case of hypothetical standard evolving binaries is also not promising for sequential mergers. If mass transfer took place from the lower-mass tertiary companion to the heavier-mass inner binary, it would likely be stable and widen the outer orbit. This would increase the GW inspiral time for the outer orbit. Therefore, we do not consider sequential mergers from subregion C to be common.

3.4. GW170729 as a Sequential Merger

3.5. GW190521: A BBH in the Mass Gap

4.1. Main Caveats to Our Method

4.2. Stellar Triples: Birth Distributions and Orbital Evolution

4.3. Uncertainties in the (PISN) Mass Gap

4.4. Multiple BH Mergers

4.5. Intermediate Mass-ratio Inspirals

4.6. Rate of Sequential Mergers

We investigated configurations of isolated massive stellar triples that lead to sequential BBH mergers. We find that triples with CHE inner binaries are good candidates for sequential mergers. Our model predicts that GW sources with one BH in the mass gap and χ_eff > 0.1 can be of sequential merger origin. We highlight two classes of triples that lead to BBHs in the mass gap. The first one has a tertiary BH above the mass gap and 0.1 ≲ χ_eff ≲ 0.27 (see red subregion A in Figures 2 and 4). The second one has a tertiary BH below the mass gap and 0.38 ≲ χ_eff ≲ 0.58 (see blue subregion B in Figures 2 and 4). We suggest that GW170729 is of triple origin and belongs to the second class. The masses and spins of mass gap event GW190521 are inconsistent with our model, which further supports a cluster origin.

From a broader point of view, we outlined a new outcome of the evolution of massive stellar triples from a proof-of-principle study. To improve upon the predictions made here, several processes should be considered: the effects of the uncertain initial orbital configurations of stellar triples, which may also lead to nonnegligible three-body dynamical effects; the nontrivial problem of mass transfer in triples; GW recoil kicks from the inner BBH merger; and their combined effects on population statistics. Higher-order multiplicity in massive stars is crucial to understanding the most energetic astronomical phenomena in the Universe.

The authors thank David R. Aguilera-Dena, Chris Belczynski, Giacomo Fragione, Ryosuke Hirai, Ilya Mandel, Chris Moore, Ataru Tanikawa, and Team COMPAS for useful discussions and suggestions. We thank Teresa Rebagliato for the graphic representation of the formation channels. A.V.-G. and E.R.-R. acknowledge funding support by the Danish National Research Foundation (DNRF132). E.R.-R. is also supported by the Heising-Simons Foundation and NSF (AST-1911206 and AST-1852393). S.T. acknowledges support from the Netherlands Research Council NWO (VENI 639.041.645 grants). N.W.C.L. gratefully acknowledges support from the Chilean government via Fondecyt Iniciación grant No. 11180005. C.-J.H. acknowledges the support of the National Science Foundation and the LIGO Laboratory. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-1764464.

Software: COMPAS¹¹ v02.11.04, publicly available at GitHub via TeamCOMPAS/COMPAS.¹² Scripts used for this study are available in GitHub via avigna/sequential-mergers.¹³

We use the publicly available data from Riley et al. (2020) for synthetic BBH formation and merger rates. That study made use of the COMPAS rapid population synthesis code (Stevenson et al. 2017; Vigna-Gómez et al. 2018; Neijssel et al. 2019), which is freely available at http://github.com/TeamCOMPAS/COMPAS. The version of COMPAS used for these simulations was v02.11.01a, built specifically for Riley et al. (2020); functionality in this release was integrated into the public COMPAS code base in v02.11.04.

Riley et al. (2020) performed, for the first time, simultaneous population synthesis of CHE and standard evolving binaries. We use the orbital distributions of BBH mergers from that study as an educated guess for the properties of the inner BBH in the sequential merger scenario. While this works well as a first approach, the initial conditions from Riley et al. (2020) are based on birth distributions from observations of isolated massive binaries (Sana et al. 2012). Riley et al. (2020) simulated binaries with metallicities $-4\leqslant {\mathrm{log}}_{10}Z\leqslant -1.825$. As there are no binary evolution models available at zero-metallicity Population III stars, all of our quantitative results from compact binaries come exclusively from CHE binaries.

The most relevant conclusions drawn from this data concerning sequential mergers are as follows.

• 1.  
  There are, overall, ∼4× more merging BBHs from standard evolving than from CHE binaries. The local merger rate (z = 0), which accounts for star formation history and galaxy mass–metallicity dependence, is 50 and 20 Gpc⁻³ yr⁻¹ for isolated binaries and the subset of CHE binaries, respectively. For CHE binaries, we include massive overcontact binaries. A massive overcontact binary can fill its inner Lagrangian point during the main sequence as long as the outer Lagrangian point is not filled (Marchant et al. 2016).
• 2.  
  The yield of merging BBHs from CHE binaries is roughly constant below ${\mathrm{log}}_{10}Z/{Z}_{\odot }\lesssim -0.5$ (see Figure 6 of Riley et al. 2020), with Z_⊙ = 0.0142. Low-metallicity CHE stars are more compact than high-metallicity ones; they also experience less mass loss and orbital widening through stellar winds.
• 3.  
  The maximum premerger total mass of BBHs for both CHE and standard evolving binaries is M_BH,1 + M_BH,2 ≈ 79 M_⊙.
• 4.  
  In this COMPAS data set, there are no BBHs above the mass gap. This holds even when including stars with initial masses up to 150 M_⊙, for which the adopted stellar evolution models from Hurley et al. (2000) are extrapolated to stars above 50 M_⊙. However, we consider that there is a possibility of BBHs above the mass gap in nature. They could come from standard evolving stars (e.g., Mangiagli et al. 2019) or CHE binaries (e.g., Marchant et al. 2016; du Buisson et al. 2020).
• 5.  
  We use the delayed supernova mechanism prescription from Fryer et al. (2012) to determine the remnant mass and natal kick distributions. This prescription predicts complete fallback for BH progenitors with carbon–oxygen core masses above 11 M_⊙. In our population, more than half of standard evolving stars and all CHE stars leading to BBHs have masses above this threshold. Complete fallback has two implications. The first is that the final baryonic mass of the remnant is the same as the presupernova mass (Equations (19) and (20) from Fryer et al. 2012) modulo neutrino emission. The second is that there are no natal kicks for heavy BHs (Equation (21) of Fryer et al. 2012).

We emphasize that a full population synthesis of massive stellar triples will help us better understand the evolution of sequential mergers and their delay time distribution and constrain their rates. We hope future software developments and observations will make this sort of study possible in the near future.

For a BBH, f_rad is the amount of energy radiated away in the form of GWs during the coalescence, expressed as a fraction of the total mass of the BBH. While the magnitude of f_rad is independent of the BBH total mass, it does depend on the binary mass ratio q_in = M_BH,2/M_BH,1 (with M_BH,1 ≥ M_BH,2) and the BH spin configuration. More specifically, it depends on ${\chi }_{1,2}^{{L}^{| | }}\equiv {\chi }_{\mathrm{1,2}}\cos ({\theta }_{\mathrm{1,2}})$, with 0 ≤ χ_1,2 ≤ 1 being the BH spin magnitude and θ_1,2 the zenith angle between the spin and orbital angular momenta at the time of merger for each BH. The f_rad is typically approximated by making use of fitting formulae that are based on numerical relativity simulations of BBHs.

In this study, we use the prescription from Equation (28) of Jiménez-Forteza et al. (2017) as implemented in LALSuite (LIGO Scientific Collaboration 2018) and show f_rad for a selection of BH spins and mass ratios in Figure 5. For simplicity, and following our assumptions, we have restricted the figure to binaries with equal ${\chi }_{1,2}^{{L}^{| | }}$ values. We use this prescription to estimate the final masses from both the inner BBH merger (Figures 1 and 2) and sequential mergers.

Zoom In Zoom Out Reset image size

In order to test the validity of the sequential merger scenario for Abbott et al. (2019), we focus on the expected delay times for our scenario. The delay time is the time it takes a system to evolve from the zero-age main sequence to the sequential BBH merger and follows from our models. We will convert it into a redshift in the following way. We use astropy's cosmology module to convert between redshift and look-back time using astropys lookback_time (Astropy Collaboration et al. 2013, 2018). We assume a flat ΛCDM cosmology with H₀ = 67.90 km s⁻¹ Mpc⁻¹ and Ω_m = 0.3065 following Planck Collaboration et al. (2016). With this setup, for a delay time of ∼5.1 Gyr, we estimate a redshift z ≈ 0.49.

In Figure 2, we mark the region of sequential mergers with masses M ≥ 124 M_⊙ and mass ratios q < 0.2 as being IMRIs. This type of system is degenerate with IMRIs of binary origin. The IMRIs have primary masses of, e.g., M_BH,1 ≥ 100 M_⊙ and mass ratios of, e.g., q < 0.1 (see, e.g., Haster et al. 2016a, 2016b). In Section 4.5, we briefly discuss how they are usually associated with a cluster origin, but we do not rule out an isolated origin. In Table 1, we present the properties of IMRIs from red subregion A and yellow subregion C.

In this Appendix, we briefly expand our analysis of sequential mergers leading to systems that are degenerate in mass with IMRIs. Figure 6 considers double compact objects with primary masses 124 ≤ M_BH,1/M_⊙ ≤ 500 and light compact object masses 1 ≤ M_light/M_⊙ ≤ 43. In this case, light compact objects include both neutron stars and BHs. It is not dependent on the maximum mass of neutron stars or the existence or absence of a mass gap between neutron stars and BHs. These limits constrain the mass ratio to 0 ≲ q ≲ 0.35. To date, GW190814 is the GW source with the most extreme mass ratio of q = M_light/M_BH,1 ≈ 2.6/23 ≈ 0.11 (Abbott et al. 2020c). However, the individual masses of GW190814 are well below those of typical IMRIs.

Zoom In Zoom Out Reset image size

For an IMRI, the inference of the effective spin might not be enough to classify its origin. For a sequential merger, if the tertiary is the most massive component, then the effective spin is dominated by this massive nonrotating BH, i.e., χ_eff ≈ 0. This value would be similar to that of a nonrotating binary. On the other hand, if the tertiary is the least massive component, that effective spin is χ_eff ≈ 0.7. This is summarized visually in Figure 4 and quantitatively in Table 1.

More detailed studies of the spin properties of BBH mergers from isolated binary, isolated triple, and cluster origin would help in the classification of future detections of IMRIs. Additionally, population studies might be helpful in constraining the merger rates from different formation channels.

In order to give an upper limit in the local (z = 0) rate of sequential mergers (${ \mathcal R }$), we make a simple estimate in the form

$$\begin{eqnarray}&&{ \mathcal R }={{ \mathcal R }}_{\mathrm{BBH},\mathrm{in}}\times ({f}_{\mathrm{triple}}/{f}_{\mathrm{binary}})\times {f}_{\mathrm{separation}}.\end{eqnarray} \tag{ E1 }$$

We estimate the BBH merger rate ${{ \mathcal R }}_{\mathrm{BBH},\mathrm{in}}=50\ {\mathrm{Gpc}}^{-3}\ {\mathrm{yr}}^{-1}$ based on the isolated binary evolution calculations of Riley et al. (2020), which account for both compact and standard evolving stars. This optimistic rate includes sequential mergers with a total mass below the mass gap. Riley et al. (2020) assumed that a fraction f_binary = 0.7 of massive star systems are binaries. The factor f_triple = 0.35 describes the fraction of massive star systems that are stellar triples (Moe & Di Stefano 2017).

The factor f_separation accounts for the fraction of systems in which the outer orbital can merge within a Hubble time due to GW emission. To estimate f_separation, we assume that the log of the outer birth orbital separation is distributed uniformly $p({a}_{\mathrm{out}})\propto {a}_{\mathrm{out}}^{-1}$ (Öpik 1924; Abt 1983). We then estimate the fraction of systems of interest in the form

$$\begin{eqnarray}&&{f}_{{\rm{s}}{\rm{e}}{\rm{p}}{\rm{a}}{\rm{r}}{\rm{a}}{\rm{t}}{\rm{i}}{\rm{o}}{\rm{n}}}=\displaystyle \frac{{\int }_{{a}_{{\rm{o}}{\rm{u}}{\rm{t}}{\textstyle \text{-}}{\rm{S}}{\rm{M}},min}}^{{a}_{{\rm{o}}{\rm{u}}{\rm{t}}{\textstyle \text{-}}{\rm{S}}{\rm{M}},max}}{a}^{-1}da}{{\int }_{{a}_{{\rm{o}}{\rm{u}}{\rm{t}},min}}^{{a}_{{\rm{o}}{\rm{u}}{\rm{t}},max}}{a}^{-1}da}=\displaystyle \frac{{\int }_{35\,{{\rm{R}}}_{\odot }}^{135\,{{\rm{R}}}_{\odot }}{a}^{-1}da}{{\int }_{28\,{{\rm{R}}}_{\odot }}^{2\times {10}^{6}\,{{\rm{R}}}_{\odot }}{a}^{-1}da}\approx 0.12,\end{eqnarray} \tag{ E2 }$$

where a_out is the outer separation of stable triples, and a_out‐SM is the outer separation of potential sequential mergers. The lower limit on a_out is given by the smallest value for a stable triple with an inner binary separation a_in ≈ 10 R_⊙, ${({a}_{\mathrm{out}}/{a}_{\mathrm{in}})}_{\mathrm{crit}}=2.8$, and therefore ${a}_{{\rm{o}}{\rm{u}}{\rm{t}},min}=28\,{{\rm{R}}}_{\odot }$. The maximum outer separation at birth that we consider is ${a}_{\mathrm{out},\max }\approx {10}^{4}\ \mathrm{au}$ ≈ 2 × 10⁶ R_⊙ (Moe & Di Stefano 2017). For the outer separation of potential sequential mergers, following the results of Table 1, as shown in top right panel of Figure 4, the limits are 35 ≲ a_out‐SM/R_⊙ ≲ 135. These limits in a_out‐SM assume that the separation does not drastically change from the zero-age main sequence until the inner BBH merger. For some triples, this assumption will not necessarily hold (Section 4). When it holds, we expect the distribution of merger times to be p(t) ∝ t⁻¹, as expected from GW-dominated binaries with a flat-in-the-log distribution at BBH formation. This assumption also neglects initially wider tertiary stars that are stripped and become potential sequential mergers (Section 3.3.1), which would increase the values of ${a}_{\mathrm{out} \mbox{-} \mathrm{SM},\max }$, f_separation, and, ultimately, the rate. The assumption of different orbital initial distributions would naturally also affect the rates. We choose these assumptions for an order-of-magnitude estimate and leave a more thorough analysis of the orbital parameter space and the distribution of merger times for a future study.

After substituting all of the estimated parameters in Equation (E1), we constrain the upper limit to ${ \mathcal R }\,\lt 3\ {\mathrm{Gpc}}^{-3}\ {\mathrm{yr}}^{-1}$ for sequential mergers.

As a proof of concept of the evolution toward sequential mergers, we have simulated the evolution of the proposed progenitor of GW170729 with two independent codes: the rapid binary population synthesis code COMPAS and the triple evolution code TRES (Toonen et al. 2016, 2020). We use these codes to model the evolution from the zero-age main sequence to the inner BBH merger.

We use COMPAS as described in Appendix A to explore a binary system that is representative of the CHE inner binary for a GW170729-like system. We follow the evolution from the zero-age main sequence until BBH formation. This system initially consists of a circular CHE binary with M₁ = M₂ = 29 M_⊙, R₁ = R₂ = 4.8 R_⊙, and a = 10.2 R_⊙. The separation, masses, and radial time evolution are shown in Figure 7. The CHEs are expected to contract throughout their lifetimes (see, e.g., Aguilera-Dena et al. 2018); however, the current implementation of them in COMPAS assumes a fixed radius during the main sequence. The low metallicity results in negligible mass loss, and the orbit barely changes throughout the full evolution. There are two milestones in the evolution. The first one is the evolution after hydrogen depletion (≈6.8 Myr), when each component becomes a 1.8 R_⊙ naked helium star. The final one is the failed supernova (≈7.1 Myr), which reduces the mass by 10% (Appendix A) and increases the separation to a = 11.6 R_⊙ and the eccentricity to e = 0.07. This BBH merges in ≈76 Myr (Peters 1964).

Zoom In Zoom Out Reset image size

Additionally, we modified the triple evolution code TRES (Toonen et al. 2016, 2020) to model the evolution from the zero-age main sequence until the BBH coalescence of a compact inner binary of a triple system. The inner binary consists of a CHE system with M₁ = M₂ = 29 M_⊙, R₁ = R₂ = 4.8 R_⊙, a_in = 10.3 R_⊙, and e_in = 10⁻⁵ (this eccentricity is the lower limit allowed in TRES for numerical reasons). Similarly to COMPAS, the radial evolution of all compact stars throughout the main sequence is kept constant. The tertiary is also assumed to be a compact star with M₃ = 33 M_⊙, R₃ = 5.2 R_⊙ in an orbit of a_out = 45 R_⊙, e_out = 10⁻⁵, and relative inclination i = 0.0. The complete orbital evolution is shown in Figure 8. There are four milestones in the evolution of this triple system. The first is the formation of an outer 33.4 M_⊙ BH at ≈6 Myr. The second is the evolution after hydrogen depletion (≈6.8 Myr), when each component of the inner binary becomes a 1.8 R_⊙ naked helium star. The third is the supernovae (≈7.1 Myr) of the stars in the inner binary, which result in two BHs of M_BH,1 = M_BH,2 = 28.5 M_⊙, a_in = 10.4 R_⊙, e_in ≈ 2 × 10⁻⁴, e_out ≈ 0.8 × 10⁻⁴, and i ≈ 10⁻⁵. The final milestone is the GW evolution from the inner binary; the inner BBH merges at ≈46 Myr. The triple remains approximately circular and coplanar throughout the evolution, validating our assumptions at the moment of triple BH formation and throughout the sequential merger.

Zoom In Zoom Out Reset image size

In summary, we used two independent codes to test the validity of our assumptions in Section 2. With COMPAS, we validate the compact inner binary evolution. We use TRES to do a dynamical calculation to corroborate that, following our initial conditions, there are no major orbital changes that would modify our assumptions during both BBH mergers. Finally, we show that our numerical results are consistent with our semianalytical formalism.