1 Introduction↩︎

Gravitational-wave LIGO-Virgo-KAGRA (LVK) science runs O1–O4 detected about 300 black hole (BH)-BH merger events (e.g., [1]). These detections promoted a huge theoretical study of the formation channels of binary BH (BBH) progenitors. Most channels fall into two broad categories (e.g., [2] for a review). (1) Isolated binary systems where the merging BBH interact without the influence of a third body, either via a common envelope evolution (CEE; e.g., [3]–[7]), or without a CEE (e.g., [8]–[12]); see [13] for a recent study. (2) The dynamical interaction channels, where a pure dynamical interaction with a third star or stars in a binary BH bound system (e.g., [14]–[19]), or with stars and ambient gas (like the disk around a supermassive black hole at the center of galaxies) of two initially unbound BHs (e.g., [20]–[24]). Note that although CEE is emphasized here, the prevalence of CEE among BBH progenitors, and indeed whether a CEE phase is required at all to form merging BBHs, remains an open question, as merger-rate constraints alone do not fix the formation pathway ([25]). Yet, there is a third channel: a triple-star system in which the main interaction is via mass transfer and CEE in a triple (or higher multiple) system, rather than a pure dynamical interaction (e.g., [26]).

In this study, we address the effective spin of the merging BBHs \[\chi_{\rm eff} \equiv \frac{M_1 \vec{\chi}_1 + M_2 \vec{\chi}_2}{M_1 + M_2} \cdot \frac{\vec{L}}{\vert \vec{L} \vert} , \label{eq:EffectiveSpin}\tag{1}\] where \(M_1\) and \(M_2\) are the masses of the two BHs, \(\chi_1=c J_1/G M^2_1\) and \(\chi_2=c J_2/GM^2_2\) are the dimensionless spins of the two BHs, \(J_1\) and \(J_2\) are the spins of the two BHs, and \(\vec{L}\) is the angular momentum of the orbital motion of the BBH. The observed distribution of \(\chi_{\rm eff}\) is asymmetric, with more systems having \(\chi_{\rm eff}>0\) than those having \(\chi_{\rm eff}<0\) (e.g., [1], [27], [28]; see Figure 1).

The asymmetric distribution shows that the orbital motion influences the final BH spins and disfavors a purely dynamical origin (e.g., [29]). The significant number of systems with \(\chi_{\rm eff}<0\) shows that other effects determine the BH spins in addition to the orbital angular momentum. The dynamical channels, where the spins of the two BHs need not be correlated with the final angular momentum of the two merging BHs, predict a large number of systems with \(\chi_{\rm eff} \leq 0\). However, they do not account for the large asymmetric distribution of \(\chi_{\rm eff}\). The expectation that dynamical assembly isotropizes the spins has, however, recently been questioned: [30] find that BBHs entering with aligned spins can retain that alignment over several encounters, though this inherited alignment adds to the positive \(\chi_{\mathrm{eff}}\) side rather than the negative tail.

One of several processes might impart negative effective spins to the BHs in the isolated binary system channels, without involving pure dynamical effects.

1. A very large natal velocity kick to the second-born BH (e.g.,

   [31]) such that the orbit is flipped. Other observations suggest that BHs do not often acquire the required large kick velocity of \(> 150 {~\rm km}{~\rm s}^{-1}\), although some might (e.g., [32]).

2. The tossing of the BH spin during its formation (e.g., [33], [34]). [34] find that BH natal kick alone cannot explain the \(\chi_{\rm eff}\) distribution. Either there is a need for spin-axis tossing, or \(72 \pm 8\%\) of merging BBHs result from dynamical interaction channels.

3. When a triple stellar system experiences two CEE phases, or one CEE and one mass transfer phase ([26]). [26] considered the entrance of the tight binary system, including one or two BHs, into the giant envelope of a red supergiant (RSG) of a third star. Although this is a triple system interaction, it is not purely dynamical (hence it does not belong to the dynamical interaction channels, and does not involve the von-Zeipel-Lidov-Kozai instability). The channel we study here differs from those studied by [26], but it does involve double-CEE evolution.

The \(\chi_{\rm eff}\) distribution and spin-orbit misalignment of BBHs are open questions that force different scenarios to make specific assumptions, either about the physical processes or about the fraction of BBHs arising from different channels. [35] conclude that only dynamical processes can form significant spin-orbit misalignment in BBH mergers. [36] studied the spin of BBH merger in the CEE channel of binary systems, assuming that the spin of the first-born BH is very low. They find most of their BBH to merge with the spin parameter of \(\chi_{\rm eff} \simeq 0\), and a few with positive \(\chi_{\rm eff}\). [37] consider that BHs are born with very low spin in the framework of the isolated binary channel. [4] find that the CEE requires that the first-born BH enter a CEE with an RSG. [38] argued that a mixture of channels is strongly preferred over any single channel of BBH formation. [20] find from their analysis that multiple channels contribute to the BBH population. They did not study the origin of spins in isolated BBH, but rather introduced a random distribution as an input parameter. Recent data-driven analyses of the GWTC-4 catalog likewise favour more than one subpopulation: [39] find the strongest evidence for two channels, with only mild support for an additional AGN-disc component, while [40] identify three subpopulations in the delay-time distribution whose local merger rates differ by about an order of magnitude.

In this study, we consider an isolated triple-star channel within the jittering-jets explosion mechanism (JJEM) of core-collapse supernovae (CCSNe). [41] study the formation of BBH taking their masses in the framework of the neutrino-driven mechanism. We, on the other hand, build our scenario within the JJEM framework. We describe this scenario for BBH formation in the framework of the JJEM and our assumptions in Section 2. We then estimate the formation rate of such systems relative to that of all CCSNe; we describe the calculation method and the relevant assumptions in Section 3, and the results in Section 4. We summarize our results and discuss the formation of misaligned BBHs in the context of the scenario we study in Section 5.

