© The Authors 2025

1. Introduction

The detection of gravitational waves (GWs) has not only confirmed their existence, as predicted by general relativity (Einstein 1916, 1918), but has also revolutionized our understanding of compact object dynamics and opened an entirely new observational window to the Universe. With over 90 confirmed detections to date (Abbott et al. 2023), GW astronomy has rapidly evolved from a theoretical construct to a data-rich observational science, profoundly impacting multiple areas of astrophysics from stellar evolution to cosmology (e.g., Miller 2016; Mapelli 2018; Mapelli et al. 2021; Mandel & Farmer 2022; Mandel & Broekgaarden 2022). However, these remarkable observations have brought into sharp focus a fundamental challenge in our understanding of compact object binary formation and evolution, commonly referred to as the “initial separation” problem or “final AU” problem (Mandel & Farmer 2022). For two black holes (BHs) to merge within the age of the Universe through GW emission alone, they must begin at an extraordinarily small separation–typically less than 50 R_⊙ for 30 M_⊙ BHs. This requirement creates a puzzling apparent contradiction: the massive stars that are the progenitors of these BHs have radii that are much larger than this critical separation during significant portions of their evolution. If such massive stars were placed at separations small enough to enable their eventual BH remnants to merge within a Hubble time, the stars themselves would merge long before they could evolve into BHs.

This contradiction has spurred intense theoretical investigation into various formation pathways that might resolve this tension. The proposed mechanisms fall broadly into two categories: stellar evolution solutions and dynamical solutions. Stellar evolution approaches include mass transfer channels. In stable mass transfer (Inayoshi et al. 2017; Gallegos-Garcia et al. 2021; van Son et al. 2022; Briel et al. 2023), one star fills its Roche lobe and undergoes controlled mass exchange with its companion, potentially leading to orbital tightening. In common-envelope evolution (Paczynski 1976; Livio & Soker 1988; Ivanova et al. 2013; Postnov & Yungelson 2014), dynamically unstable mass transfer triggers the engulfment of the mass-gaining star within the donor’s envelope, where the orbital energy is used to eject the common envelope and dramatically shrink the binary separation. Another stellar evolution pathway is chemically homogeneous evolution (Marchant et al. 2016; Mandel & de Mink 2016; de Mink & Mandel 2016; du Buisson et al. 2020; Riley et al. 2021), where strong internal mixing prevents stellar expansion, allowing initially close binaries to remain compact throughout their evolution; however, this pathway has become increasingly disfavored due to the lack of high-spin detections in the current GW catalog. Dynamical solutions, on the other hand, focus on mechanisms that can bring already-formed BHs closer together, primarily through interactions in dense stellar environments such as globular clusters (Sigurdsson & Hernquist 1993; Kulkarni et al. 1993; Portegies Zwart & McMillan 2000; Banerjee et al. 2010; Morscher et al. 2015; Rodriguez et al. 2016; Askar et al. 2017; Fragione & Kocsis 2018; Mapelli et al. 2021) or through secular interactions in hierarchical multiple systems (Antonini et al. 2016; Silsbee & Tremaine 2017; Hoang et al. 2018; Martinez et al. 2020; Dittmann et al. 2024). The hierarchical triple scenario, which is the focus of this paper, leverages the fact that a significant fraction of massive stars form in triple and higher-order multiple systems (Duchêne & Kraus 2013; Tokovinin 2021). In such configurations, the von Zeipel–Lidov–Kozai (ZLK) mechanism (von Zeipel 1910; Lidov 1962; Kozai 1962) can induce oscillations in the inner binary’s eccentricity. During periods of high eccentricity, the pericenter distance becomes extremely small, dramatically increasing GW emission and potentially leading to a merger within the age of the Universe. An important manifestation of this phenomenon is shown for instance in the recent work of Stegmann & Klencki (2025), suggesting that the present-day rate of BH–neutron star mergers may be largely driven by ZLK resonance.

While numerous aspects of the hierarchical triple pathway have been studied in recent years (Liu & Lai 2020, 2021; Su et al. 2021; Britt et al. 2021; Mannerkoski et al. 2021; Trani et al. 2022; Chandramouli & Yunes 2022; Gondán 2023; Lepp et al. 2023; Bartos et al. 2023; Liu et al. 2024; Generozov & Perets 2024; Hao et al. 2024), the literature lacks a comprehensive exploration of the full parameter space relevant to BH mergers. Previous works have typically focused on specific subsets of the parameter space, particular mass ranges, or limited samples of systems, making it difficult to develop a holistic understanding of which initial configurations are most likely to produce mergers within the age of the Universe. Our primary aim in this work is to conduct a systematic, large-scale investigation of the parameter space of hierarchical triple BH systems, with a specific focus on identifying the boundary between initial configurations that lead to mergers within 14 Gyr and those that do not. By mapping this boundary in detail, we can provide valuable insights into the viability of the hierarchical triple channel as a solution to the initial separation problem and guide future population synthesis studies seeking to predict the rates and properties of GW sources.

