New insight of time-transformed symplectic integrator I: hybrid methods for hierarchical triples

Mikkola & Aarseth (2002) developed a time-transformed LeapFrog method, which is more general than the LogH method, although it may sacrifice symplecticity. This method defines two time-transformation functions for the drift and kick steps of the LeapFrog approach. During kick steps, the function is defined as:

where $\Omega(\mathbf{r})$ is an abitrary function of $\mathbf{r}$.

In drift steps, the function is defined as¹¹1Here, $u$ represents $W$ as defined in Equation 12 of Mikkola & Aarseth (2002)

where $u$ is a new integrated variable given by

with $\mathbf{v}$ representing the velocity vectors. Since $u$ depends on $\Omega(\mathbf{r})$, it is calculated during the kick steps and used in the drift steps.

For conservative potential, we can set $\Omega(\mathbf{r})$ as $-U(\mathbf{r})$. In this case, the method (MA2002) approximates the LogH method. Here, $g_{\mathrm{k}}^{\mathrm{MA}}$ is equivalent to the function $g_{\mathrm{k}}$ (Equation 10) of the LogH method. The corrresponding $u$ can be expressed in integral form²²2The same definition of $u$ in Equation 40 of Wang et al. (2020a) is missing a negative sign and $\mathrm{d}t$.:

By setting the initial value of $u$ to the initial value of $-U(\mathbf{r})$, we have $u\equiv-U(\mathbf{r})$. Due to energy conservation, $T(\mathbf{p})+p_{\mathrm{t}}=-U(\mathbf{r})$, $u$ can be viewed as $T(\mathbf{p})+p_{\mathrm{t}}$. Thus, we can find that

which is equivalent to $g_{\mathrm{d}}$ function (Equation 10) of the LogH method. However, in numerical simulations, $T(\mathbf{p})+p_{\mathrm{t}}\approx-U(\mathbf{r})$, making $u$ an approximation of $T(\mathbf{p})+p_{\mathrm{t}}$.

Therefore, the MA2002 method can be viewed as the LogH method with $g_{\mathrm{d}}$ approximated by $g_{\mathrm{d}}^{\mathrm{MA}}$. The drift and kick steps can be described using the same Equations 8 and 9, with $g_{\mathrm{k}}$ defined as in Equation 10 and $g_{\mathrm{d}}$ replaced by Equation 14. The value of $u$ can be calculated during kick steps as:

with $\mathbf{v}|_{1}$ calculated in Equation 9. This formula approximates $u$ assuming that $\partial U(\mathbf{r})/\partial\mathbf{r}$ is constant during one step.

The original LogH method is $\Delta s$ symmetric or time symmetric because reversing $\Delta s$ for backward integration yields the same formulae of integrated quantities (Equations 8 and 9). The MA2002 method retains this feature, including the additional integration of $u$ in Equation 18. Thus, the MA2002 method also achieves good long-term behavior with respect to the cumulative energy error.

Due to approximation errors, the computation of $g_{\mathrm{d}}$ in the MA2002 method differs from the exact value in Equation 10. Consequently, the MA2002 method is less accurate than the original LogH method, especially for the Kepler orbits. The error differences between the two methods are detailed in the Appendix A. For isolated binaries, the original LogH method achieves better accuracy, with the $\Gamma$ error mainly caused by round-off errors. In contrast, the MA2002 solution closely follows the Kepler trajectory, but exhibits some secular error patterns. However, both methods show similar accuracy for triple systems. Thus, the MA2002 method serves as a good approximation to the LogH method and adequately demonstrates its limitations in evolving hierarchical triples in the following sections.

In addition, due to flexible changes in $\Omega(\mathbf{r})$ and the definition of $u$ in the MA2002 method, we use the MA2002 method to develop a new algorithm in this work. To ensure consistent comparisons among methods, the term “LogH method” will refer to the MA2002 method in the following content.

The issue discussed can be avoided in a hybrid scenario. In $N$-body codes like nbody7, sdar and petar, the LogH method is combined with other integrators, such as 4th-order Hermite with block time steps, to form a hybrid method referred to as ”H4”. In this case, the LogH method is used for the inner binary, while another integrator manages the outer orbit. The influence of the third body on the inner binary is treated as a perturbative force during the kick step in the LogH method, allowing the ability to accurately approximate Keplerian trajectories to be preserved.