2 Basic assumptions of the scenario↩︎

We conduct this study in the framework of the JJEM, where a CCSN leads to the formation of a BH when the progenitor is a massive star with a rapidly rotating core. The reason is that when the pre-collapse core rotation is fast, the accretion disk around the newly born neutron star, and then the BH, has a roughly constant axis; some small jittering or precession around the pre-collapse rotation axes is possible. Jets with a fixed axis are less efficient at expelling core material from the equatorial plane vicinity (the disk plane) than jittering jets in different directions (e.g., [42]). If, in addition, the binding energy of the core is large, which implies a massive progenitor, the falling material accumulates to form a BH. This event can lead to a very energetic CCSN, rather than to a failed CCSN (e.g., [42], [43]).

Our basic assumption, which is fully compatible with the JJEM, is that for the formation of a BH, the progenitor should be massive, a zero-age main-sequence (ZAMS) mass of \(M_{\rm ZAMS} \gtrsim 20 M_\odot\), and its inner regions should be spun-up during its post-main-sequence evolution. Close binary interaction is required to spin-up the pre-collapse core to fast rotation (e.g., [44]).

We consider a system of two very massive main-sequence stars, and a tertiary star of mass \(M_3 \lesssim 5 M_\odot\), orbiting the more massive of the two. The evolutionary scenario is depicted schematically in Figure 2, and the corresponding angular momentum geometry is illustrated in Figure 3. In the scenario we study, this more massive star with ZAMS mass of \(M_{\rm 1,ZAMS}\), which is the progenitor of the first BH to be born, engulfs the tertiary star as it evolves off the main-sequence. This forms the first CEE phase of the system. The tertiary star spirals all the way to the core vicinity and spins up the core. The core destroys the tertiary star by tidal disruption. The explosion of this star forms a BH with its spin axis \(\vec{\chi}_1\), more or less along the orbital angular momentum \(\vec{L}_{\mathrm{orb,in}}\) of the primary and tertiary (inner binary) system. Mass transfer from the RSG primary to the secondary is likely to occur.

When the secondary star, of ZAMS mass \(M_{\rm 2,ZAMS} < M_{\rm 1,ZAMS}\) evolves to become an RSG star, it engulfs the first BH, now of mass \(M_1\). This is the second CEE phase of the system. Note that because of mass transfer from the primary to the secondary, the secondary mass at this stage is likely to be \(M_2>M_{\rm 2,ZAMS}\), and might even be \(M_2>M_{\rm 1,ZAMS}\). The BH spirals close to the core of the secondary star and spins up the zones around the secondary stellar core. The massive, rapidly rotating core forms a BH according to the JJEM; this is the second-born BH. After the explosion of the secondary star, it leaves a remnant, a BH with a spin axis more or less along the orbital angular momentum axis of the binary system.

The first-born BH is likely to accrete some mass from the envelope. The angular momentum of this accretion process is aligned with the orbital angular momentum of the binary system. Therefore, the accretion process in the second CEE phase decreases the angle between \(\vec{\chi}_1\) and that of the binary orbital angular momentum \(\vec{L}_{\mathrm{orb}}\). [45], who studied binary CEE, argued that the first-born BH has a small spin (\(<0.1\)). We follow them and assume that the change in the spin of the first-born BH during the CEE is small, although in our scenario its spin does not need to be small.

Like in earlier studies of similar evolutionary routes (e.g., [45], [46]), our scenario predicts that the spin of the second-born BH, \(\vec{\chi} _2\), tends to be aligned with the orbital angular momentum of the BBH, \(\vec{L}_{\mathrm{orb}}\), and might be large. However, several processes can cause misalignment: (1) A natal kick velocity of the second-born BH (see Section 1) causes the new orbital angular momentum axis, \(\hat{L}\), to differ from the axis of the orbital angular momentum before the explosion, \(\hat{L} \ne \hat{L}_{\mathrm{orb}}\). (2) Accretion of stochastic angular momentum from the envelope, as the JJEM implies, might cause small changes in the spin axis \(\chi_1\). (3) A fourth star in the system that orbits the secondary. The fourth star must be light in order not to prevent the secondary from becoming an RSG. The three processes can cause a misalignment, but not a large one.

To summarize, the scenario we study in the framework of the JJEM requires a triple-star system to form a BBH. The scenario allows for \(\chi_{\rm eff} < 0\), and predicts that it results from the misalignment of the first-born BH, i.e., the axes of \(\vec{\chi}_1\) and \(\vec{L}_{\mathrm{orb}}\) can be highly misaligned, and even opposite. The second-born BH spin tends to be aligned with the orbital angular momentum axis. It also predicts an asymmetrical distribution of \(\chi_{\rm eff}\), as it is much more likely to have a positive than a negative value.

3 The numerical settings↩︎

3.1 Outline of the semi-analytical approach↩︎

Our goal is to assess whether the triple-star channel of Section 2 can plausibly contribute to the population of merging BBHs, rather than to quantify its precise contribution. Quantifying the contribution would require a population-synthesis treatment of triple systems, which is computationally expensive; we instead adopt a semi-analytical estimate that is sufficient to test whether the channel is a viable contributor.