To make this extensive parameter space exploration computationally feasible, we employed a secular approximation that averages the equations of motion over both orbital periods, drastically reducing the computational cost compared to directN-body integration. We validated this approach by comparing a subset of our results with fully self-consistent N-body simulations. Furthermore, we developed an innovative Markov chain Monte Carlo (MCMC) sampling strategy that preferentially explores the transition region between merging and nonmerging configurations, allowing us to efficiently map the boundary in our high-dimensional parameter space. As an additional contribution, we trained a neural network on the results of our extensive simulations to create a predictive model that can rapidly classify the merger outcome of any hierarchical triple BH system within our parameter space. The aim of this model is to achieve a trustworthy accuracy, while requiring only milliseconds of computation time per system, enabling future population synthesis studies to efficiently investigate the hierarchical triple channel without the computational burden of full dynamical simulations.

This paper is organized as follows: Section 2 provides a contextual description of the initial separation problem and the ZLK mechanism. Section 3 outlines our methodology, including the secular approximation, MCMC sampling technique, and neural network implementation. Section 4 presents our findings. We compare secular and N-body simulations, characterize the parameter space of merging and nonmerging systems, and evaluate the performance of the neural network predictor. Finally, in Sect. 5 we discuss the limitations of our approach, summarize our conclusions, and explore the implications for future research.

2. Context

The merger timescale τ_GW of a binary system due to GW emission can be derived from general relativity. For a circular binary with component masses m and M, the mergertimescale is

$$\begin{array}{c}\hfill {\tau }_{\mathrm{GW}}=\frac{5}{256}\frac{c^5}{{\mathcal{G}}^3}\frac{a^4}{(mM)(m+M)},\end{array}$$(1)

where a is the binary separation, 𝒢 the gravitational constant, and c the speed of light. For two 10 M_⊙ BHs, the maximum initial separation for which τ_GW ≤ 14 Gyr (approximately the age of the Universe) is about 12 R_⊙ (Eq. (1)). This creates an apparent contradiction: if the progenitor stars were located this close to each other, they would merge before becoming BHs, due to their much larger radii during evolution. Conversely, if located at larger separations to avoid stellar mergers, the resulting BH binary would not merge within the age of the Universe through GW emission alone.

The ZLK resonance provides a potential solution to this problem. In a hierarchical triple configuration, where an inner binary is orbited by a distant third body, this mechanism can induce large-amplitude oscillations in the inner binary’s eccentricity, while conserving angular momentum. These oscillations can temporarily drive the inner binary’s eccentricity to extremely high values, dramatically reducing the pericenter distance during certain phases of the cycle. Since GW emission is strongly enhanced at small separations, these periodic excursions to high eccentricity can accelerate the orbital decay by orders of magnitude, potentially enabling mergers within cosmological timescales even for systems that begin with large separations. The effectiveness of the hierarchical triple pathway depends on a complex interplay of multiple factors. The ZLK mechanism requires a sufficiently large mutual inclination i_mut between the inner and outer orbits, typically between 39.2° and 140.8° for circular outer orbits (Kinoshita & Nakai 1999), though this range can expand for eccentric outer orbits. The strength of the perturbation, and hence the maximum eccentricity attainable by the inner binary, depends on the mass ratios, orbital separations, and eccentricities (Naoz 2016, and references therein). This phenomenon must also be tempered by the influence of relativistic precession of the inner orbit, which acts as a competing mechanism that can quench the ZLK oscillations if sufficiently strong (e.g., Liu et al. 2015).

3. Methods

3.1. Secular evolution with the JADE code

We use a modified version of the JADE¹ code (Attia et al. 2021), which was originally developed to simulate the secular evolution of hierarchical triple systems in the context of planetary dynamics. The code treats the three-body problem using a secular approximation, where the differential equations of motion are averaged over both orbital motions. This allows much longer timesteps and faster computations compared to direct N-body integration.

The Hamiltonian of a hierarchical triple system can be expressed as (Harrington 1968)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathcal{H}& =\frac{\mathcal{G}{M}_1{M}_2}{2a_{\mathrm{inner}}}+\frac{\mathcal{G}{M}_3(M_1+M_2)}{2a_{\mathrm{outer}}}\hfill \\ \hfill +\frac{\mathcal{G}}{a_{\mathrm{outer}}}& \sum _{n=2}^{\infty }{\left(\frac{a_{\mathrm{inner}}}{a_{\mathrm{outer}}}\right)}^n{M}_n{\left(\frac{|{\mathit{r}}_{\mathbf{inner}}|}{a_{\mathrm{inner}}}\right)}^n{\left(\frac{a_{\mathrm{outer}}}{|{\mathit{r}}_{\mathbf{outer}}|}\right)}^{n+1}{P}_n(cos\mathrm{\Phi }),\hfill \end{array}\end{array}$$(2)

where M₁ and M₂ are the masses of the inner binary components, M₃ is the mass of the outer companion, a_inner and a_outer are the semimajor axes of the inner and outer orbits, r_inner and r_outer are the Jacobi coordinates (e.g., Murray & Dermott 1999), Φ is the angle between these vectors, P_n are the Legendre polynomials, and

$$\begin{array}{c}\hfill M_n=M_1{M}_2{M}_3\frac{M_1^{n-1}-{(-M_2)}^{n-1}}{{(M_1+M_2)}^n}.\end{array}$$(3)

The orbital configuration is depicted in Fig. 1. Following Beust et al. (2012), JADE truncates the infinite series of Eq. (2) at the fourth (hexadecapolar) order, which provides a good balance between accuracy and computational efficiency for hierarchical systems. To the best of our knowledge, this is the highest order that has been investigated for such secular analyses.