However, this hybrid method has a disadvantage. The primary bottleneck is that the time steps $\Delta t_{\mathrm{B,k}}$ in the LogH method are quantities to integrate, making it difficult to design $\Delta s$ to ensure integration completes at the required time $\Delta t_{\mathrm{H}}$ from the other integrator. As a result, additional iterative steps are necessary for time synchronization. If $\Delta t_{\mathrm{H}}$ is comparable to $\Delta t_{\mathrm{B,k}}$, the hybrid method becomes time-consuming and accumulates timing errors from frequent synchronizations. If $\Delta t_{\mathrm{H}}$ exceeds the period of the inner binary, the interaction between the inner binary and the outer body occurs less than once per inner orbit, resulting in inaccurate evolution.

We introduce a new hybrid method ”BlogH”, which avoids time synchronization. In this approach, all objects in the triple are integrated using Equations 8 and 9, but the definitions of $g_{\mathrm{d}}$ and $g_{\mathrm{k}}$ are modified. In the BlogH method, only the potential energy of the inner binary,

is included in $g_{\mathrm{k}}$ as follows:

where new symbol $g_{\mathrm{k,b}}$ distinguishes it from $g_{\mathrm{k}}$ in the LogH method.

Since $g_{\mathrm{k}}$ and $g_{\mathrm{d}}$ must represent the equivalent quantities to ensure consistent time steps for kick and drift, and considering total energy conservation, we have

the corresponding $g_{\mathrm{d,b}}$ in the BlogH method is given by:

which depends on coordinates, making it unsuitable for the explicit method. The solution is to adopt a similar approach to that in Equation 16 to approximate $g_{\mathrm{d,b}}$ by defining $u_{\mathrm{b}}$ as:

where $\mathbf{v}_{\mathrm{b}}$ includes only the inner binary components. This formula can be evaluated at kick steps following Equation 18. Then,

which can be used in the explicit method.

This scenario enables the inner binary to be integrated with the LogH method, while the outer orbit is integrated with the sDKD method, following the same time-step sequences $\Delta t_{\mathrm{B,k}}$ as the inner binary. This setup allows the third body to provide perturbations at the kick step without affecting the time transformation. The method retains the advantages of accurately approximating Keplerian trajectories without the need for time synchronization. The disadvantage is that the method is not symplectic due to the approximations in the MA2002 method and the unequal time steps of the sDKD loop for outer orbits. However, compared to the H4 method, BlogH is time symmetric. Figure 3 illustrates the differences in treatment in time steps among the three methods: LogH, H4 and BlogH.

The BlogH method can also be applied to hierarchical multiple systems that contain one innermost binary. For example, consider a hierarchical quadruple system with an inner binary with components 1 and 2, and two outer bodies, 3 and 4, which do not form a binary. The dynamics of the four bodies are integrated using Equations 8 and 9. However, only the innermost binary (components 1 and 2) is used to calculate $g_{\mathrm{k,b}}$ and $g_{\mathrm{d,b}}$ through Equations 21, 24 and 25.

Implementing the BlogH method in a LogH $N$-body code is straightforward. The only necessary modification is to update the functions $g_{\mathrm{d}}$ and $g_{\mathrm{k}}$, while the overall framework of the integrator remains unchanged. We have incorporated the BlogH method into the sdar code, alongside the existing LogH and H4 methods.

To show that the hybrid methods outperform the LogH method, we redo simulations of the triple system including a weakly perturbed binary by using the LogH, H4 and BLogH meothods with varying $\Delta s$. We evolve the triple system to 500 periods of the inner binary ($P_{\mathrm{i}}$) to assess the secular evolution of energy and orbital elements.