Because not every triple is born with the properties needed to begin, let alone complete, this evolutionary pathway, the calculation separates into two ingredients. (i) A draw probability \(P_{\rm draw}\): the probability that a zero-age main-sequence (ZAMS) system is born with parameters consistent with the configuration the channel requires. (ii) A success fraction \(\beta\): the conditional probability that a system satisfying these ZAMS conditions evolves through the full sequence - passes two CEE phases, two CCSNe each forming a BH, and a final orbit tight enough to merge within a Hubble time. The factor \(\beta\) absorbs the cumulative uncertainties of the scenario, such as the CEE ejection efficiency, the natal kicks, the mass loss, and the requirement of a sufficiently tight post-collapse orbit. The logical structure of the calculation is summarized schematically in Figure 4.

We normalize the calculation to the CCSN rate - an observed quantity. The local BBH merger-rate density is then \[\mathcal{R}_{\rm BBH} = \frac{\mathcal{R}_{\rm CCSN}}{f}\, P_{\rm draw}\, \beta , \label{eq:rate}\tag{2}\] where \(\mathcal{R}_{\rm CCSN}\) is the volumetric CCSN rate and \(f\) is the mean number of stars massive enough to undergo a CCSN per massive-primary system. The factor \(f\) converts the event rate \(\mathcal{R}_{\rm CCSN}\) into a rate of progenitor systems: because massive stars reside predominantly in binaries and higher-order multiples ([47], [48]), a system whose primary explodes as a CCSN typically hosts at least one further companion massive enough (\(M \gtrsim 8\,M_\odot\)) to do the same. In the channel studied here, both the primary and the secondary collapse to BHs, and averaging the companion statistics over the primary mass function gives \(f \simeq 2\). We adopt a local CCSNe rate of \(\mathcal{R}_{\rm CCSN} \simeq 6.9 \times 10^{4}~{\rm Gpc^{-3}\,yr^{-1}}\) ([49]). [28] gives a BBH merger rate of \(27.5-49.4~{\rm Gpc^{-3}\,yr^{-1}}\) at \(z=0.2\), we convert this value to \(z=0\) to ensure consistency with the local CCSNe rate (see appendix 6). This conversion yields a local BBH merger rate of \(\mathcal{R}_{\rm BBH} \simeq 17.4-31.3~{\rm Gpc^{-3}\,yr^{-1}}\), which is consistent with the local GWTC-4.0 value of \(14-26~{\rm Gpc^{-3}\,yr^{-1}}\) ([1]). Using these values results in a rate ratio of \(\mathcal{R}_{\rm BBH}/\mathcal{R}_{\rm CCSN} \simeq 3.4 \times 10^{-4}\), allowing equation 2 to be rearranged into the form: \[\beta = f\,\frac{\mathcal{R}_{\rm BBH}/\mathcal{R}_{\rm CCSN}}{P_{\rm draw}} \simeq \frac{2 \times \left(3.4 \times 10^{-4}\right)}{P_{\rm draw}} . \label{eq:beta}\tag{3}\]

This relation is the central test of the channel: for a computed \(P_{\rm draw}\), the inferred \(\beta\) must satisfy \(0 < \beta < 1\) to be interpretable as a conditional probability. A value exceeding unity would mean that the adopted ZAMS configuration is too restrictive to accommodate the observed BBH rate through this channel alone, whereas a value comfortably below unity indicates that the channel can supply a meaningful share of the rate without requiring an implausibly large fraction of progenitor systems to succeed.

3.2 ZAMS distributions and the draw probability↩︎

To evaluate the draw probability, we specify the joint ZAMS distribution \(f\) and the region of parameter space over which it is integrated, \[P_{\rm draw} = \int_{\mathcal{R}} f(M_1, \log P, q, e)\; dM_1\, d\log P\, dq\, de . \label{eq:Pdraw}\tag{4}\] where \(\mathcal{R}\) is the region picked out by the cuts on these parameters. The specific cuts and their physical justifications are presented alongside the results in section 4. For consistency, all subsequent mass and temporal quantities are reported in solar masses (\(M_\odot\)) and days, respectively.

For the joint distribution, we adopt the empirically calibrated massive-star statistics of [48], who compile and de-bias multiple observational samples and provide joint fits of primary mass, orbital period, mass ratio, and eccentricity suitable for population-synthesis initialization. Following [48] and [50], we represent the joint distribution using conditional factors that mirror the structure of the empirical fits, \[\begin{align} f(M_1, & \log P, q, e) \simeq f(M_1) f(\log P \mid M_1) \\ & \times f(q \mid \log P, M_1) f(e \mid \log P, M_1) \\ & \times \Pi_{\mathrm{occ}}(M_1, \log P) \end{align} \label{eq:prob95expression}\tag{5}\] where \(\Pi_{\mathrm{occ}}\) is a companion occurrence factor - the probability that a system with primary mass \(M_1\) has at least one companion in the period range being sampled.

The primary mass follows the standard initial mass function ([51]) described in the high mass regime by the power law \(\xi(M) \propto M^{-2.3}\). To remain consistent with the CCSN normalization of equation 2 , the parent population is defined as all stars capable of core collapse (\(M_1 \geq 8\,M_\odot\)), so every reported draw probability is conditioned on \(M_1 \geq 8\,M_\odot\).

Because we impose no eccentricity cut beyond the requirement that the system not overflow its Roche lobe at periastron, the eccentricity factor \(f(e \mid \log P, M_1)\) evaluates to unity over the sampled range. The [48] mass-ratio calibration is a continuous three-regime piecewise fit defined for \(q \geq 0.1\); since our scenario requires a low-mass tertiary (Section 2) whose mass ratio \(M_3/M_1\) can fall below this floor, we extrapolate the distribution down to \(q_{\mathrm{min}} = 0.05\) for that companion, and find the resulting change in \(P_{\rm draw}\) to be small.