Fig. 1. Generic orbital configuration of a hierarchical triple system, as investigated in this study.

In addition to gravitational interactions (Eq. (2)), JADE also accounts for relativistic precession of the inner orbit, which acts as a quenching mechanism for the ZLK resonance (e.g., Liu et al. 2015). Furthermore, since the original JADE code was developed for planetary systems, we modified it to include the effects of GW emission on the inner binary. The loss of orbital angular momentum and eccentricity damping due to GWs are given by (Peters 1964)

$$\begin{array}{c}\hfill {\frac{\mathrm{d}{L}_{\mathrm{inner}}}{\mathrm{d}t}|}_{\mathrm{GW}}=-\frac{32}{5}\frac{{\mathcal{G}}^{7/2}}{c^5}\frac{{(M_1{M}_2)}^2\sqrt{M_1+M_2}}{a_{\mathrm{inner}}^{7/2}}\frac{1+\frac{7}{8}{e}_{\mathrm{inner}}^2}{{(1-e_{\mathrm{inner}}^2)}^2}\end{array}$$(4)

and

$$\begin{array}{c}\hfill {\frac{\mathrm{d}{e}_{\mathrm{inner}}}{\mathrm{d}t}|}_{\mathrm{GW}}=-\frac{304}{15}\frac{{\mathcal{G}}^3}{c^5}\frac{M_1{M}_2(M_1+M_2)}{a_{\mathrm{inner}}^4}\frac{1+\frac{121}{304}{e}_{\mathrm{inner}}^2}{{(1-e_{\mathrm{inner}}^2)}^{5/2}}.\end{array}$$(5)

We implemented Eqs. (4), (5) into the JADE code by including their contribution to the evolution of the Runge–Lenz and reduced orbital angular momentum vectors (we refer to Attia et al. 2021, for broader details about secularization). The outer orbit is assumed to emit a negligible amount of GWs, due to its much larger separation. Consistently with any hierarchical triple problem, the perturbing orbit hence acts as an angular momentum reservoir.

3.2. Parameter space and MCMC sampling

We explored a seven-dimensional (7D) parameter space defined by the following:

• Masses of the inner binary BHs: M₁, M₂ ∈ [5,  100] M_⊙ with M₂ ≤ M₁ (as the inner masses are interchangeable);

• Mass of the outer BH: M₃ ∈ [1,  200] M_⊙;

• Semimajor axes: a_inner ∈ [1,  200] AU, a_outer ∈ [100,  10 000] AU with a_outer > 10 a_inner (hierarchical approximation);

• Outer orbit eccentricity: e_outer ∈ [0,  0.9]. We excluded the ]0.9,  1] range because we found a large discrepancy with N-body simulations for the highest eccentricities (Sect. 4.1);

• Mutual inclination: i_mut ∈ [40° , 80° ]. We restricted our study to prograde orbits based on the well-established symmetry of ZLK dynamics about 90°, which holds when the angular momentum ratio of the inner and outer orbits is low (Anderson et al. 2017). This symmetry remains valid for ratios up to ∼0.1 (Mangipudi et al. 2022), and our hierarchical condition ensures that most of our parameter space falls well below this threshold, except near the boundaries where the inner and outer angular momenta can become commensurate. By excluding the region i_mut ∈ ]80° , 90° ], where we observed significant discrepancies with N-body simulations (Sect. 4.1), we created two disconnected parameter space regions—the prograde regime we explore here and a roughly equivalent retrograde regime (i_mut ∈ [100° , 140° ]) whose dynamics can be inferred through the symmetry relationship.

We did not vary the inner eccentricity e_inner because we know that if the ZLK mechanism is active, it will oscillate through its full range regardless of the initial value.

To efficiently sample the boundary between merging and nonmerging regions, we employed an MCMC approach. First, we simulated a grid of 300 000 systems spanning the parameter space. For each system, we computed whether it merged within 14 Gyr or not, assigning a value of 1 (merger) or 0 (no merger). We then defined a scoring function for any point x in the parameter space based on its k nearest neighbors in this grid:

$$\begin{array}{c}\hfill f(\mathit{x})=\frac{1}{N}\sum _{{\mathit{x}}_k\in {\mathcal{N}}_k(\mathit{x})}\frac{M({\mathit{x}}_k)}{d(\mathit{x},{\mathit{x}}_k)}.\end{array}$$(6)

Here 𝒩_k(x) is the set of k = 100 nearest neighbors to x, M(x) is 1 if system x merges and 0 otherwise, d(x, y) is the standardized Euclidean distance between x = (x₁,x₂,…,x₇) and y = (y₁,y₂,…,y₇) such that

$$\begin{array}{c}\hfill d(\mathit{x},\mathit{y})=\sqrt{\sum _{i=1}^7\frac{{(x_i-y_i)}^2}{V(i)}},\end{array}$$(7)

with V(i) the variance of parameter i over all computed points (to avoid the distance measure being dominated by the largest scales). Finally, N is a normalization term:

$$\begin{array}{c}\hfill N=\sum _{{\mathit{x}}_k\in {\mathcal{N}}_k(\mathit{x})}\frac{1}{d(\mathit{x},{\mathit{x}}_k)}.\end{array}$$(8)