In the triple system, since the outer orbit is circular, the time step of the Hermite part ($\Delta t_{\mathrm{H}}$) can remain constant in the H4 method. According to the block time step scheme, we select $\Delta t_{\mathrm{H}}=0.5^{14}\approx 0.307P_{\mathrm{i}}$ for the maximum $\Delta s$. When $\Delta s$ is reduced by a factor of 4, $\Delta t_{\mathrm{H}}$ is halved to maintain consistent accuracy between the second-order LogH method and the fourth-order Hermite method.

Comparing the evolution of $\epsilon$ reveals a significant advantage for the H4 and BlogH methods. As shown in Figure 4, both methods demonstrate much better energy conservation than the LogH method with the same $\Delta s$. In particular, $\epsilon$ oscillates significantly in the LogH results, unlike in the other two methods. Compared to Figure 2, the H4 and BlogH methods show a behavior similar to the LogH method for isolated binaries, suggesting that hybrid methods can effectively approximate the Kepler solution for the inner binaries.

A counterintuitive behavior of the BlogH methods is that a smaller $\Delta s$ leads to a larger $\epsilon$. This occurs because the influence from round-off errors becomes significant and exceeds the approximation error. To investigate this, we performed reference simulations using a high-precision BlogH method (BlogH-HP) which applies quadruple precision with 30 digits. As shown in the last panel of Figure 4, the same $\Delta s$ results in a lower $\epsilon$, and the relationship between $\Delta s$ and $\epsilon$ aligns with the 2nd-order integration method.

For both BlogH and BlogH-HP results, secular patterns emerge in the evolution of $\epsilon$. To investigate their origin, we further explore how the evolution of orbital parameters from different methods influences the behavior of $\epsilon$. In the triple system with a weakly perturbed binary, changes in orbital parameters are minimal. To highlight the differences among the methods in the plot, we define the relative evolution of orbital parameters as

where $o$ represents orbital parameters such as $a_{\mathrm{i}}$, $a_{\mathrm{o}}$, $e_{\mathrm{i}}$ and $e_{\mathrm{o}}$, and $o[0]$ denotes the initial value at $t=0$.

The results are presented in Figure 5. In addition to the BlogH-HP models, we introduce the BlogH-Y6-HP models, which use both quadruple precision and the Y6 method to enhance integration accuracy. The BlogH-Y6-HP-S256 model is the most accurate simulation and serves as a reference. By comparing the evolution of $\Delta a_{\mathrm{i}}$ and $\Delta e_{\mathrm{i}}$ across different methods, we observe that the significant oscillation of $\epsilon$ in the LogH models arises from errors in the orbital elements of the inner binary. The LogH models exhibit substantial discrepancies compared to the reference model, with smaller $\Delta s$ resulting in reduced discrepancies. In contrast, both the H4 and BlogH methods demonstrate comparable inner orbital evolution, aligned closely with the reference model and showing a weak dependence on $\Delta s$.

The evolution of outer orbital elements shows that both LogH and BlogH are similar but significantly larger than the reference model. This explains the secular change in $\epsilon$ for the BlogH method (Figure 4), which arises from errors in the outer orbital evolution. Furthermore, for the BlogH-S256 model, the $\epsilon$ curve in Figure 4 roughly follows the lower boundary of the LogH curve, suggesting that the LogH $\epsilon$ curve combines the oscillations of inner binary errors with the secular component from outer binary errors.

In contrast, the H4 method shows improved outer orbital evolution, with smaller $\Delta a_{\mathrm{o}}$ and $\Delta e_{\mathrm{o}}$ closely aligned with the reference model. This indicates that the fourth-order H4 method is more accurate for secular evolution of the outer orbit, even though it is not symplectic.

There is no clear trend regarding how the secular patterns of outer orbital evolution appeared in LogH and BlogH results depend on $\Delta s$. One possible explanation is that different $\Delta s$ values lead to varying time phase errors in the inner orbits, resulting in different interactions between the inner and outer orbits and, consequently, different secular patterns.

Although the BlogH method is less accurate in secular evolution compared to the H4 method, it incurs significantly lower computing costs because of the absence of time synchronization. We compare the total number of integration steps, $n_{\mathrm{s}}$, reflecting these costs for the H4 and BlogH methods in Table 1. The total steps for the H4 method, $n_{\mathrm{s}}(H4,tot)$, are 2-5 times greater that those for the BlogH method, depending on $\Delta s$, The steps for time synchronization, $n_{\mathrm{s}}$(H4,syn), account for approximately $1/3-1/2$ of the total steps. This suggests that time synchronization significantly affects the performance of the H4 method.