We verified that \(P_{\rm draw}\) is not strongly sensitive to these choices: replacing the [48] period, mass-ratio, and eccentricity prescriptions with the simplified priors common in rapid population synthesis (log-uniform in period, uniform in mass ratio, and uniform or thermal in eccentricity) yields quantitatively similar values.

3.3 Extension to hierarchical triples↩︎

For triples, we extend the parameter space to include the inner orbit and the mutual inclination, \((M_1, \log P_{\rm in}, q_{\rm in}, e_{\rm in};\;\log P_{\rm out}, q_{\rm out}, e_{\rm out};\;i)\), where "in" denotes the inner (primary-tertiary) orbit and "out" the outer orbit of the secondary about the inner pair. The mutual inclination angle \(i\) is sampled uniformly in \(\cos i \in [-1,1]\), the maximum-entropy choice corresponding to isotropic relative orientations. \[P_{\rm draw}^{\rm triple} \approx P_{\rm out} \times P_{\rm in} \times P_{\rm stable}, \label{eq:triplefactor}\tag{6}\] where \(P_{\rm out}\) is the probability that a massive primary has a companion satisfying the outer-orbit cuts - identical to the binary draw probability of equation 4 . \(P_{\rm in}\) is the conditional probability that such a system additionally hosts a tertiary satisfying the inner-orbit cuts. The last factor \(P_{\rm stable}\) is the fraction of the resulting configurations that are dynamically stable.

For the stability of the hierarchy, we adopt the criterion of [52] with the inclination-dependent factor of [53], imposed as a hard cut on the ZAMS parameter space: \[\begin{align} \begin{aligned} \left(\frac{a_{\rm out}}{a_{\rm in}}\right)_{\rm crit} & = \frac{2.8}{1 - e_{\rm out}}\left(1 - \frac{0.3\,i}{\pi}\right) \\ & \times \left[\frac{(1 + q_{\rm out})(1 + e_{\rm out})}{(1 - e_{\rm out})^{1/2}}\right]^{2/5}, \label{eq:stability} \end{aligned} \end{align}\tag{7}\] a configuration being retained only if its semi-major-axis ratio exceeds this critical value. Here \(q_{\rm out} \equiv M_2/(M_1 + M_3)\) is the ratio of the outer (secondary) mass to the total mass of the primary-tertiary inner binary.

3.4 Dynamical realizability with population synthesis↩︎

Equation 7 is evaluated at ZAMS and is therefore only a snapshot: it does not guarantee that a formally stable triple remains stable as the stars evolve and the orbits change, nor that it actually proceeds through the sequence of section 2. To verify the dynamical realizability of the channel, and in particular to establish the ZAMS period ranges over which it proceeds, we evolve the sampled triples with the Multiple Stellar Evolution (MSE) population-synthesis code ([54]). MSE switches dynamically between an orbit-averaged (secular) treatment and direct \(N\)-body integration, so that the hierarchical-stability condition is re-evaluated continuously rather than only at ZAMS, and couples this dynamics to single-star evolution from the SSE tracks of [55] and to mass-transfer and CEE prescriptions adapted from BSE ([56]).

We use these runs as a realizability check rather than as a rate prediction. The essential reason is that MSE assigns compact remnants from the SSE fitting formulae on the basis of mass alone and does not follow the rotational state of the core; it therefore cannot represent the step that defines the present channel – the spin-up of the primary’s core by the inspiralling tertiary that the JJEM requires for BH formation (section 2) - and forms BHs whether or not the core was spun up, so that its BH yields are permissive relative to the full scenario. Secondary prescription choices bias the simulated yields in both directions: the mass-only remnant criterion is permissive, whereas the CEE energy budget credits only the orbital-energy term and omits additional sources (for example, energy injected by the accreting BH during the second CEE), which would ease envelope ejection for a wider range of systems. The simulations are accordingly best read as confirming that the dynamical sequence is realizable, and as localizing the inner-orbit period range over which the channel proceeds (section 4), not as predicting its absolute rate.

4 The fraction of black hole binaries↩︎

4.1 Binary calculation↩︎

Using the ZAMS prescription described in 3.2, we first evaluate the analytic draw probability for systems capable of forming merging BBHs, treating the primary and the massive companion as an isolated interacting binary. The integration region is defined by imposing the following constraints on the ZAMS parameters.

Primary mass. We require \(M_1 \geq 25 M_\odot\). BH formation from core collapse generally requires a progenitor with ZAMS mass \(M_{\rm ZAMS} \geq 20 M_\odot\) (e.g., [2], [57]). The adopted threshold of \(25 M_\odot\) lies conservatively above this minimum, ensuring that the progenitor possesses a sufficiently large core binding energy for BH formation across a range of remnant-mass prescriptions.

