Dynamical disruption timescales and chaotic behavior of hierarchical triple systems

Mardling (1995a, b); Mardling & Aarseth (1999, 2001) considered evolution of a tidally interacting binary. They defined that a binary evolution is chaotic when the direction and amount of energy flows between the tide and orbit become entirely unpredictable. In reality, the boundary between normal and chaotic behavior is not well defined, since it is extremely sensitive to small differences of orbital parameters and regions of stability and chaos are intertwined in a complicated fashion. Nevertheless, they were able to derive an approximate formula for the boundary.

Remarkably, MA01 recognized that the way of energy and angular momentum exchange between the inner and outer orbits in a hierarchical triple is very similar to that in a tidally interacting binary. On the basis of the analogy, MA01 proposed the following dynamical instability criterion for hierarchical triples:

In the above expression, $a_{\mathrm{in}}$ and $a_{\mathrm{out}}$ are the semi-major axes of the inner and outer orbits, $e_{\mathrm{out}}$ is the eccentricity of the outer orbit, $i_{\mathrm{mut}}$ is the mutual inclination between the two orbits, and $m_{12}=m_{1}+m_{2}$ is the total mass of the inner binary.

The numerical coefficient of 2.8 in inequality (1) is empirically introduced by MA01 from numerical simulations. They define the stability based on whether two orbits, differing by $1$ part in $10^{5}$ in the eccentricity, remain similar after 100 outer orbits. The additional dependence on the mutual inclination, $(1-0.3i_{\mathrm{mut}}/\pi)$, is heuristically introduced in Aarseth & Mardling (2001) based on the result of stability enhancement due to $i_{\mathrm{mut}}$ discussed by Harrington (1972). The extension of this criterion is done by Mylläri et al. (2018) to include the dependence on the mass ratio of the inner binary $m_{2}/m_{1}$, the mutual inclination $i_{\mathrm{mut}}$, and the inner eccentricity $e_{\mathrm{in}}$.

Inequality (1) indicates that the ratio between the pericenter distance of the outer orbit and the semi-major axis of the inner orbit is the key parameter dictating the chaotic behavior of the triples. Triples with $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}<1$ are unstable because their outer and inner orbits cross each other. Inequality (1) should be interpreted as a criterion for the onset of chaotic evolution, rather than for disruption, but triples with $1<r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}<(r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}})_{\mathrm{MA}}$ begin chaotic orbital evolution and are expected to be disrupted eventually. Therefore, while inequality (1) is widely known as a standard criterion for dynamically unstable triples (e.g. Perpinyà-Vallès et al., 2019; Gupta et al., 2020; Tanikawa et al., 2020), it cannot be used to predict the disruption time.

In order to estimate the disruption time of triples, MK20 apply the random walk (RW) approximation for the energy transfer between the inner and outer orbits. They focus on very eccentric outer orbits ($e_{\mathrm{out}}\sim 1$), and use an approximate analytic formula of the energy during the pericenter passages of a parabolic orbit by Roy & Haddow (2003). The transferred energy depends on the inner mean anomaly during the next outer orbital pericenter. Their RW model assumes that the inner mean anomaly changes randomly between adjacent outer pericenter passages without net energy exchange on average. The resulting energy transfer process, therefore, is supposed to be well described by a random walk. Furthermore, they assume that the disruption likely occurs when the accumulated transferred energy variance becomes comparable to the square of the outer orbital energy, and derive the disruption time of triples.

In their full RW model, they perform secular integration numerically up to the quadrupole interaction so as to include the effect of evolving orbital parameters. They also consider a simplified version of the RW model, in which the complicated evolution of orbital elements are neglected. Their simplified RW model yields the following analytic expression for the disruption time $T_{\mathrm{d}}$:

where $P_{\mathrm{in}}$ is the initial orbital period of the inner binary, and $m_{123}\equiv m_{1}+m_{2}+m_{3}$ is the total mass of the system.

