We explore the stability and fate of gravitational triple systems comprising a central massive body and a tight binary of less massive pairs. Our present purpose is two fold 1 to improve the Hilltype stability criterion for the binary in those systems and 2 to examine the fate of the triple systems after the binary breakup with particular attention to the effects of the eccentricities of the inner and outer orbits. We perform direct Newtonian Nbody simulations over much longer integration times than previous studies which is essential to determine the eventual fate of those systems statistically in a reliable fashion. We obtain an empirical fitting formula of the binary stability boundary that incorporates effects of the inner and outer eccentricities the mutual inclination of the inner and outer orbits and the mass ratios of the three bodies. We also find that those triple systems are stable for a much longer timescale after the binary breakup and that their final fates ejection of the outer body merger to the central massive body and collision of two less massive bodies are very sensitive to the initial outer eccentricity.
['Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Earth and Planetary Astrophysics']
Stability and fate of hierarchical triples comprising a central massive body and a tight binary in eccentric orbits
{'No Header': 'Draft version September 11, 2024  \nTypeset using L A T E X preprint style in AASTeX631', 'comprising a central massive body and a tight binary in eccentric orbits': 'Toshinori Hayashi , 1 Alessandro A. Trani , 2 and Yasushi Suto 3, 4, 5  \n- 1  \nYukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan 2 Niels Bohr Institute, University of Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark 3 Research Institute, Kochi University of Technology, Tosa Yamada, Kochi 782-8502, Japan 4 Department of Physics, The University of Tokyo, Tokyo 113-0033, Japan  \n5 Research Center for the Early Universe, School of Science, The University of Tokyo, Tokyo 113-0033, Japan  \n(Received August 31, 2024; Revised; Accepted)  \nSubmitted to ApJ', 'ABSTRACT': 'We explore the stability and fate of gravitational triple systems comprising a central massive body and a tight binary of less massive pairs. Our present purpose is two fold; (1) to improve the Hill-type stability criterion for the binary in those systems, and (2) to examine the fate of the triple systems after the binary break-up, with particular attention to the effects of the eccentricities of the inner and outer orbits. We perform direct Newtonian N-body simulations over much longer integration times than previous studies, which is essential to determine the eventual fate of those systems statistically in a reliable fashion. We obtain an empirical fitting formula of the binary stability boundary that incorporates effects of the inner and outer eccentricities, the mutual inclination of the inner and outer orbits, and the mass ratios of the three bodies. We also find that those triple systems are stable for a much longer timescale after the binary break-up, and that their final fates (ejection of the outer body, merger to the central massive body, and collision of two less massive bodies) are very sensitive to the initial outer eccentricity.  \nKeywords: celestial mechanics - (stars:) binaries (including multiple): close - stars: black holes', '1. INTRODUCTION': "The stability and fate of gravitational triple systems is one of the long-standing and challenging questions in mathematical physics and astronomy. Many previous authors have approached the problem using a variety of different approximations and methodologies, including Eggleton & Kiseleva  \nCorresponding author: Toshinori Hayashi  \ntoshinori.hayashi@yukawa.kyoto-u.ac.jp  \n(1995); Holman & Wiegert (1999); Mardling & Aarseth (1999, 2001); Georgakarakos (2013); Grishin et al. (2017); He & Petrovich (2018); Myllari et al. (2018); Wei et al. (2021); Lalande & Trani (2022); Tory et al. (2022); Vynatheya et al. (2022); Hayashi et al. (2022, 2023); Zhang et al. (2023), among others.  \nIn particular, a lot of attention has been paid to a hierarchical triple configuration, in which two of them form a tight binary and interact with the tertiary object. For instance, Mardling & Aarseth (2001) (hereafter, MA01) considered their stability combining the chaos theory and direct numerical simulations, and obtained the following stability criterion:  \na out (1 -e out ) a in > 2 . 8 ( 1 -0 . 3 i mut 180 · )[( 1 + m 3 m 1 + m 2 ) (1 + e out ) √ 1 -e out ] 2 / 5 . (1)  \nIn the above expression, a out and a in are the semi-major axes of the outer and inner orbits, e out is the orbital eccentricity of the outer orbit, i mut is the mutual inclination between the outer and inner orbits, m 3 is the mass of the tertiary, and m 1 and m 2 are the masses of the inner binary.  \nThe criterion (1) (implicitly) assumes that the tertiary mass m 3 is at most comparable to m 1 and m 2 , and has been frequently applied to the case of m 3 ≪ m 12 ≡ m 1 + m 2 . In what follows, we refer to such configurations as HT-P (Hierarchical Triple Planet-type) in which a tertiary orbits an inner massive binary. The classic criterion (1) has been tested and improved recently by Vynatheya et al. (2022); Hayashi et al. (2022, 2023), for instance.  \nAnother hierarchical triple configuration, which we refer to as HT-S (Hierarchical Triple Satellitetype) below, is a central massive object orbited by a less massive binary. The dynamics of HTS configurations is important in understanding the fate of a variety of physically interesting triple systems. For instance, the enhanced merging rate of Binary Black-Holes (BBHs) around supermassive BH (SMBH) (e.g. Li et al. 2015; VanLandingham et al. 2016) may be a major target of future spacebased gravitational wave detectors (e.g. Xuan et al. 2023). Extrasolar binary-planet systems also belong to this configuration, which have been theoretically predicted/discussed (e.g. Ochiai et al. 2014; Lewis et al. 2015), but not yet detected.  \nThe present paper focuses on the stability and fate of HT-S systems comprising a central massive body of mass m 0 orbited by an initially tightly bound binary of mass m 1 and m 2 . We examine a break-up condition of the binary and obtain its stability boundary as a function of the initial parameters of the systems. Then, we consider the subsequent fate of systems after the binary breakup, keeping in mind the application to black-hole triples. To avoid confusion, we use m 0 for the massive tertiary in HT-S systems, while we use m 3 for the less massive tertiary in the HT-P systems. Our current analysis assumes purely Newtonian gravity for simplicity and clarity, and the effect of general relativity will be studied in a forthcoming paper.  \nThe Hill stability condition is particularly useful in considering the stability of HT-S systems. Historically, the Hill stability was first derived for the motion of moon, applying the conservation of the Jacobi integral in restricted three-body problems (e.g. Hill 1878). Later, it was extended to more general (non-restricted) three-body problems using Sundman's inequality, which determines a required condition among momenta for a gravitational multi-body system, and the sufficient conditions for some cases were derived (e.g. Marchal & Bozis 1982).  \nWhile the Hill stability is rigorously defined in those papers, the condition is roughly understood as the competition between the gravitational tidal force of the binary due to the central object and the mutual gravitational attraction of the binary pairs.  \nConsider two objects of masses m 1 and m 2 orbiting a more massive central object of mass m 0 . If both the inner and outer orbits are circular, under the test particle limit ( m 2 ≪ m 1 < m 0 ), a binary companion m 2 orbiting m 1 is stably bound to m 1 if its semi-major axis a 12 is less than the Hill radius a Hill :  \na 12 < a Hill ≡ a 012 ( m 1 3 m 0 ) 1 / 3 . (2)  \nwhere a 012 is the semi-major axis between m 1 -m 2 binary and m 0 . The above stability condition is rewritten in terms of the inner and outer orbital periods as  \nP 012 > √ 3 P 12 . (3)  \nwhere P 012 and P 12 are the periods of outer and inner orbits, respectively. In other words, the Hillstability is roughly equivalent to the condition that the outer orbital period is longer than the inner one. Note that, for clarity, we use the labels '12' and '012' for HT-S, instead of 'in' and 'out' of HT-P, unless otherwise specified.  \n̸  \n̸  \nWhile the Hill stability condition (2) or (3) gives approximately the stability condition for the HT-S systems, it needs to be generalized to the cases of non-circular ( e 12 = 0 and/or e 012 = 0) orbits, the finite mass of m 2 , and non-coplanar orbits.  \nFor instance, Grishin et al. (2017) (hereafter GPZM17) generalized the stability condition (2) by considering the mutual inclination dependence for initially near circular inner and outer orbits:  \na 012 a 12 > 1 3 1 / 3 ( 3 m 0 m 1 + m 2 ) 1 / 3 (cos i mut + √ 3 + cos 2 i mut ) 2 / 3 ×    1 (cos 2 i mut > 3 / 5) 9 -5 cos 2 i mut 6 (cos 2 i mut ≤ 3 / 5) . (4)  \nGPZM17 derived the above criterion by replacing a 12 with a 12 (1 + 0 . 5 e 2 max ) under the von ZeipelKozai-Lidov (ZKL) oscillations (von Zeipel 1910; Kozai 1962; Lidov 1962), where e max is the maximum value of e 12 attained during the ZKL cycle. GPZM17 confirmed that the condition (4) well captures the strong destabilization of HT-S systems with near-polar orbits using numerical simulations.  \n̸  \n̸  \nThe present paper aims to further generalize the condition (4) for the HT-S systems of initially non-circular orbits ( e 12 = 0 and e 012 = 0), as pioneered by MA01 for the HT-P systems. This is shown in Figure 1, where we plot the MA01 stability condition (1) for the coplanar HT-P systems (red curves). Because GPZM17 focused on initially circular HT-S systems, their result (4) is valid for e 012 = 0 alone. Figure 1 assumes that the initial inner orbits are circular ( e in (= e 12 ) = 0), but shows the mass dependence using different line types for comparison.  \nIf the Hill stability is purely determined by the amplitude of the instantaneous tidal force due to the central massive body, one naively expects that the non-circular effect would be incorporated simply by replacing the initial semi-major axes a 12 and a 012 with their initial apocenter distance a 12 (1 + e 12 ), and pericenter distance a 012 (1 -e 012 ), respectively. In reality, however, this simple procedure completely neglects the evolution of the orbital elements. This is why we explore the non-circular effect using a series of systematic numerical simulations below.  \nAs schematically illustrated in Figure 2, the main purpose of this paper is basically two-fold. The primary purpose, derivation of a generalized stability boundary for the HT-S systems that we  \nFigure 1. Comparison of stability criteria for HT-P ( m 12 /m 3 = 2 , 10 2 , red curves) and HT-S ( m 12 /m 0 = 2 , 10 -2 , 10 -4 , black dots) in e out -a out (1 -e out ) /a in plane, where m 12 = m 1 + m 2 . For HT-S, a in , a out and e out are interpreted as a 12 , a 012 and e 012 , respectively. Note that the stability criterion of HT-P (GPZM17) focuses on circular systems, and horizontal gray lines are simple extrapolations assuming e 012 dependence is all absorbed in a 012 (1 -e 012 ).  \n<!-- image -->  \neout  \ndiscussed in the above, corresponds to the classification 1. We attempt to extend the GPZM17 stability criterion by including initial eccentricity, mutual inclination, and mass ratio dependences, using long-term N-body integrations. We emphasize that the direct N-body simulations, instead of the orbit-averaged secular simulations, are essential to explore the Hill-stability boundary between (A) and (B).  \nThe secondary purpose is to find the outcome of the Hill unstable HT-S systems. As indicated in the classification 2, there are four possible outcomes after the initial tight binary is broken. Their branching ratio is somewhat dependent on the integration time of the simulations. We practically perform the simulations up to 10 8 P 12 , and estimate the fraction of the four different outcomes evaluated at the epoch as a function of e 012 . We here note that Toonen et al. (2022) perform a pioneering study about the evolution and fate of HT-P type stellar systems, considering stellar evolution and various triple population models as well. We assume point-mass objects here and do not consider any astronomical evolution except the Newtonian dynamics.  \nThe rest of the paper is organized as follows. In section 2, we first describe initial setup of numerical simulations, and how we evaluate the breaking condition of initial binary. Then, section 3 shows the resulting Hill-stability criterion, including initial eccentricity, mutual inclination, and mass ratio  \nFigure 2. Schematic illustration of the contents in this paper. Classification 1 checks the breaking condition of SMBH-BBH triples, and Classification 2 checks the outcome fraction and timescale of Hill-unstable triples. Note that which of m 1 or m 2 becomes m in is in general time-dependent during evolution. While channel (A) can be computed using secular perturbation, channel (B) requires direct N-body simulations.  \n<!-- image -->  \ndependences. We also systematically check the initial phase dependence about the Hill stability, and show the result in the appendix. Next, in section 4, we show a demonstrative example of outcome fraction and timescale, focusing on the Hill-unstable systems with different mutual inclinations. Finally, section 5 summarizes the conclusion of this paper, and discusses possible future prospects.", '2. METHOD': 'We carry out a series of N-body simulations using TSUNAMI (see Trani & Spera 2023), and explore the parameter space for the initial values of e 12 and e 012 in particular. We also consider a set of different mass ratios for the three bodies, and four representative mutual inclinations ( i mut = 0 · , 5 · , 90 · , 180 · ).  \nWe determine the stability boundary of those systems on e 012 -a 012 (1 -e 012 ) /a 12 plane, following Hayashi et al. (2022, 2023). For a given set of values for e 12 , e 012 , a 12 , i mut , m 1 , m 2 , and m 0 , we start a simulation with a sufficiently small value of a 012 so that the m 1 -m 2 binary quickly breaks due to the Hill instability. Then, we gradually increase the value of a 012 , therefore a 012 (1 -e 012 ) /a 12 , and record the corresponding break time T break of the m 1 -m 2 binary. The procedure is repeated up to our adopted integration time limit t int = 10 8 P 12 . The stability boundary is defined by the critical value of a 012 (1 -e 012 ) /a 12 with which the binary does not break within the integration time limit. To avoid a possible spurious result due to the fluctuation, we check a larger value of a 012 (1 -e 012 ) /a 12 , and if the binary does not break again, we define this critical value as the stability boundary. The robustness of the resulting boundary is separately examined by running realizations of different initial phases later (see the appendix).  \nAll the above parameter survey is done for a given set of the initial phases (mean anomalies and pericenter arguments), but the initial-phase dependence is checked later for a sample of specific orbital configurations. Throughout this paper, we fix P 12 as 1 yr unless otherwise stated, but this does not virtually affect our result thanks to the scaling relation in Newtonian gravity.  \nIn this procedure, the definition of the binary break is the most crucial, and we basically apply the disruption condition in Hayashi et al. (2023) as follows. At each numerical timestep during the simulation, we compute the osculating Kepler orbital elements for bounded pairs of particles; ( m 1 , m 2 ), ( m 0 , m 1 ) and ( m 0 , m 2 ) at most. Then, we define that the initial binary pair ( m 1 , m 2 ) is broken if  \n0 < a 01 ( t ) < a 12 ( t ) or 0 < a 02 ( t ) < a 12 ( t ) or a 12 ( t ) < 0 . (5)  \nThe above conditions basically reflect the configuration of (B) in Figure 2. We adopt the semi-major axes, instead of the binary binding energies, in the above definition of the binary break. This is because the amplitude of binding energies is always dominated by the central massive body in the present case ( m 0 ≫ m 1 ≈ m 2 ).  \nGPZM17 adopted a binary break when either (1) the eccentricity of a binary exceeds unity or (2) the distance of binary exceeds 3 r Hill . We prefer the condition (5), rather than that adopted by GPZM17 for the following reasons. First, the osculating eccentricity may not be stable especially when the binary becomes close to a break-up configuration. Second, to use the Hill-stability radius as a measure of the break-up may lead to a circular argument, because our current purpose is to find the improved criterion for the Hill stability. In practice, however, we numerically made sure that the two different definitions change the stability boundary, i.e., the value of the y-axis in Figure 1, usually within tens of percent as shown in the next section.', '3. IMPROVED FORMULA FOR STABILITY BOUNDARIES OF THE HT-S SYSTEMS': 'In this section, we derive an empirical formula for stability boundaries of the HT-S system that improves the GPZM17 stability criterion (4). We examine the inner and outer eccentricity dependence for 4 different cases of the mutual inclination i mut = 0 · (coplanar-prograde), 5 · (near-coplanarprograde), 90 · (orthogonal), and 180 · (coplanar-retrograde). We considered i mut = 5 · in addition to 0 · because the final fate after the binary break is sensitive to the exact or near coplanarity of the orbits, but the stability boundary turned out to be almost identical. Furthermore, we confirmed the m 0 /m 12 -dependence in the GPZM17 criterion is valid for a wide range of m 0 /m 12 , and make sure that the result is statistically robust against the initial orbital phases of the three bodies (see the appendix). Adopting the i mut -dependence discovered by GPZM17, our stability criterion is written as follows:  \nΥ ≡ ˜ r 012 ˜ r 12 > Υ crit ≡ ( m 0 m 1 + m 2 ) 1 / 3 (cos i mut + √ 3 + cos 2 i mut ) 2 / 3 h ( e 012 ) , (6)  \nwhere  \n˜ r 012 ≡ a 012 (1 -e 012 ) , (7) ˜ r 12 =      a 12 (1 + 0 . 5 e 2 12 ) ( i mut = 0 · , 5 · ) a 12 [1 + 0 . 5 e 2 max ( i mut )] ( i mut = 90 · ) a 12 (1 + e 12 ) ( i mut = 180 · ) , (8)  \nh ( e 012 ) =1 . 5 + 0 . 6 e 012 ,  \n(9)  \nand e max is the maximum binary eccentricity e 12 under the ZKL oscillations. In the test-particle limit,  \ne max = √ 1 -5 3 cos 2 i mut (cos 2 i mut ≤ 3 / 5) , (10)  \nrecovering the part of the i mut -dependence in the GPZM17 criterion (4):  \n1 + 0 . 5 e 2 max ( i mut ) = 9 -5 cos 2 i mut 6 . (11)  \nHere, we note that Y is usually defined as a out (1 -e out ) /a in (1 + e in ) (e.g. Eggleton & Kiseleva 1995; Vynatheya et al. 2022), and we define Υ as the extension of this quantity.  \nThe next subsections will describe how we obtain the parameter dependence in the stability boundary formula (6). We run a series of simulations for i mut = 0 · , 5 · , 90 · , and 180 · , fixing initial phases ( ω 12 , ω 012 ,M 12 , M 012 ) as (180 · , 0 · , 30 · , 45 · ), where M 12 , M 012 , ω 12 , and ω 012 are the inner mean anomaly, outer mean anomaly, inner pericenter argument, and outer pericenter argument, respectively (see Figure 3). We use the Jacobi coordinate system, and the invariant plane, therefore fixing Ω 12 and Ω 012 as 180 · and 0 · , respectively.  \nSince our binary breaking timescale T break is evaluated in units of the initial value of P 12 , the result is also scalable with respect to the specific value of P 12 , and we set P 12 = 1 yr without loss of generality.  \nWe examine the dependence of the stability boundary on the eccentricities e 12 and e 012 in § 3.1 and § 3.2 for m 0 = 10 6 M ⊙ and m 1 = m 2 = 10 M ⊙ , with a triple system comprising a massive BH and a stellar mass binary black hole in mind. Since we neglect the effect of general relativity and do not consider the merger condition here, the result is dependent on their mass ratio alone. The dependence on their mass ratio is examined in § 3.3, and the sensitivity to the initial phases is discussed in appendix A.', '3.1. Inner eccentricity dependence': 'For the initial configuration described in the above, we examine 25 models with e 012 = 0 , 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8, and e 12 = 0 , 0 . 2 , 0 . 4 , 0 . 6 , 0 . 8. For each model, we gradually increase the initial value of a 012 (1 -e 012 ) /a 12 and find the stability boundary at which m 1 -m 2 binary does not break up before 10 8 P 12 .  \nThe result is plotted in Figure 4 in which the binary break-up timescale of each system is indicated as a color-coded filled circle. If a system does not break up within the integration time limit (10 8 P 12 ), the system is indicated by a cross symbol. The horizontal axis of the figure is e 012 , and five sequences centered at the value of e 012 correspond to the results that we slightly shifted for visual clarity according to the value of e 12 (0, 0 . 2, 0 . 4, 0 . 6 and 0 . 8 from left to right). For reference, the dashed black lines in those panels indicate the GPZM17 criterion obtained for e 012 = 0.  \nThe top four panels plot a 012 (1 -e 012 ) /a 12 (1 + 0 . 5 e 2 12 ), the pericenter distance between the binary and the central massive object in units of the orbit-averaged distance of the binary, for different i mut . While the effect of e 12 is reasonably absorbed in the above scaling for i mut = 0 · and 5 · , there remains a clear systematic trend for i mut = 90 · and 180 · .  \nInstead, we found that the residual e 12 -dependence is well absorbed by changing the vertical axis as a 012 (1 -e 012 ) /a 12 (1 + 0 . 5 e 2 max ) for i mut = 90 · and a 012 (1 -e 012 ) /a 12 (1 + e 12 ) for i mut = 180 · ; see the bottom two panels in Figure 4.  \nFigure 3. The initial configuration of the inner and outer orbits.  \n<!-- image -->  \nThe scaling with respect to e 12 for (near)-coplanar prograde cases implies that the binary break-up due to the Hill-type instability is not instantaneous generally, but happens in a somewhat cumulative fashion. This is reasonable especially around the stability/instability boundary where the outer orbital period is supposed to be much longer than the inner orbital period as indicated from the conventional Hill stability condition (3).  \nThe scaling for i mut = 90 · was already suggested by GPZM17, and can be understood as the orbitaveraged binary distance should be computed from e max , instead of the initial value of e 12 , due to the ZKL oscillation. We note that Vynatheya et al. (2022) suggest that a similar approach works well for HT-P systems.  \nIn contrast, the physical origin of the scaling for a coplanar retrograde is not clear. It seems as if the instability happened instantaneously when the temporal tidal interaction of the three bodies is strongest, but it is unlikely the case in reality. Thus, we leave the explanation to the future work, and simply show it as the empirical result at this point.', '3.2. Outer eccentricity dependence': 'Consider next the e 012 -dependence of the stability boundary. For that purpose, we replot Figure 4 in such a way that the vertical axis is now normalized by the GPZM17 criterion (4). This is equivalent to determine the functional form of h ( e 012 ) in our final expression (6). The result is shown in Figure 5, indicating that a simply linear function well captures the e 012 -dependence (red lines):  \nh ( e 012 ) = 1 . 5 + 0 . 6 e 012 . (12)  \nWe emphasize that the i mut -dependence of the GPZM17 result is in good agreement with our numerical results for i mut = 0 · , 5 · , 90 · and 180 · .  \n/  \nFigure 4. Break-up time distribution in e 012 - distance ratio plane. Six panels show the results for i mut = 0 · , 5 · , 90 · , and 180 · . Each group of five represents the same value of e 012 , but different values of e 12 , shifted horizontally in the range of 0, 0 . 2, 0 . 4, 0 . 6 and 0 . 8 from left to right. The upper four panels use the average distances of inner orbits, although lower two panels use different distance measures.  \n<!-- image -->  \nWhile our result (12) is systematically larger than that derived by GPZM17 (4) corresponding to h ( e 012 ) = 1, it is partially explained by the difference of the break definition ( § 2) and our longer integration time (10 8 P 12 instead of 100 P 012 in GPZM17). We also note that the stability boundary at T break ∼ 100 P 012 ≲ 100 P 12 is fairly independent of the value of e 012 . Therefore, we conclude that our result is consistent with the GPZM17 result.  \n%  \nFigure 5. Same as Figure 4, but the vertical axis is now normalized by the GPZM17 result. The red dashed lines are empirical fits, h ( e 012 ), introduced here to include e 012 dependence.  \n<!-- image -->  \nIncidentally, Figure 5 suggests that h ( e 012 ) has an additional i mut -dependence. If we replace 1 . 5 by 1 . 6 for i mut = 90 · and by 1 . 8 for i mut = 180 · , for instance, the fit matches better our simulation result. Nevertheless, it does not lead to a significant difference, and we prefer the simpler empirical fit here. If we consider a simple linear interpolation, for instance, the above additional i mut -dependence might be expressed as  \n˜ h ( e 012 ) = 1 . 5 + 0 . 15 i mut 90 · +0 . 6 e 012 , (13)  \nalthough systematic simulations are required for further extension. We also note that the slope of h ( e 012 ), 0 . 6, coincides with that of MA01 criterion for HT-P (see equation (1)) when e 012 ≪ 1.  \nThe break timescale distribution in Figure 4 (or Figure 5) also presents interesting features. For circular outer orbits ( e 012 = 0), especially for coplanar systems, a break occurs suddenly with a very short timescale (typically ≲ 10 P 12 ), even around the stability boundary. In contrast, highly eccentric systems ( e 012 > 0 . 4) and polar systems ( i mut = 90 · , in the strong ZKL oscillation regime) show gradually increasing break timescales towards the stability boundary at 10 8 P 12 . This behavior is similar to the result in Hayashi et al. (2023) for HT-P systems, and indicates the importance to define the stability boundary as a function of the timescale for HT-S configurations as well.  \nSo far, we defined the stability boundary at t = 10 8 P 12 . As indicated in Figures 4 and 5, however, such boundaries defined at different epochs behave differently and exhibit different dependence on e 012 and i mut , in particular. While we do not attempt in the present paper, it is feasible to obtain the improved fit by generalizing the function h ( e 012 ) to h ( e 012 , i mut ; t ) from Figures 4 and 5. Indeed, this is also the case for the stability boundary for HT-P systems shown in Figures 8 - 10 of Hayashi et al. (2022).', '3.3. Mass ratio dependence': 'Figure 6 shows how the stability boundary depends on the mass ratios, m 1 /m 0 and m 2 /m 1 , in the cases of near-coplanar prograde systems ( i mut = 5 · ) with fixed initial phases. While our simulations adopt m 0 = 10 6 M ⊙ and P 12 = 1yr for definiteness, the result are scalable with respect to those values; see equation (14) in Hayashi et al. (2022). Upper-left, upper-right, lower-left and lower-right panels correspond to ( m 1 , m 2 ) = (10 M ⊙ , 10 M ⊙ ), (19 M ⊙ , 1 M ⊙ ), (100 M ⊙ , 100 M ⊙ ), and (1 M ⊙ , 1 M ⊙ ), respectively. The vertical axis of Figure 6 is normalized by the GPZM17 result as we did in Figure 5.  \nFour panels in Figure 6 appear to be almost identical, indicating that the mass dependence of the stability boundary is described by the conventional Hill-radius scaling ( m 0 / ( m 1 + m 2 )) 1 / 3 alone, independent of m 2 /m 1 .  \nFinally we also check the mass scalability, expected from the Newtonian gravity, with different m 0 values fixing mass ratios m 1 /m 0 and m 2 /m 1 , and confirm it statistically although chaotic behavior may change individual simulation result (see Hayashi et al. 2022).', '4.1. Classification of the fate of HT-S systems': 'Now we are in a position to examine the fate of HT-S systems after the initial binary of m 1 and m 2 breaks up. As illustrated in Figure 2, those systems should form a massive inner binary orbited by a less massive tertiary m out , which is similar to the HT-P configuration. While the stability of such HT-P systems has been discussed for instance by Mardling & Aarseth (1999, 2001); Hayashi et al. (2022, 2023), we are interested in black-hole triple systems in particular. Even though we do not take into account general relativity in the present analysis, we consider approximately the merger/collision of the m i -m j pair based on their Schwarzschild radii r s ( m ) = 2 Gm/c 2 ; see (B2) and (B3) in Figure 2.  \nTo be specific, we adopt the the merger/collision criterion (see e.g. Tanikawa & Umemura 2014):  \nr ij < 10[ r s ( m i ) + r s ( m j )] , (14)  \nwhere r ij is the distance between m i and m j . This condition is not accurate and should be interpreted as qualitative. Nevertheless, it provides an insight in what circumstances general relativity plays an important role, which will be useful for future work including general relativistic effects in a postNewtonian fashion that we plan to study separately. We again note that all the result here is scalable under pure Newtonian gravity, except for this merger/collision criterion.  \nAs in Figure 2, we consider the four fates of the systems:  \n(B1) ejection of m out : This is the most common fate in conventional HT-P systems. We adopt a specific definition of the ejection following Hayashi et al. (2022), but using the semi-major axes rather than orbital energies this time (see § 2).', '𝑖 )*+ = 5 ∘ , 𝑚 -= 10 . 𝑀 ⊙ , 𝑃 01 = 1yr': 'Figure 6. Break-up time distribution on (˜ r 012 / ˜ r 12 ) norm -e 012 plane. The vertical axis is normalized by the GPZM17 stability boundary as in Figure 5. Four panels correspond to different sets of m 1 and m 2 , although fixing m 0 as 10 6 M ⊙ .  \n<!-- image -->  \n- (B2) collision of m 1 and m 2 : This fate appears very unlikely after the initial m 1 -m 2 binary breaks up. As will be shown below, such collisions occur fairly often in the coplanar cases ( i mut = 0 · and 180 · ). This branch is strongly suppressed when small mutual inclinations between the inner and outer orbits are introduced. This implies that the initial condition of the binary is imprinted somewhere; the system comes back again to a configuration close to the tightly bound m 1 -m 2 binary even after chaotic dynamical evolution, especially when the dynamics is confined in a lower-dimension parameter space.\n- (B3) m in -m 0 merger: This fate is expected to be fairly common because the HT-S systems first experience the Hill instability and form m in -m 0 binary, at least temporarily.\n- (B4) survived/undisrupted: We define survived systems that are undisrupted until the end of our integration time limit t int = 10 8 P 12 . Needless to say, this definition is quite dependent on t int that is determined due to the computational cost in practice. Nevertheless, we emphasize that our adopted value is significantly longer than the previous ones, which can identify the fraction of systems in (B1), (B2) and (B3) that may have been overlooked before.', '4.2. Branching ratio of fates of HT-S systems': "As will be shown below, the fate of HT-S systems after the initial binary break-up is determined over a much longer timescale than T break . Therefore, it is computationally demanding to survey the multi-dimensional parameter space faithfully. Instead, we perform a set of N-body simulations with the following specific initial conditions.  \nAlthough we have systematically surveyed the parameter Υ in searching for the stability boundary (6), we fix Υ = 0 . 1Υ crit to check the evolution after the initial binary breaks. The outer eccentricity is varied as e 012 = 0 . 1 , 0 . 2 , 0 . 3 , · · · , 0 . 9, and we consider four different mutual inclinations; i mut = 0 · , 5 · , 90 · and 180 · . For each model, we generate 100 realizations with random initial phases ( ω 12 , ω 012 , M 12 , M 012 ). For simplicity, we fix e 12 = 0 . 5 assuming that e 12 does not significantly change the fate when the initial phases are statistically randomized. Again, we specifically adopt m 0 = 10 6 M ⊙ and m 1 = m 2 = 10 M ⊙ .  \nSuch systems are very Hill-unstable, and the initial m 1 -m 2 binary breaks up typically within 10 P 12 . We continue to follow their orbital evolution up to t int = 10 8 P 12 , and obtain the statistics of the timescales when those systems reach the three different fates (B1), (B2) and (B3). The remaining undisrupted systems at t int are classified as (B4).  \nFigures 7 to 10 plot the resulting timescales of the initial binary break (red), and the different fates (black) for i mut = 0 · , 5 · , 90 · and 180 · , respectively. In each figure, upper-left, upper-right, lower-left, and lower-right panels correspond to the fates (B1), (B2), (B3), and (B4), respectively.  \nThose figures imply that the binary break timescales are much less than 10 P 12 , although increasing with e 012 . Instead of log 10 T break /P 12 , if we plot log 10 T break /P 012 , it is almost independent of e 012 , implying T break is basically determined by P 012 . On the other hand, the subsequent disruption events ( m out -ejection, m 1 -m 2 collision, m in -m 0 merger) occur with significantly longer timescales and are sensitive to e 012 .  \nIncidentally, We found one exceptional case in which the original binary collides without experiencing the break; a black asterisk labeled 'unbreak' in Figure 7. Although this occurrence is quite rare under pure Newtonian gravity, the fraction may significantly increase including gravitational wave emissions, depending on the initial parameters. Therefore, this kind of fate may be important astrophysically in practice.  \nThe strong dependence of the fates on the outer eccentricity e 012 is clearly shown in Figure 11, which plots the fraction of each fate against e 012 ; i mut = 0 · (upper-left), 5 · (upper-right), 90 · (lower-left), and 180 · (lower-right).  \nExcept i mut = 0 · , the dependence of the fate fraction on e 012 is very similar for the other mutual inclinations. If the initial outer orbits are significantly eccentric ( e 012 > 0 . 7), m in -m 0 merger (B3) dominates. For near-circular outer orbits ( e 012 < 0 . 4), those systems are mostly undisrupted (B4) with small fraction experiencing the m 1 -m 2 collision (B2). For the intermediate value of the outer eccentricity (0 . 4 < e 012 < 0 . 7), the ejection of m out (B1) becomes sub-dominant. The trends of the corresponding timescales for those fates are consistent with the above picture.  \nFigure 7. The distributions for timescales of binary break (red) and each outcome (black) for i mut = 0 · . Upper-left, upper-right, lower-left, lower right panels correspond to the results for (B1) m out -ejection, (B2) m 1 -m 2 collision, (B3) m in -m 0 merger, and (B4) survival at 10 8 P 12 , respectively.  \n<!-- image -->  \n∘  \nFigure 8. Same as Figure 7, but for i mut = 5 · .  \n<!-- image -->  \n𝑖  \nFigure 10. Same as Figure 7, but for i mut = 180 · .  \n<!-- image -->  \n= 90  \n∘  \nFigure 9. Same as Figure 7, but for i mut = 90 · .  \n𝑖  \n= 180  \n<!-- image -->  \n∘  \nFigure 11. Fraction of each outcome for i mut = 0 · (upper-left), 5 · (upper-right), 90 · (lower-left), and 180 · (lower-right), respectively.  \n<!-- image -->  \nThe case with exactly coplanar prograde orbits is exceptional. Our interpretation is that the initial condition of the tightly bound m 1 -m 2 binary is somehow imprinted in the phase space of the three bodies. After long-term dynamical evolution, the system sometimes reaches the configuration close to the initial condition, significantly increasing the probability of the m 1 -m 2 collision that is unlikely in more generic initial conditions. In particular, for coplanar systems, the trajectories of three bodies are bounded on a two-dimensional orbital plane, and their physical degree-of-freedom is expected to be more limited by that of inclined triples. Indeed, Ford et al. (2001) consider nearly coplanar star-two giant planet systems, where two planets initially have near orbits, and show that a fairly significant fraction of systems experience two planet collisions. We suggest that similar orbital configurations can be realized from coplanar HT-S systems after breaking an initial binary.  \nNevertheless, the difference between coplanar prograde and retrograde systems is hardly understood from this interpretation. One possible physical interpretation may be the difference of energy and/or angular momentum exchanges affect the evolution (see Figure 14 and 15 in Hayashi et al. (2022), for example), and changes the initial condition after breaks. This may point to an interesting mathematical problem, but is beyond the scope of the present paper.", '5. SUMMARY AND CONCLUSION': "We have studied the long-term dynamical evolution of hierarchical triple systems comprising a central massive body and a tight binary in eccentric orbits, which we referred to as HT-S systems. The two major purposes are to find the Hill-type stability condition for the binary due to the tidal effect of the central body, and to classify the fate of the systems after the binary break-up as a function  \nFigure 12. Stability boundaries for HT-P (red curves, e in = 0, m 12 /m 3 , m 2 /m 1 = 20 / 2, 10 / 10, respectively.), and HT-S (blue curves, e 12 = 0, m 12 /m 0 , m 2 /m 1 = 2000 / 10 6 , 1000 / 1000, respectively.) in e out -a out (1 -e out ) /a in plane. We use equation (16) in Hayashi et al. (2022) for HT-P stability boundaries, and equation (6) in this paper for HT-S stability boundaries, respectively.  \n<!-- image -->  \neout  \nof the initial outer eccentricity e 012 . Using a series of direct N-body simulations, we obtained the stability criterion for the binary, inequality (6), which generalizes the formula (4) derived by Grishin et al. (2017).  \nWe compare the binary stability boundary (6) for HT-S systems and the disruption criterion for HT-P systems, equation (16) in Hayashi et al. (2022) in Figure 12. Those systems above the boundary lines are stable up to t int = 10 8 P 12 (HT-S; this work) and t int = 10 9 P in (HT-P; Hayashi et al. 2022). Note that the vertical axis of the figure is a out (1 -e out ) /a in so as to plot both HT-S and HT-P systems simultaneously. Thus, the scaling with respect to e 12 , equation (8), is not incorporated in Figure 12.  \nOur major findings are summarized as follows.  \n- (A) Outer eccentricity: The most important parameter that determines the stability of hierarchical triple systems is the ratio of the semi-major axes of the outer and inner orbits, a 012 /a 12 . In general, eccentricities of both orbits, e 012 and e 12 , tend to destabilize the systems. The destabilizing effect of e 012 is mostly incorporated by replacing a 012 with a 012 (1 -e 012 ), but the stability boundaries after the scaling still weakly increase as a function of e 012 . We find that the residual effect for HT-S systems is well approximated by a linear function of e 012 , equation (12), for e 012 < 0 . 8, almost independent of the mutual inclination i mut . This is consistent with  \nthe behavior of HT-P systems. The latter systems exhibit further stronger e out -dependence for e out > 0 . 9, but we did not check the regime for HT-S systems.  \n- (B) Inner eccentricity: The effect of e 12 on the stability is weaker but more subtle than that of e 012 . We find empirical scaling expressions for the e 12 -dependence that is sensitive to i mut .\n- (C) Mutual inclination of the inner and outer orbits: Our simulations with i mut = 0 · , 5 · , 90 · and 180 · confirm the i mut -scaling of the binary stability criterion for HT-S systems that was discovered before by Grishin et al. (2017). After applying this scaling, we were able to separate and find the above scaling relations for e 12 and e 012 .\n- (D) Mass ratio: The mass dependence of the binary stability criterion is well described by the simply factor of ( m 0 / ( m 1 + m 2 )) 1 / 3 , at least for near-coplanar prograde systems. It may be a bit surprising that the criterion is determined by the total mass of the initial binary m 1 + m 2 , and insensitive to their ratio m 1 /m 2 .\n- (E) Fate of HT-S and HT-P systems: Stability criteria plotted in Figure 12 exhibit some similarity between HT-S and HT-P systems, but they have very different implications for their fates. For HT-P systems, the criterion is time-dependent and corresponds to their disruption event practically (almost always, the ejection of the tertiary body from the system); see detailed discussion in Hayashi et al. (2022). In contrast, the criterion for HT-S systems that we derived in the present paper concerns the Hill-type binary break-up, and still the triple system itself is undisrupted. In a sense, the criterion should be interpreted to indicate the transition from HT-S to HT-P systems; see Figure 2. Such timescales for the binary break-up are typically very short ( ≲ 10 P 12 ), and the final fates of those systems are determined through their orbital evolution over much longer timescales (more than several order-of-magnitudes longer than the binary breaking timescale, similar to those for the stability criterion for the HT-P systems plotted in Figure 12). Those investigations were made possible by performing long-time numerical integration of the direct N-body systems for the first time. We find that the HT-S systems reach different fates (ejection, collision, and merger) and that their corresponding timescales are decreasing functions of e 012 in general. Also the fraction of each fate is sensitive to the initial value of e 012 ; see discussion in § 4.2. In long-timescale calculations, stellar evolution should play an important role for stellar triples in practice. For instance, Toonen et al. (2022) study the evolution and fates of stellar triples, including stellar evolution.  \nWhile we do not consider general relativity in the present work, the Hill stability including GR effects becomes important for some parameter ranges. For instance, Suzuki et al. (2020) derived a criterion using analytic treatments concerning Sundman's inequality and first-order post-Newtonian (1PN) approximate numerical simulations for coplanar near-circular systems. They suggested that the PN corrections tend to stabilize a system compared with a purely Newtonian case. It is also wellknown that the first-order GR correction (apsidal precession term) suppresses the ZKL oscillations if two timescales become comparable. Besides, gravitational wave emissions reduce the energy for a long timescale, and therefore shrink the semi-major axis of a binary. Close scatterings before collisions may also be significantly affected by the GR corrections. Therefore, the GR corrections affect the Hill stability of triple systems with relevant parameter values.  \nFinally, we would like to emphasize that the result in this paper is basically scalable satisfying so called Kepler's third law, under the assumption that a system is completely dominated by the Newtonian gravity. Therefore, the result is applicable, not only for black hole triples, but also for planetary (or satellite) systems.", 'ACKNOWLEDGMENTS': 'T.H. gratefully acknowledges Atsushi Taruya for fruitful discussions about dynamics in three-body systems. The numerical simulations were carried out on the local computer cluster awamori , and the general calculation server from Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan (NAOJ). T.H. gratefully acknowledges the fellowship by Japan Society for the Promotion of Science (JSPS). This work is supported partly by the JSPS KAKENHI grant No. JP23K25908 (Y.S.), JP21J11378 and JP23KJ1153 (T.H.), and JP21K13914 (A.A.T.).', 'A. EFFECT OF DIFFERENT INITIAL PHASES ON THE STABILITY BOUNDARY': 'So far, all the simulation runs have been performed for a given set of initial phases of three bodies. We vary the initial phases so as to test if our stability boundary is robust against those changes. For that purpose, we focus on coplanar prograde systems ( i mut = 0 · ) in circular ( e 12 = e 012 = 0) and highly eccentric ( e 12 = e 012 = 0 . 8) orbits, and vary the initial phases as follows.  \nFor circular systems, the pericenter arguments ( ω 12 , ω 012 ) are irrelevant, and we arbitrarily fix them as 0 · that simply define the zero point of the locations of each body for coplanar systems. In addition, due to symmetry, we set M 12 = 0 · without loss of generality. Thus, the initial phases are parameterized by the mean anomaly difference ∆ M ≡ M 012 -M 12 = M 012 alone, and we vary it as 2 · i ( i = 1 , 2 , · · · , 180).  \nFor eccentric systems, however, the initial phases are specified by three independent parameters; the relative location of pericenters (∆ ω ≡ ω 012 -ω 12 = ω 012 ), and mean anomalies ( M 12 and M 012 ). Therefore, we vary them as ∆ ω = ω 012 = 30 · j , M 12 = 60 · k , and M 012 = 60 · l , where the integers j , k , and l run from 1 to 6; see Figure 3.  \nFigure 13 plots how T break /P 12 for circular systems is affected by the mean anomaly difference ∆ M . There is a clear systematic trend in the stability boundary, up to ∼ 20%, as a function of ∆ M . This result may seem somewhat counter-intuitive, because the initial phases are expected to be well mixed up after the orbital evolution. However, the binary break in the circular systems happens so quickly, T break ≲ 10 P 12 , before the memory of the initial location is lost.  \nThe variation of the stability boundary for circular systems, however, turned out to be less than 20%, and thus not important in practice. While the origin of the above behavior of circular systems may be possibly related to the chaos theory (e.g. Mardling 1995a,b), it is beyond the scope of this paper and we do not consider further.  \nFigure 14 shows the results for eccentric systems. Each panel corresponds to different values of M 12 , and the horizontal axis denotes M 012 here. Similarly to the visualization treatment for e 12 in the previous sections, we plot the results for different ∆ ω by shifting data sequences horizontally. Contrary to the circular systems, Figure 14 shows no systematic dependence on initial phases for eccentric systems. This is because eccentric systems around the stability boundary are relatively  \n∘  \n(  \nFigure 13. Break-up time distribution in ∆ M - ˜ r 012 / ˜ r 12 plane.  \n<!-- image -->  \nmore stable and the longer timescale for the orbital evolution before the instability loses the memory of their initial phases. So, we conclude that the initial phase dependence is negligible for stability boundaries of the eccentric systems.  \n<!-- image -->  \n∘  \n∘  \n∘  \n∘  \n∘  \n∘  \nFigure 14. Break-up time distribution in M 012 - ˜ r 012 / ˜ r 12 plane.  \n<!-- image -->'}