Orbital period. The orbital period is restricted to \(\log(P/\text{day}) \in [2.3, 3.5]\); The upper bound lies well within the interaction regime identified by [48], who find that massive binaries with \(\log(P/\text{day}) \lesssim 3.7\) undergo Roche-lobe overflow during post-main-sequence evolution. This ensures that the binary undergoes mass transfer and/or CEE as required by the evolutionary scenario described in section 2. The lower bound of \(\log(P/\text{day}) = 2.3\) is adopted to coincide with the outer-orbit period range used in the triple calculation of section 4.2, which is itself set by the hierarchical-stability requirement of [53]; see Figure 5. Matching the binary period range to the triple outer-orbit range makes the binary and triple draw probabilities directly comparable; in the conditional-probability decomposition adopted in this paper, the binary calculation can then be read as the outer-orbit sub-case of the triple calculation. We note that this lower bound is significantly more conservative than the threshold for main-sequence interaction in massive binaries (e.g., [9], [58], [59]), which becomes important only at \(\log(P/\text{day}) \lesssim 1\); the systems excluded by our cut at \(1 \lesssim \log(P/\text{day}) < 2.3\) are therefore physically viable for post-MS interaction but are removed for the comparability reason stated above. As a cross-check, we examined the ZAMS properties of all systems that produce BBH mergers in the compas (Compact Object Mergers: Population Astrophysics and Statistics; [60]–[62] population synthesis dataset, version 02.31.06 [63]¹), shown in Figure 6. The orbital period distribution of merging systems is concentrated at \(\log(P/\text{day}) \sim 1.5-3\), with the dominant peak near \(\log(P/\text{day}) \sim 2-2.5\) and very few systems below \(\log(P/\text{day}) = 1\). The adopted range [2.3,3.5] thus lies on the high-period flank of the merging-system distribution and is conservative relative to its peak. We note that compas employs a BH formation prescription based on pre-collapse core mass and remnant-mass mappings, in which a sufficiently massive star can form a BH without the rapid core rotation required by the JJEM. The compas - derived distributions are therefore not directly applicable to the present scenario, and are used here only to verify that the physically motivated ZAMS cuts are not implausibly narrow or wide.

Secondary mass. We require \(M_2 \geq 18 M_\odot\) at ZAMS. BH formation in standard remnant-mass prescriptions requires a ZAMS-equivalent mass of approximately \(20 M_\odot\) or more (e.g., [2], [57]). A ZAMS threshold below this value is adopted because the evolutionary scenario 2 includes mass transfer from the primary RSG to the secondary prior to the first core collapse. As noted by [26], this accretion can substantially increase the secondary mass above its birth value, potentially even exceeding \(M_{1,\rm ZAMS}\). The relevance of this channel rests on the well-established observation that massive stars rarely evolve in isolation: [47] find that more than \(70 \%\) of O-type stars are members of close interacting binaries that exchange mass with a companion during their evolution. Detailed population studies further show that a substantial fraction of these massive binaries produce accretors whose mass is significantly increased through stable Roche-lobe overflow from the primary, with the gain scaling with the donor’s envelope mass and the assumed accretion efficiency ([64]). On this basis, accretion of at least \(2 M_\odot\) onto a \(18 M_\odot\) secondary - sufficient to bring its pre-collapse mass above the BH-formation threshold - is a conservative expectation. As illustrated in the left and middle panels of Figure 6, approximately 26% of the BBH-merging systems within the compas population feature \(M_{1,\mathrm{ZAMS}} < 18 M_\odot\), while 40% have secondary masses below this same threshold (\(M_{2,\mathrm{ZAMS}} < 18 M_\odot\)). This demonstrates that our adopted thresholds of \(25 M_\odot\) for the primary and \(18 M_\odot\) for the secondary are conservative: a substantial fraction of compas mergers begin below them, so relaxing the thresholds would only increase \(P_{\mathrm{draw}}\) (and lower the required \(\beta\)). The adopted threshold thus serves as a ZAMS proxy for the physical requirement on the pre-collapse mass.

Eccentricity. The eccentricity is drawn from the power-law distribution described in section 3.1 ([48]), without additional range restriction beyond the Roche-lobe-filling criterion in the periastron. Because this is the only constraint imposed on \(e\), the eccentricity is sampled freely over the physically allowed interval \([0, e_{\max}(P)]\) for each configuration \((M_1, \log P, \;q)\), and the integration of the joint probability density function over this interval contributes only a normalization factor. Eccentricity, therefore, does not impose a conditional restriction on \(P_{\mathrm{draw}}\) beyond what is absorbed in this normalization.

Integrating the empirically calibrated joint distribution (equation 5 ) over this region yields a fiducial binary draw probability \(P_{\mathrm{draw}} \simeq 3 \times 10^{-2}\). The constituent conditional probabilities that contribute to this aggregate value are detailed below. For consistency, all subsequent mass and temporal quantities are reported in solar masses (\(M_\odot\)) and days, respectively.

\[\tag{8} \begin{align} & P\!\left(M_1 \in [25,150] \mid M_1 \ge 8\right) && \simeq 0.199, \tag{9}\\[7pt] & P\!\left(n > 0 \mid M_1 \in [25,150]\right) && \simeq 0.915, \tag{10}\\[7pt] & \begin{aligned} &P\bigl(\log P \in [2.3,3.5] \mid n>0, \\[1pt] &\quad M_1 \in [25,150]\bigr) \end{aligned} && \simeq 0.375, \tag{11}\\[7pt] & \begin{align} &P\bigl(M_2 \in [18,M_1] \mid \log P \in [2.3,3.5], \\[1pt] &\quad n>0,\, M_1 \in [25,150]\bigr) \end{align} && \simeq 0.441, \tag{12} \end{align}\]

When inserted into the CCSN-normalized rate relation (equation 3 ), this corresponds to a success fraction \(\beta_{\mathrm{binary}} \simeq 0.023\). This value is physically plausible: it is consistent with the expectation that only a small fraction of initially favorable systems survive the combined uncertainties of CE ejection efficiency, natal kick disruption, and the requirement of merger within a Hubble time.

A further effect that reduces the number of merging-BBH systems is the requirement of a tertiary star (Section 2), which we examine next.

4.2 Triple calculation↩︎

We extend the calculation to hierarchical triples by introducing an additional inner companion to the primary - the tertiary star of mass \(M_{3}\) (Figures 2 and 3).