The factor $2$ in equation (2) is heuristically determined from their $N$-body simulations. They found that the RW model estimations can reasonably reproduce the disruption times of numerical simulations within one order-of-magnitude for high-$e_{\mathrm{out}}$ triples; they also perform simulations for $e_{\mathrm{out}}=0.1$, 0.3, 0.5, 0.7, and $0.9$, and confirm that their models agree well with the simulations down to $e_{\mathrm{out}}=0.7$ although the accuracy of equation (2) significantly decreases at lower eccentricities.

The energy transfer model adopted by MK20 is based on a binary-single scattering process for a parabolic encounter by Roy & Haddow (2003). Indeed, this model turns out to be important to understand the dependence of the disruption time on the mutual inclination, we briefly summarize the result of Roy & Haddow (2003) in this subsection.

They apply the perturbation theory and compute the energy exchange during an encounter in which the collision parameter is much larger than the semi-major axis of a binary, neglecting rapidly oscillating terms. If the eccentricity of outer orbit is close to unity, the energy transfer during a pericenter passage is well approximated as a binary-single scattering for a parabolic encounter. The energy transfer during one pericenter passage for a triple system is written as (Mushkin & Katz, 2020):

where $r_{\mathrm{p,\mathrm{out}}}$, $\phi$, and $\Omega$ are the pericenter distance of the outer orbit, the characteristic phase angle, and the relative longitude of the ascending nodes, respectively.

The characteristic phase $\phi$ is defined as

where $\omega$ and $M^{*}_{\mathrm{in}}$ are the pericentre argument, and the inner mean anomaly during the passage, respectively. In equation (3),

where

The function $F$ in equation (3) is given by

and the dependence on $e_{\mathrm{in}}$ and $i_{\mathrm{mut}}$ is explicitly included in the three coefficients:

In the above expressions, $f_{i}$ are written in terms of the Bessel function of the $i$-th order $J_{i}$:

In section 6, we show that the energy transfer is suppressed for retrograde orbits using these expressions, and discuss a possible origin of the dependence of orbital configurations on the disruption timescales.

As illustrated in Figure 1, there are a number of parameters characterizing triples, and it is not realistic to cover the parameter space systematically. Thus, we adopt a fiducial set of parameters, which is listed in Table 1. The specific choice of $m_{12}$ is partly motivated by our previous papers (Hayashi et al., 2020; Hayashi & Suto, 2020) that consider a triple comprising an inner BBH and a stellar tertiary. Due to the scaling in the basic equations, however, the result can be easily applied for other cases. More specifically, mass and length, and time scales in purely gravitational problem admit the following scaling relation:

The relation (14) implies that the disruption time $T_{\mathrm{d}}$ computed from our set of simulation runs can be scaled as $T^{\prime}_{\mathrm{d}}=\sqrt{\beta^{3}/\alpha}~{}T_{\mathrm{d}}$ Indeed, this is why we mostly express our results below in terms of dimensionless variables like $T_{\mathrm{d}}/P_{\mathrm{in}}$.

Therefore, more relevant parameters are the mass ratio of the inner binary $q_{21}=m_{2}/m_{1}$, the mass ratio of secondary and tertiary $q_{23}=m_{2}/m_{3}$, the outer pericenter distance over inner semi-major axis $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$, the outer eccentricity $e_{\mathrm{out}}$ and the mutual inclination $i_{\mathrm{mut}}$.

In addition to those parameters, the variation of initial phase angles ($\omega_{\mathrm{in}}$, $\omega_{\mathrm{out}}$, $M_{\mathrm{in}}$, and $M_{\mathrm{out}}$) would introduce statistical variance on $T_{\mathrm{d}}/P_{\mathrm{in}}$. As will be discussed in appendix, we find that the resulting stochastic variance amounts to one or two orders of magnitudes, which is indeed comparable to the level of fluctuations due to the intrinsic chaotic nature of the triple systems that we consider in this paper.

We define the primary as the most massive body in the inner binary (so that $m_{1}\geq m_{2}$), and thus $0<q_{21}\leq 1$. While $q_{23}>1$ for a hierarchical triple consisting of inner massive binary and a tertiary, we consider cases of $q_{23}<1$ as well that are relevant for planetary systems in general. Figure 2 illustrates schematic pictures of triples in different $q_{21}$ and $q_{23}$ regimes.

