Lagrange versus Lyapunov Stability of Hierarchical Triple Systems: Dependence on the Mutual Inclination between Inner and Outer Orbits

While there have been many studies examining the stability of hierarchical triple systems, the meaning of “stability” is somewhat vague and has been interpreted differently in previous literatures. The present paper focuses on “Lagrange stability,” which roughly refers to the stability against the escape of a body from the system, or “disruption” of the triple system, in contrast to “Lyapunov-like stability,” which is related to the chaotic nature of the system dynamics. We compute the evolution of triple systems using direct N-body simulations up to 107 P out, which is significantly longer than previous studies (with P out being the initial orbital period of the outer body). We obtain the resulting disruption timescale T d as a function of the triple orbital parameters with particular attention to the dependence on the mutual inclination between the inner and outer orbits, i mut. By doing so, we have clarified explicitly the difference between Lagrange and Lyapunov stabilities in astronomical triples. Furthermore, we find that the von Zeipel–Kozai–Lidov oscillations significantly destabilize inclined triples (roughly with 60° < i mut < 150°) relative to those with i mut = 0°. On the other hand, retrograde triples with i mut > 160° become strongly stabilized with much longer disruption timescales. We show the sensitivity of the normalized disruption timescale T d/P out to the orbital parameters of triple system. The resulting T d/P out distribution is practically more useful in a broad range of astronomical applications than the stability criterion based on the Lyapunov divergence.