This weighted average f(x) (Eq. (6)) ranges from 0 to 1, representing the likelihood of merger for a system with parameters x. We are most interested in regions where f(x)≃0.5, which corresponds to the boundary between merging and nonmerging configurations.

We convert f(x) to a log-likelihood function g(x) for the MCMC,

$$\begin{array}{c}\hfill g(\mathit{x})=log[1-{\left(\frac{|f(\mathit{x})-0.5|}{0.5}\right)}^2],\end{array}$$(9)

which gives the highest scores to systems with f(x)≃0.5 and strongly penalizes systems with f(x) near 0 or 1. We note that any strictly positive exponent inside the log function would make Eq. (9) a valid log-likelihood function. We chose the value of 2 because it is the lowest exponent for which g is C^∞.

We ran the MCMC with 8000 parallel workers for approximately 3000 iterations each. Importantly, our approach deviates from the classical MCMC approach in a substantial way: as new samples are simulated during the MCMC process, they are added to the grid of precomputed systems used for scoring. This means the log-likelihood function g(x) evolves throughout the process (since 𝒩_k(x) is refined, Eq. (6)) as more data become available, making our method self-improving. When a proposed sample is rejected, we still simulate its evolution and add the result to our database, thereby continuously enhancing the accuracy of our nearest-neighbor estimator. This adaptive approach is more effective than standard MCMC, as the k nearest neighbors of a new sample become progressively closer and more representative as the grid density increases. In total, we simulated nearly 15 million different initial configurations; when duplicates are removed, we have 14 456 391 unique systems across all MCMC chains, including burn-in samples. After discarding the first 50% of each chain as burn-in (to remove samples from the initial random exploration phase), we retained 14 006 724 samples in our final dataset, preserving duplicates. It is important to retain duplicates as they represent high-scoring regions that were sampled multiple times, and their frequency provides statistical weight in our analysis.

3.3. Neural network for outcome prediction

To enable the rapid prediction of merger outcomes without full dynamical simulations, we trained a readily usable open-source neural network² on the results of our MCMC exploration. After removing duplicates, we had 14 456 391 unique systems with known outcomes, which we split into training (60%), validation (20%), and test (20%) sets.

Our neural network architecture, typical of binary classification schemes, consists of

• an input layer with seven neurons (one for each parameter);

• three hidden layers with 128 neurons each using ReLU activation functions;

• an output layer with a single neuron using a sigmoid activation function.

The input parameters were normalized to the range [0,  1]. The network was trained for 15 epochs, at which point the validation accuracy no longer improved.

Although the model is trained on binary data, since the real outcome y = M(x) is either 0 or 1, its predictions $\widehat{y}$ are values between 0 and 1, which can be interpreted as the probability of merger. For binary classification, we use a threshold of 0.5, with outputs ≥0.5 classified as mergers. We define in this framework a confidence measure

$$\begin{array}{c}\hfill c(\widehat{y})=2|\widehat{y}-0.5|,\end{array}$$(10)

which ranges from 0 (the model is completely uncertain of its prediction) to 1 (completely certain).

4. Results

4.1. Comparison with N-body simulations

To validate our secular approach, we performed an extensive comparative analysis between the modified JADE code and the TSUNAMI direct N-body integrator (Trani & Spera 2023). To this end, our validation set comprised 1000 systems, providing a statistically robust assessment of the secular approximation’s reliability across different parameter regimes. We divided our validation sample into four categories:

• Predictable outcomes (300 systems). We selected 150 systems expected to merge and 150 expected not to merge based on their location in well-characterized regions of parameter space. The secular and N-body codes agreed perfectly on the outcome (merger or no merger) for all 300 systems, confirming that our approach reliably captures the dynamics in regions far from the merger boundary.

• Uncertain outcomes (500 systems). For systems sampled from regions where the posterior merger probability was uncertain (0.45 ≤ f(x)≤0.55), we found good agreement between the two approaches. Of the 500 systems tested, 18 were dynamically unstable (the triple disintegrated during evolution). Among the remaining 482 stable systems, the secular code correctly predicted 382 outcomes (177 true positives and 205 true negatives), yielding a 79% accuracy rate. The disagreements consisted of 69 false positives and 31 false negatives, with no strong bias toward either type of error.

• High mutual inclination regime (100 systems). Systems with 80° < i_mut < 90° showed poor agreement between the two methods. Of 100 randomly sampled systems in this range, 12 were dynamically unstable, and among the remaining 88, only 45 (51%) showed consistent outcomes, which motivated our decision to exclude this parameter region from our main analysis.

• High outer eccentricity regime (100 systems). Systems with e_outer > 0.9 similarly showed poor agreement. Of the 100 systems tested, 37 were unstable; among the 63 stable systems, only 17 (27%) showed consistent outcomes. The breakdown of the secular approximation at very high eccentricities likely occurs because the outer companion’s close pericenter passages violate the assumption of well-separated orbital timescales.