The values of $n_{\mathrm{s}}$ can indicate the performance of the algorithm; however, the actual wall clock time ($t_{\mathrm{w}}$) per inner orbit of the simulation can be influenced by various aspects of implementation. Therefore, we compare the maximum $\epsilon$ within the evolution of $5~{}P_{\mathrm{i}}$ versus $t_{\mathrm{w}}$ for the three methods–LogH, H4 and BlogH–across different $\Delta s$ values, both with and without the Y6 method. The results are shown in Figure 6.

The LogH results show an increasing trend in $t_{\mathrm{w}}$ as both $\Delta s$ and $\epsilon$ decrease. However, the H4 and BlogH methods have $\epsilon$ values that are already dominated by round-off errors, so smaller $\Delta s$ does not significantly reduce $\epsilon$ and may even lead to an increase. Therefore, it is preferable to use a larger $\Delta s$ for the H4 and BlogH methods, such as S8, for better performance.

When only the second-order method is employed, both H4 and BlogH exhibit a few orders of magnitude less $t_{\mathrm{w}}$ compared to the LogH method. However, when the Y6 method is applied, LogH-Y6 becomes comparable to the H4 method at $\epsilon\approx 10^{-11}$. In contrast, the BlogH method still demonstrates at least one order of magnitude less $t_{\mathrm{w}}$, indicating its significant advantage in computational efficiency.

Please note that the absolute values of $t_{\mathrm{w}}$ for each method can vary depending on the computing hardware and the implementation of the algorithm. The H4 method implemented in the code sdar is significantly more complex than the LogH or BlogH methods; therefore, although the values of $n_{\mathrm{s}}$ for the H4 method may appear better than those for the LogH method at $\epsilon\approx 10^{-11}$, the actual $t_{\mathrm{w}}$ remains comparable.

The Hybrid and BLogH methods are accurate not only for triples including a weakly perturbed binary but also significantly enhance accuracy in solving triple systems experiencing Kozai-Lidov oscillations with a highly eccentric inner binary.

When a triple system meets the orbtial criterion for Kozai-Lidov oscillation, the inner eccentricity $e_{\mathrm{i}}$ and the inclination angle $\theta$ between the inner and outer orbits do oscillate on a timescale given by (e.g. Antognini, 2015):

In the following example, we set an initial condition for a Kozai-Lidov triple shown in Figure 7. The outer body has $m_{3}=10$, significantly more massive than the inner binary with $m_{1}=0.009$ and $m_{2}=0.001$. The initial inclination $\theta$ is $90^{\circ}$ and the corresponding $t_{\mathrm{KL}}$ is approximately $1138$. We can scale the triple to different astrophysical systems. For instance, with a mass unit of $100~{}M_{\odot}$, the triple represents a stellar-mass binary of $0.1$ and $0.9~{}M_{\odot}$ orbiting an intermediate-mass black hole of $1000~{}M_{\odot}$. Alternatively, with a mass unit of $0.1~{}M_{\odot}$, the triple represents a solar system with binary Jupyter-mass planets.

Using the three methods, we integrate the evolution of the triple for about 2 $t_{\mathrm{KL}}$. The orbital parameters and $\epsilon$ evolution are shown in Figure 8. We compare the results of all the methods with a step size of $\Delta s=\Delta s_{0}/32$. The LogH-Y6-S32 model exhibits a markedly different evolution of $e_{\mathrm{i}}$ compared to the H4-Y6-S32 and BlogH-Y6-S32 models, indicating an unphysical result of the LogH method.