Outer orbit. The outer orbit is retained within the same intervals \(\log(P_{\text{in}}/\text{day}) \in [2.3, 3.5]\) as in the binary case, and the outer companion mass is required to satisfy \(M_{2} \in [18,150] M_\odot\), (with the implicit constraint \(M_{2} < M_{1}\)), for the same physical reasons discussed in section 4.1. The outer-orbit parameters therefore contribute the same conditional probabilities as in the binary calculation.

Inner orbit. The inner orbital period is restricted to \(\log(P_{\text{in}}/\text{day}) \in [0.7,1.5]\). This choice is motivated by several considerations. This interval lies within the empirically established interaction regime for massive binaries (e.g., [48]): the tertiary must be engulfed once the primary expands after the MS, which requires a sufficiently tight inner orbit. Both bounds are motivated by the time evolution runs with the Multiple Stellar Evolution (MSE) code ([54]) and have distinct physical origins (Figure 7). The lower bound is set by the onset of the inner CEE: at \(\log(P_{\text{in}}/\text{day}) < 0.7\) the inner orbit is so tight that the primary fills its Roche lobe, or otherwise interacts with the tertiary, while both stars are still on the MS - that is, before the primary expands after the MS and drives the inner CEE of the channel. In the pilot runs, no system below \(\log(P_{\text{in}}/\text{day}) = 0.7\) initiates the inner CEE at all (Figure 7a). The upper bound is set, by contrast, by the failure of the second CEE. The ZAMS stability criterion alone (equation 7 ) would permit inner periods up to \(\log(P_{\text{in}}/\text{day}) \simeq 1.85\), but the pilot runs show the channel ceasing to complete well below this value: the fraction of systems that proceed to the second CEE between the first-born BH and the secondary declines with increasing inner period and drops to zero above \(\log(P_{\text{in}}/\text{day}) \simeq 1.4-1.5\) (Figure 7b). Such wide-inner-orbit triples can still undergo the first CEE but do not complete the remainder of the sequence, in part because a wider inner orbit leaves the system more susceptible to dynamical instability as the stars evolve. The value \(\log(P_{\text{in}}/\text{day}) = 1.5\) is therefore adopted as the period above which the channel is no longer realized in the simulations, rather than as the formal ZAMS-stability limit. Both bounds thus reflect the evolution of the triple and not merely its zero-age configuration: the ZAMS stability criterion is necessary but not sufficient, and the time-evolution runs localize the inner-period range over which the full channel proceeds.

Inner companion mass. The inner companion (tertiary) mass is required to satisfy \(M_3 \in [1.25,4.5] M_\odot\). The upper bound reflects the requirement \(M_3 < 5 M_\odot\) from the evolutionary scenario (Section 2): the tertiary must be a low-mass star that can be tidally disrupted by the primary’s core during the first CEE, rather than a massive star that would itself significantly affect the system’s evolution prior to being engulfed. The lower bound is not an independent mass cut but a consequence of the sampling: it follows from the mass-ratio floor \(q_{\min}=0.05\) the extrapolation limit of the [48] distribution adopted in section 3.2 - together with the minimum primary mass \(M_1 = 25 M_\odot\), so that \(M_{\mathrm{in}} \geq q_{\min} \times M_{1,\min} = 1.25 M_\odot\). The corresponding inner mass ratios \(q_{\mathrm{in}} = M_{\mathrm{in}}/M_{1}\) thus span approximately \(0.05-0.3\) for primaries \(M_1 \geq 25 M_\odot\), lying at the low end of - and partially below - the calibrated domain of [48] (\(q \geq0.1\)). Tertiaries near this lower bound carry only a small orbital angular momentum relative to the primary; the spin-up of the primary’s core required by the JJEM (Section 2) is therefore expected to be most effective for tertiaries toward the upper part of the adopted mass range.

Dynamical stability. The [53] stability criterion (equation 7 ) is imposed as a hard cut in the ZAMS parameter space. The fraction of Monte Carlo draws that satisfy the stability criterion, among systems that already meet all other constraints, is \(P_{\rm stable} \simeq 0.58\).

Within this fiducial triple region, the conditional probabilities contributing to the draw probability are: (the following mass and temporal quantities are in units of solar masses (\(M_\odot\)) and days, respectively).

\[\tag{13} \begin{alignat}{2} & P\!\left(M_1 \in [25,150] \mid M_1 \ge 8\right) && \simeq 0.199, \tag{14}\\[7pt] & P\!\left(n_{\mathrm{in}} > 0 \mid M_1 \in [25,150]\right) && \simeq 0.961, \tag{15}\\[7pt] & P\!\left(n_{\mathrm{out}} > 0 \mid M_1 \in [25,150]\right) && \simeq 0.915, \tag{16}\\[7pt] & \begin{align} &P\bigl(\log P_{\mathrm{in}} \in [0.7,1.5] \mid n_{\mathrm{in}}>0, \\[1pt] &\quad M_1 \in [25,150]\bigr) \end{align} && \simeq 0.283, \tag{17}\\[7pt] & \begin{aligned} &P\bigl(\log P_{\mathrm{out}} \in [2.3,3.5] \mid n_{\mathrm{out}}>0, \\[1pt] &\quad M_1 \in [25,150]\bigr) \end{aligned} && \simeq 0.375, \tag{18}\\[7pt] & \begin{align} &P\bigl(M_3 \in [1.25,4.5] \mid \log P_{\mathrm{in}} \in [0.7,1.5], \\[1pt] &\quad n_{\mathrm{in}}>0,\, M_1 \in [25,150]\bigr) \end{align} && \simeq 0.221, \tag{19}\\[7pt] & \begin{align} &P\bigl(M_2 \in [18,M_1] \mid \log P_{\mathrm{out}} \in [2.3,3.5], \\[1pt] &\quad n_{\mathrm{out}}>0,\, M_1 \in [25,150]\bigr) \end{align} && \simeq 0.441, \tag{20} \end{alignat}\]