Excluding the two problematic parameter regions (high i_mut and high e_outer), our overall validation accuracy reaches 87% for stable systems. This high agreement rate validates our methodology for large-scale parameter space exploration, while clearly delineating the regimes where the secular approximation should not be trusted. The dynamical instabilities we observed highlight a limitation of our simple hierarchical criterion (a_outer > 10 a_inner). While conservative, this criterion does not account for the destabilizing effects of mass ratios and eccentricities. More sophisticated stability criteria exist in the literature (Mushkin & Katz 2020; Hayashi et al. 2022, 2023), including recent machine learning approaches (Lalande & Trani 2022; Vynatheya et al. 2022), which could better predict unstable configurations. However, our testing indicates that instabilities are relatively rare (∼ 2% of systems) outside the high i_mut and e_outer regimes, primarily occurring when M₃ ≫ (M₁ + M₂) (Mardling & Aarseth 2001). Given our focus on mapping the merger boundary rather than predicting individual system evolution, we deemed this level of contamination acceptable.

For systems predicted to merge by both codes, there were notable quantitative differences in the merger timescales. The secular code typically predicted merger times that differed from N-body results by a factor of 2 − 3, occasionally up to an order of magnitude, with the secular prediction generally being longer. This discrepancy arises primarily from differences in the predicted ZLK oscillation periods, with the secular code typically estimating slightly longer periods. Figure 2 illustrates the comparative evolution of two similar systems as predicted by the two codes. While both codes agree that these systems will merge, the secular code predicts a slightly longer ZLK period. This is visible in the Fourier transform of the eccentricity oscillations, where the peak frequency is slightly lower for the secular code. Despite these phase differences, the qualitative behavior is consistent between the twoapproaches.

Fig. 2. Comparative evolution of inner semimajor axis (top) and inner eccentricity (bottom left), with Fourier transform of inner eccentricity (bottom right) including position of peak frequencies (dotted lines), between N-body (black) and secular (purple) codes for two similar initial configurations: M1 = 68.7 M⊙, M2 = 17.0 M⊙, M3 = 92.0 M⊙, ainner = 88.0 AU, aouter = 6.82 kAU, eouter = 0.670, imut = 78.3°.

4.2. Characterizing the parameter space

Our MCMC exploration efficiently identified the regions in parameter space where systems transition from certainly merging to certainly not merging. We monitored the Gelman–Rubin statistic for convergence and found $\widehat{R}<1.1$ for all parameters after burn-in. In regions with rapid transitions between outcomes (particularly near i_mut ≃ 60°, as developed below), we observed slower mixing, suggesting possible complex structure in the merger boundary. Our adaptive sampling approach partially mitigates this by continuously refining the likelihood surface. After analyzing approximately 15 million different initial configurations, we were able to characterize the merger probability distribution across the 7D parameter space and to identify key trends. Figure 3 shows the fraction of systems that merge as a function of individual parameters and parameter pairs. Merger-favorable regions are quite restricted, primarily concentrated around the specific parameter combinations that maximize the effectiveness of the ZLK mechanism or that enable direct merger via GW emission.

Fig. 3. Fraction of systems that merge in each cell for each pair of parameters, or in each bin for each single parameter. The color ranges from purple (no systems merge) through white (∼50% merger rate) to yellow (all systems merge).

Examining the distribution of merger probability as a function of individual parameters reveals several important trends. For the inner binary masses (M₁ and M₂), we observe that the merger fraction is relatively flat with respect to M₁ for values above approximately 25 M_⊙. However, when M₁ < 25 M_⊙, the merger fraction decreases, likely because both masses of the inner binary are low (due to our constraint that M₂ ≤ M₁), resulting in weaker GW emission that cannot overcome the advantages of a stronger ZLK effect for lower masses. Notably, for a given total inner binary mass, a more asymmetric mass configuration (M₁ ≫ M₂) shows a significantly higher merger fraction than a symmetric one. This can be traced to the octupolar contribution of the ZLK mechanism, which depends on the mass difference M₁ − M₂ (e.g., Naoz et al. 2013). A moreasymmetric configuration can reach a higher maximum eccentricity, leading to a smaller pericenter distance and consequently a stronger GW emission. This finding challenges the intuition from isolated binary merger timescales, where equal masses typically lead to faster mergers at a given separation and total mass (Eq. (1)). Somewhat surprisingly, the tertiary mass (M₃) has minimal impact on merger probability across most of its range. Intuition might suggest that a more massive perturber would induce stronger ZLK oscillations, but our results indicate that this parameter is less critical than others in determining merger outcomes.

The semimajor axes (a_inner and a_outer) show the most dramatic influences on merger probability. Systems with very small inner separations (a_inner < 10 AU) can merge through GW emission alone, without significant assistance from ZLK oscillations. As a_inner increases beyond this threshold, systems increasingly require ZLK-induced high-eccentricity episodes to merge within the age of the Universe. The outer separation exhibits a strong negative correlation with merger probability: smaller a_outer values lead to much higher merger fractions, due to stronger perturbations on the inner binary. Systems with a_outer ≳ 10 a_inner (near the stability limit of hierarchical triples) almost universally merge. Furthermore, the outer eccentricity e_outer exhibits a strong positive correlation with merger probability. Higher eccentricities bring the tertiary closer to the inner binary at pericenter, enhancing the strength of ZLK oscillations and leading to higher inner binary eccentricities. This effect becomes particularly pronounced for e_outer > 0.6.