To verify this, we add two simulations using the LogH method with much smaller step sizes: S512 and S4096. These models show a consistent evolution of $e_{\mathrm{i}}$ with H4-Y6-S32 and BlogH-Y6-S32, confirming that LogH-Y6-S32 presents an unphysical result, despite its lower $\epsilon$. Moreover, even with smaller steps, the LogH models display jumps in $a_{\mathrm{i}}$, $e_{\mathrm{o}}$ and $\epsilon$ when $e_{\mathrm{i}}$ approaches unity (approximately $0.9999999$), This behavior is less pronounced in the BlogH-Y6-S32 model and absent in the H4-Y6-S32 method, both at the same maximum $e_{\mathrm{i}}$. This suggests that the LogH method is significantly less accurate than the other two when dealing with high eccentric Kozai-Lidov triples.

We further compare the step counts of the methods in Table 2. The H4-Y6-S32 model has approximately twice as many as the BlogH-Y6-S32 model, suggesting that time synchronization affects the efficiency of the H4 method, consistent with the findings for triples including a weakly perturbed binary in Section 4.3. With the same $\Delta s$, the LogH-Y6-S32 model requires about 7.4 times more steps than the BlogH-Y6-S32 model. This indicates that the relationship between $\Delta s$ and $\Delta E$ in Equation 11 for the Kepler orbit is also disrupted in the LogH-Y6-S32 model. The LogH-Y6-S8192 model, which exhibits similar accuracy in orbital evolution to the BlogH-Y6-S32 model, needs around 256 times more steps per orbit. Therefore, the BlogH method demonstrates a significant performance advantage.

Two major issues with the low accuracy of the LogH method are the unphysical evolution of $e_{\mathrm{i}}$ and the significant jump in $\epsilon$ when $e_{\mathrm{i}}$ is close to 1. To investigate how these behaviors depend on $\Delta s$, we perform a series of models using the LogH method with varying $\Delta s$. The evolutions of $e_{\mathrm{i}}$ and $\epsilon$ are compared in Figure 9.

The S16, S32 and S64 models show significantly different secular evolutions of $e_{\mathrm{i}}$. In these simulations, the $n_{\mathrm{s}}$ per $P_{\mathrm{i}}$ are 135, 269 and 538, respectively. Despite using the high-accuracy 6th-order Yoshida method, $e_{\mathrm{i}}$ still evolves unphysically. In addition, repeated peaks appear in the evolution of $\epsilon$ for the S32 and S64 models. These peaks occur at the periapsis of the inner binary, as illustrated by the LogH results for the triple, including a weakly perturbed binary (Figure 2). This suggests that the low accuracy at periapsis from the LogH method introduces an artifical perturbation that repeatedly impacts $e_{\mathrm{i}}$ and leads to an unphysical secular evolution.

Despite the unphysical evolution of $e_{\mathrm{i}}$, $\epsilon$ can still recover to a low value, aligned with the symplectic method. The maximum $\epsilon$ reaches approximately $10^{-8}$, corresponding to a relative energy error of $10^{-6}$. Although this relative error may be acceptable for $N$-body simulations of star clusters, it is insufficient for this triple system. We also investigated whether the artificial perturbation originates from backward steps in the 6th-order Yoshida method and found no correlation.

As $\Delta s$ increases from S128 to S8192, the evolution trend of $e_{\mathrm{i}}$ converges with the H4 and BlogH results, and jumps in $\epsilon$ begin to appear. While the jump value decreases, the average value of $\epsilon$ increases as $\Delta s$ decreases. This suggests that a smaller $\Delta s$ can help minimize a large jump but increases cumulative round-off errors. This inverse relationship between $\Delta s$ and the average $\epsilon$ also occurs in the simulations for the triple including a weakly perturbed binary, as illustrated in Figure 4.

An important feature is that the S32 and S64 models show a similar average $\epsilon$. Because $e_{\mathrm{i}}$ is far from 1, there is no jump in $\epsilon$. In contrast, for the S128 models, the jump occurs at $e_{\mathrm{i}}\approx 1$, resulting in a final $\epsilon$ that is significantly larger than those of the S32 and S64 models. If we rely solely on the final $\epsilon$ to assess integration accuracy without considering the secular evolution of $e_{\mathrm{i}}$, we may mistakenly conclude that the S32 and S64 models are more accurate.