In order to obtain the disruption time distributions of triples, we perform a series of numerical simulations with the $N$-body integrator TSUNAMI (Trani et al., 2019, Trani et al., in prep.). TSUNAMI is a direct $N$-body integrator specifically designed to accurately simulate few-body systems. The code solves the Newtonian equations of motion derived from a time-transformed, extended Hamiltonian (Mikkola & Tanikawa, 2013a, b). This numerical procedure, also called algorithmic regularization, serves to avoid the well known issue of gravitational integrators when two particles get very close together. When no regularization is employed, the acceleration grows quickly, and the timestep becomes extremely small, sometimes even halting the integration. Specifically, we employ the logarithmic Hamiltonian regularization scheme, which corresponds to setting $\alpha=1$, $\beta=0$, and $\gamma=0$ in the notation of Aarseth, Tout, & Mardling (2008). Furthermore, TSUNAMI is suited to evolve hierarchical systems of particles like the ones examined here, thanks to its chain-coordinate system that reduces round-off errors when calculating distances between close particles far from the center of mass of the system (Mikkola & Aarseth, 1993). Finally, TSUNAMI uses Bulirsch-Stoer extrapolation to improve the accuracy of the integration, ensuring accuracy and adaptability over a wide range of dynamical scales (Stoer & Bulirsch, 2002). As a result, TSUNAMI can perform accurate three-body simulations roughly $10-100$ times faster than other integrators. For more details on the integration scheme, we refer to Trani & Spera (2022).

The definition of the dynamical stability is not unique. Indeed, it may depend on the maximum integration time that we adopt in numerical simulations to some extent. First, we define the disruption time $T_{\mathrm{d}}$ for each triple following Manwadkar et al. (2020, 2021). At each timestep, the integrator evaluates the binding energy for each of the three pairs of bodies, i.e. $(m_{1},m_{2})$, $(m_{1},m_{3})$ and $(m_{2},m_{3})$. The pair with the highest (negative) binding energy is considered as the inner binary, and we call the pair of the inner binary and the remaining tertiary as the outer pair. When the binding energy of the outer pair becomes positive and the radial velocity of the tertiary becomes positive (i.e. is moving away from the inner binary), we tentatively regard it as a candidate of a disrupted triple, and assign the epoch as the disruption time $T_{\mathrm{d}}$. Some of them, however, represent a transient unbound state, and subsequently become gravitationally bound again. Therefore, we continue the run until its binary-single distance exceeds 20 times the binary semi-major axis so as to make sure that it is a truly disrupted triple. In that case, we stop the run, and the triple is classified as disrupted with the previously assigned value of $T_{\mathrm{d}}$. If a system becomes bound again, on the other hand, we reset the temporarily assigned value of $T_{\mathrm{d}}$, and keep running the integration of the system.

Finally, we define stable triples as those that do not disrupt and survive until our maximum integration time $t_{\mathrm{int}}$. We set $t_{\mathrm{int}}=10^{9}P_{\mathrm{in}}$ unless otherwise specified.

Due to the chaotic nature of the system, strict numerical convergence cannot be achieved, and the disruption timescales obtained from our simulation runs should be interpreted in a statistical manner. In §4.2 and appendix A, we attempt a numerical convergence test by running the same set of runs on different machines allowing for a tiny difference of the initial orbital parameter, as well as varying the initial phases of the inner and outer orbits and of the three bodies. In all cases, the disruption timescales agree within one or two orders of magnitudes, indicating that our results are statistically robust within those uncertainties.

This subsection presents our simulation results of the disruption timescales, focusing on triples with an equal-mass inner binary. More general cases will be described in section 5.

Figure 3 shows the normalized disruption time $T_{\mathrm{d}}/P_{\mathrm{in}}$ distribution on $e_{\mathrm{out}}$ - $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ plane for coplanar prograde triples with an equal-mass inner binary ($i_{\mathrm{mut}}=0^{\circ},q_{21}=1$). The secondary to tertiary mass ratio $q_{23}$ is fixed as $5.0$, i.e., the tertiary mass $m_{3}$ is $2M_{\odot}$ for $m_{1}=m_{2}=10M_{\odot}$. We adopt the fiducial values for the other parameters that are listed in Table 1. As discussed in section 4.2, $T_{\mathrm{d}}$ does not exactly scales with $P_{\mathrm{in}}$ because of the scatters inevitably induced by the chaoticity. Nevertheless, $T_{\mathrm{d}}/P_{\mathrm{in}}$ is independent of $P_{\mathrm{in}}$ statistically, and we consider the normalized disruption time $T_{\mathrm{d}}/P_{\mathrm{in}}$ in what follows.