Similarly, the mutual inclination i_mut shows a generally positive correlation with merger probability, consistent with theoretical expectations that higher inclinations produce stronger ZLK effects (e.g., Katz et al. 2011). Interestingly, we observe a nonmonotonic feature with peaks in merger probability centered around i_mut = 54° and 64°, separated by a trough at i_mut = 60°. This feature is statistically robust across different values of other parameters, but a clear physical explanation for this specific pattern remains elusive. It may in fact reflect the complex, potentially fractal nature of the boundary between merging and nonmerging configurations in the three-body problem (in which regular and chaotic regions coexist at all scales, as shown by Trani et al. 2024). Additionally, as discussed by Hamers (2021), departing from the test-particle limit introduces asymmetries in eccentricity excitation, which can manifest as such irregularfeatures.

To further characterize the parameter space, we classified nonmerging systems based on their final state after evolving for 14 Gyr. Specifically, we examined whether the inner binary’s separation decreased (indicating GW emission) and whether the inner eccentricity changed from its initial value (indicating ZLK oscillations). This classification resulted in four distinct categories of nonmergingsystems:

• GW, ZLK (26% of nonmerging systems): Systems that emit GWs and experience ZLK oscillations, but do not merge within 14 Gyr. These systems are typically found in regions with moderate inner separations (10 − 50 AU) and outer separations smaller than about 2000 AU.

• GW, no ZLK (10%): Systems that emit GWs without significant ZLK oscillations. These are predominantly systems with very small inner separations (< 10 AU), where the binary is already close enough for GW emission to be significant, but where the ZLK mechanism is suppressed due to relativistic precession.

• No GW, ZLK (64%): Systems that experience ZLK oscillations, but do not emit significant GWs. As expected, these systems have very low inner masses, or alternatively large inner separations where even the maximum eccentricity reached during ZLK cycles does not bring the pericenter close enough for substantial GW emission.

• No GW, no ZLK (< 1%): Systems that neither experience ZLK oscillations nor emit significant GWs. These rare cases typically occur when the mutual inclination is near the lower threshold for ZLK oscillations and other parameters are unfavorable.

Figure 4 visualizes the distribution of these categories across the parameter space, revealing clear regional trends. For each parameter, distinct patterns emerge in how they influence the distribution of nonmerging system types. The inner separation a_inner is particularly decisive. The “GW, no ZLK” systems dominate at very small separations (a_inner < 10 AU), especially forfar-orbiting perturbers (triggering minimal ZLK resonances) and for systems with more massive inner binaries, where GW emission is strong enough to drive orbital decay without requiring ZLK-induced eccentricity. As a_inner increases, we observe a transition to “GW, ZLK” systems, which peak at moderate separations (10 − 50 AU) where relativistic precession is weaker and allows ZLK oscillations to operate effectively, while still permitting significant GW radiation. At larger separations (> 50 AU), “no GW, ZLK” systems become prevalent, as the separations are too large for efficient GW emission regardless of ZLK-induced eccentricity. The outer separation a_outer also exhibits a clear gradient effect. When it is close to the stability limit (a_outer ≳ 10 a_inner), systems predominantly merge (Fig. 3). As a_outer increases, we first encounter a transition zone of mostly “GW, ZLK” systems, followed by a region with “No GW, ZLK” systems catching up at larger separations where the perturbing influence diminishes, reducing the maximum eccentricity achieved during ZLK cycles.

Fig. 4. Distribution of three categories of nonmerging systems in the parameter space. Each category is represented by a base color: cyan (GW, ZLK), magenta (no GW, ZLK), and yellow (GW, no ZLK). For each pair of parameters, cells are color-coded according to the combination of the three base colors resulting from the fraction of systems corresponding to the three categories.

The mutual inclination i_mut strongly influences the system category, with a noticeable pattern: “no GW, ZLK” systems prevail for the lowest inclinations, except when conditions are too unfavorable for ZLK (typically for circular perturbers and short inner orbits), and gradually fade out as i_mut is increased. This is related to the fact that the highest ZLK-induced eccentricities, required for efficient GW emission, can only be achieved with orbits that are tilted enough. Around i_mut ≳ 40°, “GW, no ZLK” systems are more common than elsewhere, as these inclinations are near the threshold required for ZLK oscillations. As inclination increases toward 90°, “GW, ZLK” systems become increasingly dominant, reaching over 50% of nonmerging systems at the highest inclinations. Again, this reflects the greater efficiency of ZLK oscillations at architectures closer to perpendicular orbital configurations. For e_outer, higher eccentricities show fewer “No GW, ZLK” and more “GW, ZLK” systems because the closer pericenter passages of the tertiary body enhance the GW radiation.

The masses of the three BHs undergo more subtle changes. For the inner binary masses M₁ and M₂, systems with higher masses show increased fractions of “GW, no ZLK” and “GW, ZLK” categories, as expected from the stronger GW emission from more massive binaries. However, the mass asymmetry in the inner binary is also important; when M₁ ≫ M₂, the octupolar contributions to ZLK oscillations are enhanced, leading to higher fractions of systems experiencing ZLK oscillations. The mass imbalance effect is clear for the lowest mutual inclinations: such regions are usually populated by “no GW, ZLK” systems, but they transition to “GW, no ZLK” for the highest values of M₂. The latter correspond to zones where the octupolar effect is the most modest (due to our constraint that M₂ ≤ M₁), hence suppressing ZLK cycles. The tertiary mass M₃ shows a relatively uniform distribution of categories across its range, confirming its secondary importance with respect to geometric parameters.