Expressed as a fraction of the massive-primary population, the calculation implies that roughly \(0.9 \%\) of primaries with \(M_1 \geq 25 M_\odot\) host companions at both the inner and outer separations of interest, of which roughly half form dynamically stable configurations - a modest but non-negligible fraction.

Inserting \(P_{\rm draw}^{\rm triple}\) into the CCSN-normalized rate relation (equation 3 ) gives the success fraction for the triple channel \(\beta_{\mathrm{triple}} \lesssim 0.67\). As in the binary case, this value is best read as an upper bound: it is the per-system success fraction the channel would require if it produced the entire observed BBH merger rate, whereas in practice it contributes alongside other formation channels. This bound should be interpreted cautiously, because \(\beta_{\mathrm{triple}}\) is high - about thirty times the corresponding binary value (\(\beta_{\mathrm{binary}} \simeq 0.023\)), the difference reflecting the multiplicity penalty that makes suitable triples rarer. Producing the entire observed rate through this channel would require roughly four in five of the suitable triples to complete the full sequence of section 2 - the inner CEE and tidal disruption of the tertiary, two core collapses each forming a BH, a second CEE that leaves a bound and sufficiently tight orbit, and coalescence within a Hubble time - despite the large uncertainties that \(\beta\) absorbs, including CE ejection efficiency, natal kicks, mass loss, and the JJEM core-rotation requirement that the simulations do not test (section 3.4). A success fraction this large is demanding for so long and uncertain an evolutionary chain. We therefore do not read \(\beta_{\mathrm{triple}} \simeq 0.67\) as evidence that the triple channel produces the whole BBH population; the more natural reading is that it is a viable contributor supplying some fraction of the observed rate, in which case the required per-system success fraction falls proportionally to a more comfortable value. The bound nonetheless lies below unity, so the channel is not excluded on rate grounds and is promising as a contributor to the merging-BBH population - but the present work does not establish that it dominates that population.

4.3 Expanding the boundaries of the constraints↩︎

The fiducial primary-mass cut, \(M_1 \geq 25\,M_\odot\), was placed conservatively above the \(\simeq 20\,M_\odot\) minimum generally required for BH formation (section 4.1). To gauge how sensitive the inferred success fraction is to this choice, we repeat the triple calculation with the primary mass relaxed to the BH-formation floor, \(M_1 \in [20,150]\,M_\odot\). Admitting lower-mass primaries enlarges the qualifying population and thus raises the draw probability; since \(\beta \propto P_{\rm draw}^{-1}\) (equation 3 ), this can only reduce the required success fraction. The conditional probabilities for the relaxed cut are: (the following mass and temporal quantities are in units of solar masses (\(M_\odot\)) and days, respectively).

\[\tag{21} \begin{alignat}{2} & P\!\left(M_1 \in [20,150] \mid M_1 \ge 8\right) && \simeq 0.276, \tag{22}\\[7pt] & P\!\left(n_{\mathrm{in}} > 0 \mid M_1 \in [20,150]\right) && \simeq 0.953, \tag{23}\\[7pt] & P\!\left(n_{\mathrm{out}} > 0 \mid M_1 \in [20,150]\right) && \simeq 0.9, \tag{24}\\[7pt] & \begin{align} &P\bigl(\log P_{\mathrm{in}} \in [0.7,1.5] \mid n_{\mathrm{in}}>0, \\[1pt] &\quad M_1 \in [20,150]\bigr) \end{align} && \simeq 0.267, \tag{25}\\[7pt] & \begin{aligned} &P\bigl(\log P_{\mathrm{out}} \in [2.3,3.5] \mid n_{\mathrm{out}}>0, \\[1pt] &\quad M_1 \in [20,150]\bigr) \end{aligned} && \simeq 0.367, \tag{26}\\[7pt] & \begin{align} &P\bigl(M_3 \in [1.25,4.5] \mid \log P_{\mathrm{in}} \in [0.7,1.5], \\[1pt] &\quad n_{\mathrm{in}}>0,\, M_1 \in [20,150]\bigr) \end{align} && \simeq 0.307, \tag{27}\\[7pt] & \begin{align} &P\bigl(M_2 \in [18,M_1] \mid \log P_{\mathrm{out}} \in [2.3,3.5], \\[1pt] &\quad n_{\mathrm{out}}>0,\, M_1 \in [20,150]\bigr) \end{align} && \simeq 0.355, \tag{28} \end{alignat}\]

These factors give an outer draw probability \(P_{\rm out} \simeq 3.3\times10^{-2}\) and an inner conditional probability \(P_{\rm in} \simeq 7.8\times10^{-2}\), while the fraction of stability-satisfying configurations is \(P_{\rm stable} \simeq 0.57\). Expressed as a fraction of the massive-primary population, roughly \(0.9\%\) of primaries with \(M_1 \geq 20\,M_\odot\) host companions at both the inner and outer separations of interest, of which slightly more than half remain dynamically stable. The total draw probability is therefore \(P_{\rm draw}^{\rm triple} \simeq 1.5\times10^{-3}\) (equation 6 ), which, inserted into the CCSN-normalized relation (equation 3 ), gives \(\beta_{\rm triple} \simeq 0.47\).

Relaxing the primary-mass cut to the BH-formation floor thus lowers the required success fraction from the fiducial \(\beta_{\rm triple} \simeq 0.67\) (section 4.2) to \(\simeq 0.47\). The channel’s viability does not hinge on the precise primary-mass threshold: a wider primary range only moves the inferred \(\beta\) further below unity.