Triples located above the magenta curves in Figure 3 are expected to be stable according to the MA01 dynamical stability condition. The disruption timescales predicted by the MK20 model are plotted in black contours (for $T_{\mathrm{d}}/P_{\mathrm{in}}=10^{3},10^{4},10^{5},10^{6},10^{7},10^{8}$ and $10^{9}$). The contours are nearly straight lines on $e_{\mathrm{out}}$ - $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ plane because the exponential term in equation (2) dominates the parameter dependence. The green shaded region in Figure 3 corresponds to the orbit-crossing triples whose outer orbit intersects the inner binary orbit at the pericentre, $a_{\mathrm{out}}(1-e_{\mathrm{out}})<a_{\mathrm{in}}$. Those triples are supposed to be very unstable in general, and are not our main interest in the present paper.

The result of simulated disruption timescales are plotted as filled circles with different colors in Figure 3 according to the side color-scale. We start the simulation with the initial parameters of $(r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}},e_{\mathrm{out}})=(1.0,0.04+0.04m)$ with the integer $m$ running over $0-23$. For a given value of $m$, we sequentially increase the value of $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ by 0.2. Due to the limitation of the CPU time, we stop each simulation at $t_{\mathrm{int}}=1.0\times 10^{9}P_{\mathrm{in}}$ even if the system is not disrupted, which is approximately equal to $2.7$ Gyrs for our fiducial parameters. For any given $m$, we stop increasing the value of $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ when the previous two realizations at lower $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ do not break before the final integration time. In other words, we scan the $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ - $e_{\mathrm{out}}$ of Figure 3 going from bottom to top, stopping only when two realizations in a row fail to break up before $t_{\mathrm{int}}$. Thus those triples located above the circle-cross symbols in Figure 3 correspond to $T_{\mathrm{d}}>1.0\times 10^{9}P_{\mathrm{in}}$. In what follows, we sometimes refer to the boundary as the simulated stability boundary just for convenience.

Figure 3 shows that the MA01 dynamical stability condition, inequality (1), roughly explains the stability boundary of the simulated coplanar prograde triples with an equal-mass inner binary; all triples with $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}>({r_{\mathrm{p,\mathrm{out}}}}/{a_{\mathrm{in}}})_{\mathrm{MA}}$ are stable at least until $t_{\mathrm{int}}$.

The predicted disruption timescales based on the MK20 RW model agree well with the simulation result only for triples with $e_{\mathrm{out}}\gtrsim 0.8$. This is indeed expected because their model assumes that the energy transfer between the inner and outer orbits is approximated by a random walk for a highly eccentric outer orbit. MK20 also confirmed that their RW model reproduces their simulations reasonably well down to $e_{\mathrm{out}}\sim 0.7$, which is consistent with our present result. The right panel in Figure 3, furthermore, implies that the lower limit of $e_{\mathrm{out}}$ where the RW model approximation is valid becomes larger as $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ increases.

In reality, however, the above rough agreement between the theoretical predictions and simulations does not hold for orthogonal ($i_{\mathrm{mut}}=90^{\circ}$) and coplanar retrograde ($i_{\mathrm{mut}}=180^{\circ}$) cases; see Figures 4 and 5.

Orthogonal triples are less stable than the MA01 criterion, which may be due to their drastic orbital evolution. In particular, the effect of the ZKL oscillation is not properly taken into account in the MA01 model. In contrast, the MK20 RW model seems to agree better than the prograde simulations for even lower values of $e_{\mathrm{out}}$. Finally, Figure 5 clearly indicates that coplanar retrograde triples are significantly more stable than MA01 and MK20 predict. We will discuss this point later based on the energy transfer model by Roy & Haddow (2003).