These detailed characterizations of the parameter space provide valuable insights into which initial configurations are most likely to produce BH mergers via the hierarchical triple channel. The clear regional patterns we identify can inform models of GW progenitor formation and help interpret the growing catalog of GW detections. We note that while we present here the results found using the natural mass and separation parameters, we verified that alternative parameterizations (e.g., M_inner = M₁ + M₂, q₂₁ = M₂/M₁, q₃ = M₃/M_inner, α = a_inner/a_outer) yield consistent conclusions about the key trends governing merger probability.

4.3. Neural network performance

The neural network we developed achieves an overall accuracy of 94.7% on the test set of 2 632 159 samples (20% of the full dataset), demonstrating its ability to reliably predict merger outcomes without performing full dynamical simulations. This high accuracy is particularly remarkable given the complex, nonlinear nature of the ZLK mechanism and its interaction with GW emission. Table 1 presents the confusion matrix of our model’s predictions. The test set contains a somewhat unbalanced distribution, with approximately 70% nonmerging systems, reflecting the natural distribution found while exploring the parameter space, including the rejected samples in the MCMC. Despite this imbalance, the model performs well on both classes, showing no critical differences between false positives (50 828) and false negatives (89 798).

The model’s performance, illustrated by the receiver operating characteristic (ROC) curve in Fig. 5, can be characterized by several standard metrics:

Fig. 5. ROC curve of the neural network. The AUC is shown, as is the diagonal (pure chance).

• Positive predictive value (precision): 93.4% – when the model predicts a merger, it is correct 93.4% of the time;

• Negative predictive value: 95.2% – when the model predicts no merger, it is correct 95.2% of the time;

• True positive rate (recall): 88.9% – the model correctly identifies 88.9% of all merging systems;

• True negative rate (specificity): 97.2% – the model correctly identifies 97.2% of all nonmerging systems;

• F1 score: 91.1% – the harmonic mean of precision and recall, providing a balanced measure of the model’s performance;

• Area under the ROC curve (AUC): 99.0% – this metric reflects the model’s ability to distinguish between merging and nonmerging classes across all classification thresholds.

A particularly valuable aspect of our neural network implementation is the confidence measure (Eq. (10)), which provides a robust assessment of prediction reliability. This measure allowed us to analyze the relationship between the model’s confidence and its accuracy, revealing a strong positive correlation. Figure 6 illustrates this relationship by showing the distribution of true and false predictions across different confidence levels. Analysis of the confidence distribution reveals that 69.6% of all test samples have a confidence measure above 0.99, and 79.5% have confidence above 0.9. More importantly, the accuracy is strongly correlated with confidence: for predictions with $c(\widehat{y})>0.9$, the accuracy increases dramatically to 99.7%. This includes 69.3% of true positive and 89.4% of true negative predictions having confidence above 0.9, while only 3.0% of false positives and 6.5% of false negatives have such high confidence. For these high-confidence systems, the model achieves a positive predictive value of 99.7% and a negative predictive value of 99.6%. In practical terms, this means that for approximately 80% of all systems, we can predict the merger outcome with near certainty. For the remaining 20% of systems with lower confidence predictions, the model still achieves approximately 75% accuracy.

Fig. 6. Histograms of true positives, true negatives, false positives, and false negatives as functions of prediction confidence $c(\widehat{y})$.

The singular performance of our neural network, particularly its ability to provide reliable confidence measures, makes it an ideal tool for rapid population synthesis studies. Researchers can use this model to quickly classify millions of systems, focusing detailed simulations only on the subset of cases where the neural network indicates uncertainty. The prediction phase is remarkably efficient, requiring only milliseconds per system or a few seconds for millions of systems on standard computational hardware, representing an improvement of many orders of magnitude compared to full simulations.

5. Discussion and conclusion

Our study has several caveats that should be acknowledged. While computationally efficient, our secular approach has inherent limitations. First, we do not include semisecular corrections that become important in the most compact hierarchical configurations (as reviewed by Tremaine 2023). These corrections, particularly the single-averaged terms that account for variations on the inner orbital timescale, can significantly enhance eccentricity excitation and potentially reveal additional mergingsystems (Mangipudi et al. 2022). Second, our validation against N-body simulations revealed poor accuracy of the secular formalism for near-perpendicular mutual inclinations on the one hand, and for near-unity outer eccentricities where the fundamental assumption of well-separated orbital timescales is violated on the other hand. The discrepancies in these regimes may partly stem from the missing semisecular corrections becoming crucial when orbital timescales approach commensurability. However, these same regions also exhibit high rates of dynamical instability (12% for i_mut ∈ [80° , 90° ] and 37% for e_outer > 0.9), making their exclusion from our parameter space prudent regardless of the theoretical framework employed. Third, our simple hierarchical stability criterion permits some (∼2%) dynamically unstable configurations in the remaining parameter space. While more sophisticated stability criteria incorporating masses and eccentricities are available, we opted for computational efficiency given that our primary goal is mapping broad merger probability trends rather than predicting individual systemevolution.

Despite these limitations, our 87% agreement rate with N-body simulations (excluding the problematic high i_mut/e_outer regimes) confirms that the secular approximation captures the essential dynamics for the vast majority of our parameter space. Nevertheless, it should be noted that we left several parameters unexplored. We actually did not vary the arguments of pericenter or longitudes of ascending node. BH spins were also not included in our dynamical model, though spin–orbit and spin–spin couplings could influence evolution in some systems (Kidder et al. 1993; Hartl & Buonanno 2005; Bern et al. 2021). In addition, we started our simulations after BH formation, ignoring the effects of stellar evolution and the ZLK mechanism during the progenitor phase. Finally, our simulations did not follow the final inspiral phase in detail, thus preventing accurate gravitational waveform predictions.