We now compare the maximum $\epsilon$ and $t_{\mathrm{w}}$ for the LogH-Y6 models with different $\Delta s$ values, as well as for the BlogH-Y6-S32 and H4-Y6-S32 models. The results are presented in Figure 10. For the LogH-Y6 method, below S256, the evolution of $e_{\mathrm{i}}$ is unphysical; consequently, $\epsilon$ appears low and increases with $\Delta s$, along with $t_{\mathrm{w}}$, except for S16. Above S512, energy jumps when $e_{\mathrm{i}}$ approaches 1 dominate the maximum $\epsilon$, leading to a more reasonable trend where increasing $\Delta s$ leads to a lower maximum $\epsilon$. Comparing this with Figure 6, we observe that the H4 and BlogH methods demonstrate greater efficiency at equivalent $\epsilon$ values when solving the Kozai-Lidov triples with high eccentricities, and the H4 method being notably more efficient than the LogH method.

Note that for the LogH method, the maximum $\epsilon$ is primarily attributed to the energy jumps. It is also possible to reduce the step sizes when $e_{\mathrm{i}}$ is close to 1 to avoid significant energy jumps, although this breaks the symplecity of the (original) LogH method. Consequently, with larger $\Delta s$ values, such as S512, it becomes feasible to obtain a reasonable physical result with much lower maximum $\epsilon$.

The absence of time synchronization offers advantages to the BlogH method, but also limits its applicability. When two binaries exist in a system, such as a binary-binary quadruple, the LogH method must be applied to each binary with different time transformation functions. In this case, the same $\Delta s$ results in different $\Delta t$ values for the two binaries, making time synchronization unavoidable and preventing the use of the BlogH method. In future studies, we will explore the extension of the BlogH approach to systems with multiple binaries.

To address the low accuracy of the second order sDKD method in the LogH approach (for both original and MA2002 methods), Mikkola & Tanikawa (1999) applied the Bulirsch-Stoer extropolation integrator to create the LogH-BS method. This method uses polynomial extropolation on a sequence of LogH results obtained through the sDKD loop with varying step sizes of $\Delta s$, predicting an accurate result as the step size approaches zero.

The Bulirsch-Stoer method is efficient when paired with a classic integrator. However, the LogH-BS method is less efficient than the pure LogH method for binary systems. The LogH method can precisely follow the trajectory of a Kepler orbit with only a time-phase error, whereas the extropolated results from the LogH-BS method do not guarantee exact alignment with the trajectory. To achieve the same accuracy as the LogH method, the LogH-BS method requires more integration steps.

The LogH-BS method is more beneficial for systems with $N>2$. However, as illustrated in Figure 9, for Kozai-Lidov triples, a low $\epsilon$ does not ensure accurate secular evolution of $e_{\mathrm{i}}$. Therefore, only with a strict criterion of $\epsilon<10^{-11}$, can the LogH-BS method potentially provide a physical evolution of $e_{\mathrm{i}}$ while avoiding large errors at the pariapsis of the inner binary. In contrast, hybrid methods such as the H4 or BlogH method are much more efficient.

The hybrid methods can be combined with the Bulirsch-Stoer method for faster convergence of hierarchical triples. Further investigation of this possibility is needed in future work.

For all three methods, the time phase has errors due to the approximation in time integration, which manifests itself as an artificial phase shift when a large $\Delta s$ is used. This error also affects time synchronization in the H4 method. Depending on the systems being studied, this error can be significant or negligible. In planetary systems where orbital resonances are important, such phase errors may impact long-term evolution and should be carefully investigated. However, for multiple systems in star clusters, which is the focus of this work, the error in the time phase is less critical than the errors in energy and angular momentum. This is because studies of star clusters primarily concentrate on the statistical evolution of semi-major axes and eccentricities of binaries, which are essential for estimating merger rates and orbital properties of events like gravitational waves. Thus, the orbital (time) phase is not important. Consequently, orbital average methods that may produce incorrect time phases, such as the slowdown method (Mikkola & Aarseth, 1996; Wang et al., 2020b), are frequently employed in star cluster simulations to enhance efficiency, including codes like nbody6, petar and birfost. More research is needed to explore how orbital average methods can be integrated with the hybrid method.