The results in Figures 3 to 5 are computed initially for a set of fixed phase parameters including mean anomalies ($M_{\mathrm{in}}$, $M_{\mathrm{out}}$) and arguments of pericenter ($\omega_{\mathrm{in}}$, $\omega_{\mathrm{out}}$). While all those initial parameters are mixed up in the course of the dynamical evolution, their different values at the initial epoch are likely to affect the disruption timescales to some extent, and to introduce statistical scatters in the estimated timescales. This is quite expected, and we show in appendix A that those phase angles change the disruption timescales by one or two orders of magnitudes.

Instead, this subsection addresses the question how the disruption timescale is sensitive to the (tiny) difference of the initial orbital parameters, given the chaotic nature of the gravitational many-body systems in general (e.g., Suto, 1991). For this purpose, Figure 6 examines quantitatively to what extent the disruption timescale is affected by the tiny difference of the initial period of the inner binary.

We select three sets of runs from our fiducial model simulation, and change the initial value of $P_{\mathrm{in}}$ (and therefore that of $a_{\mathrm{in}}$) alone while keeping the other parameters intact. In practice, we change the initial input value of $P_{\mathrm{in}}=10^{3}$ days only by its 15-th digit (left panels) and 4-th digit (right panels), which are the horizontal axes of Figure 6, i.e., $n=10^{15}\Delta P_{\mathrm{in}}/P_{\mathrm{in}}$ and $n=10^{4}\Delta P_{\mathrm{in}}/P_{\mathrm{in}}$. We also repeat the identical run on three different machines; AMD Ryzen7 Core in CfCA calculation cluster at NAOJ, Intel Core i9 on 16-inch Macbook Pro 2019, and AMD EPYC 7502 in our local computer cluster. In this subsection, we stop the runs at $t_{\mathrm{int}}=4.0\times 10^{7}P_{\mathrm{in}}$ to save computational time.

Let us look at the left panels of Figure 6 first. The bottom panel corresponds to a case that the disruption time is insensitive to such a small change. The cases plotted in the top and middle panels, in contrast, indicate that $T_{\rm d}$ varies significantly among different machines; tiny differences in $P_{\mathrm{in}}$ lead to even by a few orders of magnitudes. This amazing sensitivity to the initial condition originate from the intrinsically chaotic nature of the gravitational dynamics of triples. The variability among different machines is instead due to the limits of floating-point arithmetic and the differences in optimizing compilers that result in slightly different machine code being generated. Together with the extreme sensitivity to initial conditions, these factors cause the introduction of “numerical chaos”, when small but different accumulations of integration errors result in a completely different evolution of the same system.

We repeated those runs except that we modify only the 4-th digit of the initial $P_{\mathrm{in}}$. The result is shown in the right panels of Figure 6, which are very similar to the left panels. For a triple whose $T_{\rm d}$ is insensitive to the 15-th digit of the input $P_{\mathrm{in}}$ (the bottom panels), the value is identical in the three different machines even if the 4-th digit is modified, implying that the disruption of the triple seems to proceed in a regular, instead of chaotic, fashion. We could speculate that the Lyapunov time of such system is so long that even a relatively large change in the initial conditions does not result in a different evolution. In contrast, the range of variance of $T_{\rm d}$ is statistically the same between the left and right panels for the top and middle examples, suggesting that their chaotic nature is fairly independent of the fluctuation amplitude of the input $P_{\mathrm{in}}$. Such a chaotic degree of the unstable triples is extremely sensitive to the set of orbital parameters in a complicated fashion, and is difficult to predict. Nevertheless, Figure 6 indicates that the disruption timescales from our current simulations should be understood to vary by one or two orders of magnitudes, in practice.

The disruption timescale $T_{\mathrm{d}}$ should admit the scaling relation, equation (14), with respect to $P_{\mathrm{in}}$ and $(m_{1},m_{2},m_{3})$, as long as non-Newtonian effects are neglected. For instance, the RW model prediction, equation (2), respects the scaling. In reality, however, the numerical results do not necessarily satisfy the scaling due to the intrinsic chaoticity described in the above.