Still, we conducted in this work the most comprehensive exploration to date of the hierarchical triple BH parameter space, analyzing the conditions under which such systems can produce inner binary mergers within the age of the Universe through an intensive study scrutinizing around 15 million configurations. Our investigation reveals that the ZLK mechanism provides a viable and efficient pathway for overcoming the initial separation problem in compact object binaries, particularly under specific parameter combinations that maximize the mechanism’s effectiveness. The general patterns we were able to pinpoint are that systems with asymmetric inner binary masses benefit from enhanced octupolar contributions to the ZLK mechanism, leading to stronger eccentricity oscillations and consequently more efficient GW emission. When combined with moderate inner separations (where relativistic precession does not overwhelm the ZLK effect), small outer separations (providing stronger perturbations), and high outer eccentricities (bringing the perturber closer at pericenter), these configurations consistently produce mergers within cosmological timescales.

The limit between merging and nonmerging configurations in our 7D parameter space exhibits complex structures. While we can identify broad statistical trends, the three-body problem’s inherently chaotic nature means this boundary is not smooth, but likely exhibits fractal characteristics. As demonstrated by Trani et al. (2024), regular and chaotic trajectories are interwoven at all scales in a multifractal pattern, with regular regions occupying a substantial portion of the phase space depending on the initial configuration. In the hierarchical limit with quadrupole approximation, the system can be integrable (e.g., Kinoshita & Nakai 1999), but relaxing these assumptions (as done in this work) introduces chaos that may manifest as fractal boundaries between different outcomes.

Our methodological approach, coupling secular approximation with an adaptive MCMC technique, has globally proven effective for navigating this high-dimensional parameter space. By targeting the transition regions between certainly merging and certainly nonmerging configurations, we mapped the probability landscape in fine detail. This mapping not only characterizes where mergers occur, but also identifies distinct categories of nonmerging systems, providing additional insights into the interplay between ZLK oscillations and GW emission. Nonetheless, the true boundary likely has a fractal dimension smaller than the embedding dimension of our parameter space. This means that arbitrarily fine sampling would continue to reveal new structure, with pockets of regular behavior embedded within chaotic regions and vice versa. Future work could investigate the fractal dimension of these boundaries and their impact on merger rate predictions.

Perhaps the most practical outcome of our work is the neural network model³, which provides very fast predictions of merger outcomes with good accuracy. This tool transforms what would otherwise be computationally prohibitive population synthesis studies into feasible investigations. With an area under the ROC curve of 99%, 95% overall accuracy, and nearly perfect predictions (99.7%) for the majority of systems, researchers can now efficiently identify promising configurations for further detailed study without the computational burden of full dynamical simulations.

As GW detectors continue to increase in sensitivity and the catalog of detected signals grows, the patterns we have identified can help constrain the astrophysical origins of observed mergers. Our results suggest that hierarchical triple systems contribute a distinct population of merging BHs, potentially distinguishable by their mass ratios and orbital characteristics from those produced by other formation channels. Future work connecting our parameter space exploration to realistic initial conditions from stellar evolution models will further refine our understanding of the relative importance of this pathway in nature. Ultimately, this study advances our understanding of the complex gravitational dynamics that shape the Universe’s most extreme environments, bringing us closer to deciphering the origin stories written in the GWs that ripple through spacetime.

Acknowledgments

We are grateful to our referee, Alessandro Alberto Trani, for improving the quality of our article through a thoughtful review. We offer our thanks to Georges Meynet for his support and for fruitful conversations. We made use of the Claude AI assistant (Anthropic, 2024) for code optimization. We thank the University of Geneva HPC team for computational resources.

References

All Tables

All Figures

	Fig. 1. Generic orbital configuration of a hierarchical triple system, as investigated in this study.
In the text

Fig. 2. Comparative evolution of inner semimajor axis (top) and inner eccentricity (bottom left), with Fourier transform of inner eccentricity (bottom right) including position of peak frequencies (dotted lines), between N-body (black) and secular (purple) codes for two similar initial configurations: M1 = 68.7 M⊙, M2 = 17.0 M⊙, M3 = 92.0 M⊙, ainner = 88.0 AU, aouter = 6.82 kAU, eouter = 0.670, imut = 78.3°.

In the text

	Fig. 3. Fraction of systems that merge in each cell for each pair of parameters, or in each bin for each single parameter. The color ranges from purple (no systems merge) through white (∼50% merger rate) to yellow (all systems merge).
In the text

	Fig. 4. Distribution of three categories of nonmerging systems in the parameter space. Each category is represented by a base color: cyan (GW, ZLK), magenta (no GW, ZLK), and yellow (GW, no ZLK). For each pair of parameters, cells are color-coded according to the combination of the three base colors resulting from the fraction of systems corresponding to the three categories.
In the text

	Fig. 5. ROC curve of the neural network. The AUC is shown, as is the diagonal (pure chance).
In the text

	Fig. 6. Histograms of true positives, true negatives, false positives, and false negatives as functions of prediction confidence $c(\widehat{y})$.
In the text