5 Summary↩︎

GW observations show that the distribution of the effective inspiral spin parameter \(\chi_{\mathrm{eff}}\) (equation 1 ) of merging BBHs is asymmetric about zero: it peaks slightly positive but retains a significant negative tail (Figure 1), with an inferred negative fraction of \(0.3-0.46\) ([28]). Neither isolated binary evolution, which predicts \(\chi_{\mathrm{eff}}\) narrowly peaked near zero, nor gas-free dynamical assembly, which predicts a distribution symmetric about zero, reproduces this asymmetry on its own (Section 1).

In this paper, we examined whether a hierarchical triple-star channel, which we formulated within the framework of the JJEM, can contribute a non-negligible fraction of the misaligned population. In this channel (Figures 2 and 3), the primary engulfs a low-mass tertiary in a first CEE, whose inspiral spins up the primary’s core and enables BH formation according to the JJEM. The secondary subsequently engulfs the first-born BH in a second CEE and collapses to the second-born BH (Section 2). The spin–orbit misalignment arises from the inclination between the inner and outer orbital planes, rather than from a natal kick or from purely dynamical forcing. Namely, the misalignment is inherited from the formation of the triple system.

We set constraints on the parameters of the triple system (equation 13 ): the masses of the three stars, the inner and outer orbital periods, and the eccentricity. We demanded that during the main-sequence phase of the primary star, the triple system be dynamically stable. We semi-analytically estimated the expected fraction of merging BBHs under our proposed scenario, and determined its ratio to the CCSN explosion rate (equation 2 ). Comparing this ratio to the observed one, we found the fraction \(\beta\) of systems obeying our initial criteria that should end as observed merging BBH. For our fiducial triple-star parameters, we found that a fraction of \(\beta_{\mathrm{triple}} \simeq 0.67\) of our initial systems should end as merging BBH. As this value lies below unity, the evolutionary channel we propose might contribute significantly to merging BBHs. Namely, the rate does not exclude this channel from forming most or all of the merging BBHs. Relaxing somewhat the constraints, like allowing a lower minimum mass for the primary star (equation 21 ), eases the formation of systems, and reduces the value of \(\beta_{\mathrm{triple}}\).

The main advantage of our proposed channel is that it produces both a positive average \(\chi_{\mathrm{eff}}\) and a significant negative tail (Section 2). Therefore, it can qualitatively reproduce the main features of the observed distribution and is therefore a promising non-negligible contributor to the misaligned population of merging BBHs. Whether it dominates this population or contributes less than half of it is a subject for a future study. As well, there is a need to derive the quantitative distribution of \(\chi_{\mathrm{eff}}\).

Acknowledgements↩︎

NS thanks the Charles Wolfson Academic Chair at the Technion for the support.

6 Conversion of the BBH merger rate to z=0↩︎

[28] reports the density of the BBH merger-rate at \(z=0.2\) rather than at \(z=0\). Because our normalization compares this quantity with the local (\(z \simeq 0\)) CCSNe rate ([49]), the two must be evaluated in the same redshift; therefore, we convert the GWTC-5.0 value to \(z=0\).

GWTC-5.0 models the redshift distribution of mergers as \[\pi(z\mid\kappa_z) \;\propto\; \frac{1}{1+z}\,\frac{dV_c}{dz}\,(1+z)^{\kappa_z}, \label{eq:piz}\tag{29}\] where \(dV_c/dz\) is the comoving volume element and the factor \(1/(1+z)\) accounts for the time dilation between the source and detector frames. The astrophysical content of this model is carried out by the comoving merger-rate density, which scales as \(\mathcal{R}(z) \propto (1+z)^{\kappa_z}\); this is the quantity tabulated by GWTC-5.0 and the one relevant to our normalization. The remaining factors in equation 29 convert this rate density into the observed distribution of detections and do not enter the rescaling of the rate density between redshifts. Because \(\mathcal{R}(z)\propto(1+z)^{\kappa_z}\), the ratio of the rate at two redshifts is independent of its overall normalization, \[\frac{\mathcal{R}(z=0)}{\mathcal{R}(z=0.2)} = \left(\frac{1+0}{1+0.2}\right)^{\kappa_z} = (1.2)^{-\kappa_z},\] so that \[\mathcal{R}(z=0) = \mathcal{R}(z=0.2)\,(1.2)^{-\kappa_z}, \label{eq:conv}\tag{30}\]

GWTC-5.0 infers \(\kappa_z = 2.5\pm0.7\) for the BBH population, consistent with the redshift evolution of the star-formation rate ([65]). With \(\kappa_z=2.5\) the conversion factor is \((1.2)^{-2.5} \simeq 0.633\) and equation 30 maps the reported interval \(27.5-49.4 \;\mathrm{Gpc^{-3}\,yr^{-1}}\) to \[\mathcal{R}_{\rm BBH}(z=0) \simeq 17.4-31.3 \;\mathrm{Gpc^{-3}\,yr^{-1}}.\] The uncertainty on \(\kappa_z\) propagates only weakly: over \(\kappa_z\in[1.8,3.2]\) the factor ranges from \((1.2)^{-1.8} \simeq 0.72\) to \((1.2)^{-3.2} \simeq 0.55\), a \(\sim 15\%\) spread that is subdominant to the width of the rate interval itself. As an independent check, the converted value is consistent with the \(z=0\) BBH rate reported directly by the earlier GWTC-4.0 analysis, \(\mathcal{R}_{\mathrm{BBH}} = 14-26 \;\mathrm{Gpc^{-3}{~\rm yr}^{-1}}\) ([1]).

References↩︎