In order to see this, we run simulations for a subset of fiducial triple systems corresponding to Figures 3 to 5 except that we change the value of $P_{\mathrm{in}}$ from $10^{3}$ days to $10^{4}$ days (note that $P_{\mathrm{out}}$ is derived from the other fixed parameters). We define the ratio of the normalized disruption timescales for individual triples of $P_{\mathrm{in}}=10^{3}$ and $10^{4}$ days:

Figure 7 plots $R$ against $\log_{10}(T_{\mathrm{d}}/P_{\mathrm{in}})_{P_{\mathrm{in}}=1000{\rm days}}$ for prograde (top), orthogonal (middle), and retrograde (bottom) triples with $q_{21}=1$ and $q_{23}=5$. In order to save the CPU time, we adopt $t_{\mathrm{int}}=4\times 10^{7}P_{\mathrm{in}}$ in those examples.

The ratios derived from simulations are not exactly equal to unity due to the chaotic behavior of the triple dynamics, but they distribute around unity in a statistical sense. For instance, we find that their distribution is well approximated by log-normal functions as indicated in the right panels.

In summary, we find that the numerically estimated disruption timescales are inevitably affected by the chaotic behavior of triple dynamics. The timescales may vary one or two orders of magnitudes due to the differences of the initial phase angles and also to the tiny changes in the other orbital parameters. Nevertheless, we conclude that the disruption timescales are statistically robust, and can be predicted from the orbital parameters within one or two order-of-magnitude uncertainty.

Figures 3 to 5 suggest that there is a fairly sharp boundary of the triple stability on $e_{\mathrm{out}}-r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}$ plane. The boundary seems to be roughly approximated by shifting inequality (1) along the vertical direction. Thus, we empirically introduce the coefficient $K$ to the MA01 model, and fit the simulated stability boundary as

The results are shown in Figures 8 to 10, together with the value of $K$.

Although equation (16) is a simple extrapolation of the MA01 model, its $e_{\mathrm{out}}$ dependence reproduces the numerical results reasonably well by fitting the value of $K$, especially given the expected scatters discussed in subsection 4.2. The boundaries between stable and unstable regions appear to be very sharp especially for $e_{\mathrm{out}}<0.3$ and becomes gentle for $e_{\mathrm{out}}>0.8$, at least for $i_{\mathrm{mut}}=0^{\circ}$ and $180^{\circ}$ (Figures 8 and 10). Unlike inequality (1), however, the stability boundary does not monotonically depend on the mutual inclination $i_{\mathrm{mut}}$, and the orthogonal triples are more unstable than the other two coplanar configurations; see Figure 9. In all cases, the disruption timescales for those triples with ${r_{\mathrm{p,\mathrm{out}}}}/{a_{\mathrm{in}}}<({r_{\mathrm{p,\mathrm{out}}}}/{a_{\mathrm{in}}})_{\rm SB}$ are mainly determined by the value of ${r_{\mathrm{p,\mathrm{out}}}}/{a_{\mathrm{in}}}$, in agreement with the MK20 model.

Figure 11 plots the best-fit values of $K(i_{\mathrm{mut}},q_{21},q_{23})$, which clearly shows the parameter dependence of the stability boundaries in Figures 8 to 10. The fit is not so good in some cases, especially for orthogonal triples. Nevertheless, we perform the least-square fitting for them, and plot the best-fit values of $K$ in Figure 11, so as to indicate the approximate behavior of the stability boundary in a qualitative sense. The poor fit for orthogonal triples may imply the importance of their complicated orbital evolution due to the ZKL oscillation; see section 6. Relative to the MA01 criterion, orthogonal triples are more unstable and coplanar retrograde triples are more stable. The dependence on $q_{23}$ is not well represented by the combination of $(1+m_{3}/m_{12})$ alone. These differences may likely come from the fact that their model is based on the simulations up to $10^{2}P_{\mathrm{out}}$, which is much shorter than our $t_{\rm int}=10^{9}P_{\mathrm{in}}$. Of course, it simply reflects the progress of computational resources over the last two decades, and we would like to emphasize the amazing insight and quality of the pioneering work by MA01. We plan to explore the dependence of the dynamical stability boundary on $i_{\mathrm{mut}}$, $q_{21}$ and $q_{23}$ in detail and report the result elsewhere.