Among all the literature, the dynamical stability criterion for hierarchical triples proposed by Mardling & Aarseth (2001, hereafter MA01) is well known and widely used in various studies of the triple dynamics. They examine the stability of triples as follows (R. Mardling 2022, private communication). First, they derive the analytic expression for the stability criterion based on the chaotic evolution boundary (see also Mardling 1995aMardling , 1995b. In order to determine the numerical factor, they integrate two almost identical systems (the given system and its "ghost") except for their different inner eccentricities e in of Δe in = 10 −7 , and monitor the difference in the inner semimajor axes a in at the time of outer apocenter passage. This variable is expected to be fairly insensitive to the initial difference for a stable system, while it is expected to grow exponentially for chaotic and unstable systems. This consideration motivated MA01 to judge and classify the stable and unstable systems from the departure of a in evaluated at 100 times the outer orbital period, P out (see also Mardling 2008).
The Lyapunov stability is mathematically defined as the stability of solutions of dynamical systems near their fixed (equilibrium) points, but has been also interpreted as the stability around arbitrary solutions (or trajectories) of dynamical systems against the tiny changes of system parameters (see, for instance, Lichtenberg & Lieberman 1983;Kandrup 1990;Suto 1991;Lichtenberg & Lieberman 1992). The methodology of MA01 is devised to identify the chaotic nature of triple systems in the spirit similar to the Lyapunov stability in the latter sense. We note that their methodology is close to the Lyapunov exponent method, which is widely used to characterize the chaoticity of systems. Therefore, we refer their stability condition as the Lyapunov stability below in a broader sense.
In any case, the Lyapunov stability or the local chaoticity may not be directly used in judging the global fate of the astronomical triple systems in reality. Even if triple systems are locally chaotic, their orbital configurations may be bounded, i.e., Lagrange stable. Indeed, the latter stability would correspond to a more relevant and practical definition in most applications for aircraft and nuclear plants (Gyftopoulos 1963), as well as for astronomical triple systems of our interest.
It seems that the stability criterion of triples derived by MA01 is sometimes misunderstood, and used as a criterion of triple disruption in the sense of Lagrange instability. Recently (Hayashi et al. 2022, hereafter Paper I) computed the disruption timescales T d of triples by directly integrating the three-body dynamics up to 10 9 P in , where P in is the initial orbital period of the inner binary. The distribution of the normalized disruption timescale T d /P in follows an approximate scaling similar to the MA01 criterion for coplanar-prograde triples in which the mutual inclination between inner and outer orbits i mut is 0°. The orthogonal and coplanar-retrograde triples, however, exhibit very different behavior, illustrating the clear difference between Lagrange and Lyapunov stabilities. In particular, Paper I finds that the von Zeipel- Kozai-Lidov (ZKL) oscillations (von Zeipel 1910;Kozai 1962;Lidov 1962) play a major role in disrupting significantly inclined triples. Because the characteristic ZKL timescale τ ZKL is comparable or even longer than 100 P out in general, their Lagrange stability cannot be determined from the short-term local behavior that characterizes Lyapunov stability.
The relation between the chaoticity and Lagrange stability is subtle. Figure 6 of Paper I illustrated such examples; the lower panels in the figure show that a Lagrange-unstable triple with T d ≈ 3000P in is robust against the change of the input initial conditions, unlike the other two cases (the upper and middle panels). Therefore, Paper I concluded that a Lagrange-unstable triple is not necessarily sensitive to the initial conditions. More strictly, a Lagrange-unstable triple may become chaotic eventually (just at the last moment of the disruption), but it is not clear how long it takes for them to exhibit noticeable chaoticity. Indeed, this is what we intend to examine quantitatively in Paper I and the present paper; we record the disruption timescales during the simulations, not merely the stable/unstable outcome alone.
In addition, we would like to emphasize that a chaotic triple does not always necessarily become Lagrange unstable. Thus, the Lagrange instability cannot be judged simply from the first few hundred (or fewer) orbits of a trajectory to its nearby ghost as MA01 assumed. This is exactly why we have performed a significantly longer time integration (up to 10 9 P in in Paper I and 10 7 P out in the present paper) to check the Lagrange instability of various triples. Only with such a direct confirmation of the disruption timescales can one understand the relation between the Lagrange instability and the chaotic behavior of the triple systems.
This point is also discussed in a different approach recently by Gajdoš & Vaňko (2023), who found that several observed exoplanetary systems have relatively short Lyapunov times (calculated from the Lyapunov exponent), even down to O(10) yr. This fact implies that the Lyapunov time is not necessarily related to the Lagrange stability of the systems, in good agreement with our finding mentioned above.
This paper examines the disruption timescale distribution for hierarchical triples using the same methodology of Paper I, and addresses the following questions more specifically: (i) the relation between Lyapunov stability criterion and Lagrange stability as a function of disruption timescales of triples with different orbital configurations, (ii) the dependence of disruption timescales on the mutual inclination, i mut , between the inner and outer orbits, and (iii) the effect of the relative orbital phases among the three bodies that has not been explicitly studied in most previous literatures.

Method
We consider hierarchical triple systems comprising an inner binary (m 1 and m 2 ) and a tertiary (m 3 ). Since one of the main purposes of the present paper is to identify the mutual inclination dependence on the disruption timescale, we simply fix the mass ratio as m 1 = m 2 = 5m 3 , corresponding to q 21 ≡ m 2 /m 1 = 1 and q 23 ≡ m 2 /m 3 = 5 in Paper I. The numerical simulation employs the N-body integrator TSUNAMI (Trani et al. 2019 that is specifically designed to accurately simulate few-body systems; further details of the code can be found in Trani et al. (2019), Trani & Spera (2022), and Paper I.
We consider an initially circular inner orbit with fixing the inner mean anomaly M in and the inner pericenter argument ω in . We specify their explicit values for each result below. In the barycentric reference frame whose x−y plane is chosen as the invariant plane perpendicular to its angular momentum vector, the inner and outer longitudes of the ascending node differ by 180°(see also Figure 1 of Paper I). Without a loss of generality, we set them to Ω in = 180°and Ω out = 0°. Therefore, the phase information of the initial conditions is then fully specified by the values of the outer pericenter argument and mean anomaly, ω out and M out . We fix these values first, and vary them randomly to check the initial phase dependence in Section 3.3.
We adopt the definition of the disruption time T d following Paper I and similar to Manwadkar et al. (2020Manwadkar et al. ( , 2021; at each time step, 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 consider 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 record the system as a candidate of the disrupted system tentatively. We found that the disruption normally occurs very fast after one body is ejected from a system (see also Figure 14 in Paper I).
A fraction of them, however, just represent a transient unbound state, and become gravitationally bound again. Thus, in order to make sure the system is truly disrupted, we continue the run until its binary−single distance exceeds 20 times the binary semimajor axis. In that case, we stop the run, and the triple is classified as disrupted with recording the epoch as T d . Otherwise, we reset the disruption candidate condition, and keep running the integration of the system.
Finally, we define a triple as (Lagrange) stable when T d > t int , i.e., when it does not disrupt and instead survives until our maximum integration time t int . Nevertheless, we would like to emphasize again that the above definition is just for convenience, and the main purpose of the present paper is to study the disruption timescales related to the Lagrange stability. In the present paper, we set t int = 10 7 P out unless otherwise specified, in contrast to Paper I that adopts t int = 10 9 P in , because the current paper expresses the normalized disruption time T d /P out instead of T d /P in . In reality, however, the result is identical except for the boundary of stable and unstable triples. We also note that under Newtonian gravity alone, the distribution of T d /P out is scale free (Paper I). So, if we consider an outer orbit of P out = 10 3 yr, for instance, the stopping time of our simulations corresponds to t int = 10 Gyr, i.e., the present age of the universe. Even for P out = 10 yr, t int = 10 8 yr would be a practically reasonable timescale to consider the fate of astronomical triples.
We emphasize that the above definition of the disruption is expected to correspond to the Lagrange stability, instead of the Lyapunov stability studied by MA01. The Lyapunov instability may be a sufficient condition to the Lagrange instability (for infinite future). In practice, however, Lagrange instability is difficult to be identified with a limited computational time and its numerical accuracy. Thus we decided to examine the corresponding disruption timescales up to 10 7 P out , and not to repeat the Lyapunov-type analysis in the present paper. We simply clarify that the two stabilities are different concepts, and their relation is not trivial in general. As we show below, Lyapunov stability, at least according to the criterion by MA01 that involves an empirical extrapolation of the inclination dependence, is neither a necessary nor sufficient condition for Lagrange stability up to 10 7 P out .

Previous Models of the Disruption Timescale and Lyapunov Stability Criteria
In order to illustrate the basic idea behind the present analysis, we show in Figure 1 an example of the distribution of the disruption timescale T d for the simulated orthogonal triples on e out -r p,out /a in plane. Here, e out is the eccentricity of outer orbit, and r p,out /a in is the ratio of the pericenter distance of outer orbit and the semimajor axis of inner orbit.
Following Paper I, we sequentially perform the simulations at each e out from lower to higher values of r p,out /a in . Once the previous two realizations at lower r p,out /a in do not disrupt before t int , we stop further simulations for one sequence, and define the boundary to save computational cost. The detail of this procedure is described in Paper I. The left and right panels in Figure 1 plot T d /P in and T d /P out , respectively. While the information context of the two panels is essentially the same, their comparison provides a useful insight into the purpose of the present paper.
For reference, we overlay the disruption timescale contours predicted from the random walk (RW) model by Mushkin & Katz (2020) where m 123 ≡ m 1 + m 2 + m 3 is the total mass of the system, and m 12 ≡ m 1 + m 2 is the total mass of inner binary.
Equations (1) and (2) serve as useful analytical estimates for the disruption timescale of Lagrange-unstable systems.
In Figure 1, we also include the dynamical stability criteria by MA01 and Vynatheya et al. (2022) as magenta solid and dashed lines, respectively. Vynatheya et al. (2022) defined the boundary from their simulations at 100P out by selecting the triples whose inner and outer semimajor axes do not change by more than 10% of the initial values. This is not exactly the same condition adopted by MA01, but it is expected to be similar from the viewpoint of Lyapunov stability.
According to them, those triples with r p,out /a in exceeding the critical thresholds below are Lyapunov stable, respectively:  (1) and (2)). Cross symbols indicate the triples that have not been broken until t int . Left and right plots adopt t int = 10 9 P in and t int = 10 7 P out , respectively. We note that the left plot is adapted from Paper I.  Figure 1 is adapted from the simulation results with t int = 10 9 P in in Paper I. Due to the chaotic nature of three-body dynamics, the distribution of T d exhibits relatively large scatters; see Figures 6, 7, 16, 17, and 18 in Paper I. Nevertheless, the values of T d /P in are basically determined by the value of r p,out /a in . This is consistent with the RW model prediction, Equation (1), which suggests that T d /P in depends on e out very weakly (µ e 1 out ). The right panel shows T d /P out from our new simulations with t int = 10 7 P out .
Since Equations (3) and (4) represent the thresholds for Lyapunov stability, there is no reason why their criteria are directly related to the disruption timescale characterizing Lagrange stability. Figure 1 indicates that their Lyapunov stability boundaries indeed correspond to a broad range of disruption timescales, 10 2 < T d /P out < 10 5 in our simulation runs. In other words, their criteria should not be interpreted as a black-and-white boundary of the disruption of triples. This seems to be a common misinterpretation in works studying triple stability (R. Mardling 2022, private communication). This is why we have repeatedly stressed the conceptual difference between Lagrange and Lyapunov stabilities.
The comparison of the two panels in Figure 1 motivates us to study more systematically the dependence of T d /P out on the mutual inclination i mut of the inner and outer orbits. We note again that Paper I adopted the normalized disruption timescales of T d /P in , while the present paper uses T d /P out . This is not essential and just a matter of definition, but the different normalization may be useful in providing a complementary view to the problem relative to Paper I. Furthermore, if the outer binary is highly eccentric, most of the energy transfer between the two orbits, which is responsible for triggering the instability, should occur once every P out when the outer body is at pericenter. Therefore, the scaling of T d /P out may be more relevant from the physical perspective.

Mutual Inclination Dependence of the Disruption Timescales
Figure 2 plots T d /P out from our simulations on the e out − r p,out /a in plane for different values of i mut . This generalizes Paper I that focused on three representative mutual inclinations (i mut = 0°, 90°, and 180°) alone, and quantitatively estimates the disruption timescales associated with Lagrange stability by systematically varying the initial inclinations i mut in an interval of 15°. We adopt the same initial parameter sets in Paper I and Figure 1, and fix the initial phases as (ω in , M in , ω out , M out ) = (180°, 30°, 0°, 45°). The estimated timescales are supposed to vary by 1 or 2 orders of magnitude when the initial phases are randomly assigned, in addition to the similar amount of scatters due to the intrinsic chaotic nature of the triple dynamics (Paper I). Thus, the distribution of T d /P out shown in Figure 2 should be understood to have a scatter of 1 or 2 orders of magnitude in general. The variation due to the different initial phases will be examined in Figure 3 below.
For reference, we plot the stability boundary proposed by Vynatheya et al. (2022) alone in Figure 2; as Figure 1 indicates, the MA01 boundary is roughly identical. The top three panels indicate that triples with < ( ) r a r a p,out in p,out in V are disrupted approximately within T d < 10 3 P out , while those with > ( ) r a r a p,out in p,out in V mostly survive more than 10 7 P out . In that sense, Lyapunov stable triples with i mut < 45°(nearly coplanar prograde) are Lagrange stable as well.
However, it is not the case in general. In particular, strong ZKL oscillations for i mut > 60°tend to destabilize the triple, and the fate of those triples can be revealed only by integrating them much longer than the corresponding quadrupole ZKL timescale (e.g., Antognini 2015). Indeed, the disruption timescales for triples with 60°< i mut < 150°monotonically increase as r p,out /a in without any discontinuous change around = ( ) r a r a p,out in p,out in V . This fact implies that the Lagrange instability is mainly triggered via the repeated ZKL oscillations over long timescales.
Furthermore, the significant suppression of the energy transfer between the inner and outer orbits for i mut ≈ 180°s trongly stabilize those triples in coplanar-retrograde orbits (Paper I) from the viewpoint of disruption timescales. Thus an accurate prediction of their fate requires a much longer time integration as well. Figure 2 clearly illustrates the impact of i mut on the disruption timescale of triples. In particular, it is worth noting that the Lagrange stability boundary appears to be highly asymmetric with respect to i mut = 90°. This is not expected from the quadrupole ZKL interaction mechanism alone, because it should be symmetric around 90°. The broken symmetry might arise due to semisecular effects that excite the eccentricity of the inner binary more efficiently than the simple quadrupole interaction mechanism (Grishin et al. 2018;Mangipudi et al. 2022), and thus increase the energy diffusion between inner and outer orbits (Paper I).

Dependence of Disruption Timescales on r p,out /a in , e out and i mut with Randomized Initial Phases
Due to the chaotic nature of the triple dynamics, the result of Figure 2 may change significantly by the difference of the initial phases (see Paper I). Therefore, we examine the sensitivity of the disruption timescales to the initial phases by computing sections of Figure 2 along the constant e out direction. To be specific, we fix M in = 0°and ω in = 0°, and change r p,out /a in by a step of 0.05, which is 4 times denser sampling than those adopted in Figure 2. Furthermore, we compute 10 different realizations for the same set of parameters, (e out , r p,out /a in , i mut ) using randomly selected values of ω out and M out . We emphasize that the tiny difference of other orbital parameters induces the scatters of T d /P out by a similar amount (see Paper I). Figure 3 shows the resulting T d /P out distributions for coplanar-prograde (top), orthogonal (middle), and coplanarretrograde (bottom) triples. Just for reference, we plot the expected stability boundary, Equation (4), in dashed lines. The top panel shows that T d /P out is a very sensitive function of r p,out /a in for low-inclined systems, even taking into account the scatters due to initial phase differences; T d /P out changes significantly around the transition region corresponding to the value of 10 4 -10 5 . The top panel have very few data points between 10 5 and 10 7 , while a majority of triples are simply Figure 2. The normalized disruption timescales T d /P out on e out − r p,out /a in plane for different values of i mut . The disruption timescales evaluated from simulations are indicated according to the side color scales. Magenta curves represent the dynamical stability criterion from Vynatheya et al. (2022;see Equation (4)). Cross symbols indicate the triples that have not been broken until t int = 10 7 P out . While the result for i mut = 30°is not shown here, it is very similar to a simple interpolation of those for i mut = 15°and 45°. located at the upper limit of integration time (T d /P out = 10 7 ). This result agrees qualitatively with that of Vynatheya et al. (2022) except for highly eccentric cases. For a highly eccentric case of e out = 0.9, the transition becomes less abrupt; T d /P out increases gradually with r p,out /a in up to 4.5, and then shows significant jump.
The bottom panel, for coplanar-retrograde triples, shows a similar behavior, although the systems are significantly more stable than prograde cases, and the transition looks less abrupt.
In contrast, for orthogonal triples (middle panel); T d /P out increases gradually with r p,out /a in until r p,out /a in ≈ 3 and 4 for e out = 0.1 and 0.5, respectively, and then exhibits significant increase only when r p,out /a in exceeds those values. Even after exceeding those values, the transition is still not so abrupt for orthogonal triples, unlike the case for coplanar-prograde triples. For e out = 0.9, T d /P out continues to increase gradually with r p,out /a in up to r p,out /a in = 6 where T d /P out reaches our upper limit of the integration time (t int = 10 7 P out ). Those trends are statistically robust against the different choices of the initial phases as illustrated by the scatters of the data points.

Summary and Conclusion
We have computed the distribution of disruption timescales for hierarchical triples with an equal-mass and initially circular inner binary, using a series of N-body simulations. We extend Paper I and systematically examine the dependence on the mutual inclination between the inner and outer orbits in particular. By integrating the triple systems much longer than the previous studies, we are able to reveal the fate of those triples more reliably than before. Our main findings are summarized as follows.
1. In studying the dynamical stability of hierarchical triple systems, it is important to distinguish between Lyapunov and Lagrange stabilities. A widely used stability criterion by Mardling & Aarseth (1999, 2001 and Mardling (2008) is derived from the local divergence of the orbits, and thus corresponds to the former. Paper I and the present paper compute the disruption timescale of triple that is related to the Lagrange stability. Those concepts are complementary but very different, which should be kept in mind when applying those results for specific purposes. 2. Dynamical disruption timescales of triples are very sensitive to the mutual inclination i mut of the inner and outer orbits. They change drastically around the Lyapunov stability boundary if i mut < 45°. Thus, the Lyapunov stability condition also guarantees the long-term Lagrange stability for moderately inclined triples. In contrast, for triples of significantly inclined inner and outer orbits (60°< i mut < 150°), the disruption timescales vary smoothly against r p,out /a in due to the ZKL oscillations. The difference of the initial phases adds 1 or 2 orders of magnitudes scatters to the disruption timescales, but does not change the behavior in a statistical sense. 3. The Lagrange stability of triples does not change monotonically with the mutual inclination. For triples with moderate inclinations (i mut < 45°), the stability is not so sensitive to i mut . By increasing i mut up to ∼100°, those triples become less stable due to the ZKL oscillations.
Retrograde triples with i mut > 120°become stabilized again as i mut increases, and the coplanar-retrograde triple (i mut = 180°) is the most stable configuration. Grishin et al. (2017) found the similar behavior about the mutual inclination dependence in the Hill stability, and pointed out that the ZKL oscillations affect the stability. Interestingly, the above behavior is not symmetric around i mut = 90°, in contrast to a naive expectation from the quadrupole ZKL oscillations, indicating the importance of quasi-secular effects (e.g., Grishin et al. 2018;Mangipudi et al. 2022) beyond the simple ZKL quadrupole interaction. Figure 3. Behavior of the stable and unstable transition for triples with i mut = 0°, 90°, and 180°. In each panel, the different colored symbols plot T d / P out against r p,out /a in for e out = 0.1 (red), 0.5 (blue), and 0.9 (black). For each set of (e out , r p,out /a in , i mut ), we run 10 realizations with different initial phases (ω out , M out ) randomly sampled from uniform distribution. The value of N in the panels denotes the total number of the corresponding runs (data points).