We find that there is a relatively clear boundary between dynamically stable and unstable regions on $r_{\mathrm{p,\mathrm{out}}}/a_{\mathrm{in}}-e_{\mathrm{out}}$ plane. The boundary is reasonably well approximated by equation (16), which is a simple generalization of the MA01 model. In the unstable regions, the disruption timescales seem to be roughly consistent with the MK20 model. Thus, following equation (2), we introduce a parameter $x$:

Then, equation (2) is rewritten in terms of $x$ as follows, which clearly implies that the value of $x$ basically determines the disruption time:

Figure 12 compares the $T_{\mathrm{d}}/P_{\mathrm{in}}$ distribution from our simulation against the MK20 model; top, middle, and bottom panels show the prograde ($i_{\mathrm{mut}}=0^{\circ}$), orthogonal ($i_{\mathrm{mut}}=90^{\circ}$), and retrograde cases ($i_{\mathrm{mut}}=180^{\circ}$), while left and right panels are for equal-mass ($q_{21}=1.0$) and unequal-mass inner binaries. The results with different $q_{23}$ are plotted in different colors: $q_{23}=5.0$ (magenta), $1.0$ (blue), and $0.5$ (black). Thin solid and dotted lines indicate the MK20 model predictions, equation (19), for $e_{\mathrm{out}}=0.04$ and $0.96$, respectively, so as to bracket the dependence on $e_{\mathrm{out}}$. Filled circles in different colors correspond to our simulation data in Figures 8 to 10. The sequence of data located at $T_{\mathrm{d}}/P_{\mathrm{in}}=10^{9}$ corresponds to triples that are dynamically stable up to our integration time $t_{\mathrm{int}}=10^{9}P_{\mathrm{in}}$, and thus are not considered in the discussion below.

Let us examine more quantitatively Figure 12. Coplanar prograde triples (top panels) do not exhibit noticeable dependence on $q_{23}$ in the range $0.5$ - $5.0$, which is consistent with the very weak dependence on $q_{23}$ that the MK20 model predicts. Furthermore, unequal-mass inner binaries ($q_{21}=0.1$) tend to be slightly more stable than equal-mass ones, which is in a qualitative agreement as well with the MK model.

Orthogonal triples (middle panels) exhibit no strong dependence on $q_{23}$ either, given the large scatters of simulated data, while unequal-mass inner binaries seem to be more stable for less massive tertiaries (with $q_{23}=5$) in contrast to the MK20 model. This may indicate the importance of other processes than the RW energy transfer, such as the ZKL oscillation up to octupole interaction. In section 6, we show one typical orbital evolution under the ZKL oscillation.

Coplanar retrograde triples (bottom panels), on the other hand, are significantly more stable than the above two cases, and behave very differently relative to the MK20 model, as already indicated in Figures 8 to 10. The fact that simply reversing the motion of the outer body changes significantly the disruption time distribution is not so intuitive, and we will discuss further in the next section.

Since we have not yet identified physical processes that explain the orbital configuration dependence on the disruption time, we introduce a simple empirical model, and determine the best-fit parameters to summarizes our current results. For that purpose, we consider the following parameterization:

where the dominant term in the MK20 RW model corresponds to $S=1$.

We first attempted to fit the results of Figure 12 by varying $S$ and $C$, and found that most of the results are fitted with $C$ close to unity. Given the scatters of the disruption timescales due to the initial phase angles and chaotic behavior of the triples, therefore, we decided to adopt $C=1$ and fit with $S$ alone. Thick-dashed lines in Figure 12 represent the best-fit models corresponding to equation (20) with $C=1$. Figure 13 shows the best-fit value and the associated 1$\sigma$-confidence for $S$. While equation (20) is a very rough approximation, it depicts our simulation results at least qualitatively, and is helpful to understand their basic behavior in a simple way.

Evolution of orbital elements, energies and angular momenta for hierarchical triples with $q_{21}=1$ and $q_{23}=5$ is plotted in Figure 14. Left, middle and right panels are for the coplanar prograde, orthogonal, and coplanar retrograde triples. In each panel, magenta and green curves represent the evolution of systems that disrupt around $10^{6}P_{\mathrm{in}}$ and remain stable at $4\times 10^{7}P_{\mathrm{in}}$, respectively.

From top to bottom panels, we plot the semi-major axis ratio of the inner and outer orbits $a_{\mathrm{in}}/a_{\mathrm{out}}$, the ratio of $(r_{\mathrm{p,out}}/a_{\mathrm{in}})$ and the MA01 threshold $(r_{\mathrm{p,out}}/a_{\mathrm{in}})_{\mathrm{MA}}$, the outer eccentricity $e_{\mathrm{out}}$, the inner eccentricity $e_{\mathrm{in}}$, the mutual inclination $i_{\mathrm{mut}}$, the fractional variation of the total energy and angular momentum of the outer orbit $|\Delta E_{\mathrm{out}}/E^{(0)}_{\mathrm{out}}|$ and $|\Delta L_{\mathrm{out}}/L^{(0)}_{\mathrm{out}}|$, $|W/E_{\mathrm{out}}|$, and finally $A_{i}$ ($i=1,2,$ equations 8–9) and $|A_{3}|$ (equation 10).

Both coplanar triples remain coplanar for the entire evolution, and their inner eccentricities are relatively small. In contrast, the initially orthogonal triple changes the mutual inclination between $40^{\circ}$ and $90^{\circ}$ due to the ZKL effect, which is accompanied by the large periodic oscillation of $e_{\mathrm{in}}$. As we will show in the next subsection, this is why the initially orthogonal triples are more unstable than the coplanar cases.

After a certain period of the energy transfer between the inner and outer orbits, $a_{\mathrm{in}}/a_{\mathrm{out}}$ drops and $e_{\mathrm{out}}$ is enhanced suddenly for unstable cases, leading to an ejection of the tertiary.

The fact that the disruption timescales are sensitive to the mutual inclination angle between the inner and outer orbits may be understood from the model by Roy & Haddow (2003), which is briefly summarized in section 2.3. Specifically, equations (3) and (7) indicate that the energy exchange between the inner orbit and a parabolic perturber is approximately estimated as

where $W$ is given by equation (5).

Let us assume that the energy transfer between the inner and outer orbits in our hierarchical triples is approximated by such an encounter, which happens periodically at a pericenter of the eccentric outer orbit.

The coefficients $A_{1}$ to $A_{3}$ are explicitly given by equations (8) to (13), and plotted in Figure 15 against $e_{\mathrm{in}}$ for $i_{\mathrm{mut}}=0^{\circ}$ (solid), $90^{\circ}$ (dotted), and $180^{\circ}$ (dashed). Expanding those equations with respect to $e_{\mathrm{in}}$, one obtains their leading-order expressions as

Equation (22) and three curves in Figure 15 seem to suggest that orthogonal triples are more stable than coplanar prograde triples if compared at the same value of $e_{\mathrm{in}}$. In reality, however, our simulation runs indicate that the orthogonal triples are the most unstable. This is due to the fact that the inner eccentricity significantly varies for orthogonal triples due to the ZKL oscillation, thus we have to take account of the evolution of $e_{\mathrm{in}}$ as well. From the conservation of the total angular momentum in the ZKL effect, one can derive the relation (e.g. Naoz, 2016):

for triples initially with $e_{\mathrm{in}}\approx 0$ and $i_{\mathrm{mut}}=90^{\circ}$. In the quadrupole approximation (or in the case of equal-mass inner binaries that do not have the octupole term), $e_{\mathrm{out}}$ and $a_{\mathrm{in}}/a_{\mathrm{out}}$ are constant, and equation (23) directly predicts $i_{\mathrm{mut}}$ as a function of $e_{\mathrm{in}}$. Dots in Figure 15 are sampled from the orbital evolution of the models shown in Figure 14, and are in good agreement with the prediction based on equation (23); orthogonal triples exhibit a wide range of $e_{\mathrm{in}}$, while in the other two coplanar cases, $i_{\mathrm{mut}}$ is constant and $e_{\mathrm{in}}$ stays smaller than 0.3 at most. Thus, the behavior of $A_{i}$ along the real orbital evolution of the three models nicely explains why the orthogonal triples are the most unstable.